0
Research Papers

A Geometrically Exact Contact Model for Polytopes in Multirigid-Body Simulation OPEN ACCESS

[+] Author and Article Information
Jedediyah Williams

CS Robotics Lab,
Department of Computer Science,
Rensselaer Polytechnic Institute,
Troy, NY 12180
e-mail: jedediyah@gmail.com

Ying Lu

CS Robotics Lab,
Department of Computer Science,
Rensselaer Polytechnic Institute,
Troy, NY 12180
e-mail: rosebudflyaway@gmail.com

J. C. Trinkle

Director
CS Robotics Lab,
Department of Computer Science,
Rensselaer Polytechnic Institute,
Troy, NY 12180
e-mail: trink@cs.rpi.edu

Contributed by the Design Engineering Division of ASME for publication in the JOURNAL OF COMPUTATIONAL AND NONLINEAR DYNAMICS. Manuscript received January 3, 2015; final manuscript received July 25, 2016; published online December 2, 2016. Assoc. Editor: Andreas Mueller.

J. Comput. Nonlinear Dynam 12(2), 021001 (Dec 02, 2016) (14 pages) Paper No: CND-15-1003; doi: 10.1115/1.4034390 History: Received January 03, 2015; Revised July 25, 2016

We present a formulation of nonpenetration constraint between pairs of polytopes which accounts for all possible combinations of active contact between geometric features. This is the first formulation that exactly models the body geometries near points of potential contact, preventing interpenetration while not overconstraining body motions. Unlike many popular methods, ours does not wait for penetrations to occur as a way to identify which contact constraints to enforce. Nor do we overconstrain by representing the free space between pairs of bodies as convex, when it is in fact nonconvex. Instead, each contact constraint incorporates all feasible potential contacts in a way that represents the true geometry of the bodies. This ensures penetration-free, physically correct configurations at the end of each time step while allowing bodies to accurately traverse the free space surrounding other bodies. The new formulation improves accuracy, dramatically reduces the need for ad hoc corrections of constraint violations, and avoids many of the inevitable instabilities consequent of other contact models. Although the dynamics problem at each time step is larger, the inherent stability of our method means that much larger time steps can be taken without loss of physical fidelity. As will be seen, the results obtained with our method demonstrate the effective elimination of interpenetration, and as a result, correction-induced instabilities, in multibody simulations.

There is a wide spectrum of applications for simulation, but an understandable bias toward real-time interactivity has bred a culture that leaves concerns of physical accuracy in a dark corner. As a result, applications which require some level of physical fidelity from simulation are left with a relatively underdeveloped set of tools. The major focus of the work presented herein is on improving the simulation accuracy, and consequently stability, without significantly increasing the computational cost.

The tradeoff of speed versus accuracy in simulation is a classic one, and many applications simply do not gain from improved physical accuracy. The gaming community is naturally concerned with methods that are fast and reliable, but not necessarily accurate [13], inspiring the term “game physics.” Many games, in fact, intentionally incorporate nonrealistic physics into their play. The graphics community, which is constantly pushing the boundaries of realistic-looking visual effects, does not offer many improvements to physically accurate simulation. Although there are many examples in a long history of “physics-based” methods in use in animation [48], there is really no need for such methods to do more than approximate some convenient alternative version of reality. Indeed, the goal in many graphics applications is to present the viewer with something that “feels” correct.

Particularly in the field of robotics, physical fidelity is a pressing concern and is especially challenging when simulating systems experiencing intermittent frictional contact, which is the case when robots work in unstructured, cluttered environments. Robotics experiments tend to be exceedingly costly and time consuming, so that the ability to perform realistic experiments in simulation would be a great benefit. Although work has been done regarding validation of simulation [912], the problem of measuring and quantifying the physical accuracy of simulation is an open issue that is challenging to solve. Comparison of a particular simulator's results with those obtained from a commercial simulator is meaningless without proof of the commercial simulator's accuracy. Even validation of a simulator with a physical experiment can be misinterpreted when it takes the tuning of several parameters to match a single or small number of experiments. Such an attempt at validation says little, if anything, of the simulator's ability to predict a different experiment.

Simulation can be a powerfully useful tool if shown to be accurate. A simulation requires accuracy in order to make any claims as to its ability to predict behavior in the physical world, yet the growing number of variants on old simulation methods all have at least this in common: they sacrifice physical fidelity at multiple stages in order to achieve computational speed-up. The ideal simulation tool would be both accurate and fast, but when we sacrifice accuracy, we are necessarily sacrificing the usefulness of that tool.

The polytopal exact geometry (PEG) contact model presented herein is a step toward improving the physical fidelity of simulations of polytopal rigid bodies in intermittent contact. In most popular simulation tools, in order to make it easy to formulate the nonpenetration constraints, the shapes of the bodies are effectively altered. Specifically, facets of polytopes are represented as hyperplanes (e.g., faces of polyhedra in 3D are represented as planes). As will be illustrated below, the price for this convenience is very high, as it leads to nonphysical simulation artifacts, such as huge constraint forces that eject bodies from the scene, or worse, cause the simulation's solvers to fail completely. By contrast, PEG exactly represents the boundaries of polytopal bodies, including their finite extent, and thus avoids these problems. To achieve this, the PEG contact model is more complicated than the standard contact model. With the same time-step sizes, simulation with PEG takes approximately twice as long as simulation with the standard contact model. However, its improved accuracy avoids nonphysical artifacts, including penetration-induced instabilities that can crash a simulation. In our opinion, this speed reduction is a penalty worth paying.

In the remainder of this section, we review multibody dynamics and some of the approaches used in the simulation. In Sec. 2, we detail the deficiencies of standard contact modeling practices and begin to derive the PEG contact model by introducing new approaches to cases that standard contact models fail to accurately represent. In Sec. 3, we define the general formulations of three classes of nonpenetration constraints and provide abstractions of them in order to model contact in a clean way. In Sec. 4, we precisely define contact constraint applicability and feasibility, and establish geometric tests to determine which contacts to include in PEG. In Sec. 5, we utilize the abstractions from Sec. 4 in order to define four fundamental contact constraints for spatial problems, which reduce to two when specialized to the planar case. We then demonstrate how to formulate a time-stepping subproblem that incorporates PEG. In Sec. 6, we compare the performance of multibody simulations using a popular time-stepping method with the standard contact model against the same time stepper with the PEG contact model. Finally, Sec. 7 discusses our results and possibilities for future work.

Previous Work.

The two main classes of simulation methods for multibody systems are known as event-driven and time-stepping. In event-driven methods, the simulation is monitored for impact and separation events as well as changes between sticking and slipping (see Ref. [13]). For systems with small numbers of contacts, this method is preferred, because it uses an adaptive step size to accurately pinpoint and process the contact state transitions. Its downfall is that in systems with many bodies and contact events, the size of the time step becomes so small that computation times become impractically large. On the other hand, this is precisely the circumstance in which time-stepping methods have their greatest value. Time-stepping methods use large time steps that can advance the simulation across a number of contact events in one step, such that at the end of the time step, Newton's laws are satisfied and bodies do not overlap [14].

Time-stepping simulation methods utilizing the complementarity problem (defined formally below) were introduced by Moreau and Panagiotopoulos [15]. In these methods, the data needed to advance the simulation by one time step are obtained as the solution of a complementarity problem, the so-called time-stepping subproblem. Stewart and Trinkle introduced a modification of this problem that incorporated stabilization of nonpenetration constraints [16], and they provided a proof that body trajectories approximated by the solutions to a sequence of these subproblems converged to a solution of the instantaneous-time formulation as the step size goes to zero [17]. Anitescu and Potra [18] included bilateral constraints and showed that the subproblem with linearized friction cones is always solvable.

Physics engines such as DART [19] use semi-implicit, velocity-based time-stepping subproblems, formulated as linear complementarity problems with nonpenetration constraint stabilization. However, because it uses standard contact modeling, it can suffer from inaccuracy and instability as discussed above. Chrono::Engine [20] is built on top of the Bullet collision detection engine [21], which reports penetration depth and gap distance, implying that penetration is not avoided. Across the gamut of available physics engines, simulations with bodies modeled as spheres or unions of spheres are popular for benchmarking, especially when there are large numbers of contacts. This is because when simulating polytopes in contact using the standard contact model, simulations can fail due to nonphysical, overconstraint imposed on the polygons. These overconstrained situations, henceforth referred to as standard-model traps (SM traps), can cause bodies to obtain huge nonphysical velocities or result in complete failure of the simulation.

Background.

In this section, we will introduce an instantaneous-time model of polyhedral bodies in intermittent contact, and then discretize it to obtain the discrete-time counterpart. To simulate the motions of the bodies, the discrete-time model (i.e., the time-stepping subproblem) must be formulated and solved repeatedly, each time advancing the simulation forward in time by the chosen time step. Note that the development herein excludes contact friction in order to maintain a sharp focus on the new contact model, which depends only on body geometry and not friction. Friction can, however, be easily incorporated by augmenting the time-stepping subproblem with one's favorite friction model, several of which are detailed in Refs. [22,23].

Before formulating PEG and a time-stepping subproblem that employs it, it is important to introduce the complementarity problem [24].

Complementarity.

The complementarity problem offers a convenient way to model contact because it naturally represents disjunctive conditions. The two key cases in multibody dynamics are: two bodies are in contact OR not in contact and an existing contact is sticking OR sliding. There is a vast body of research on complementarity problems and methods for their solution [25,26].

Let n denote n-dimensional Euclidean space. Given a mapping f(z):nn, the complementarity problem is to find the unknown vector zn that satisfies the three conditions Display Formula

(1)z0f(z)0zTf(z)=0

where T is the matrix transpose operator. For the sake of brevity, we may write all three conditions as Display Formula

(2)0f(z)z0

where represents the orthogonality of condition 3 in Eq. (1). It is common to refer to all three of these conditions jointly as a complementarity condition.

There are many general-purpose solvers available for complementarity problems [27] and in particular for linear complementarity problems (LCPs) [24,28]. In an LCP, f(z) is a linear function of z, i.e., f(z)=Az+b, where A is a known constant square matrix, and b is a known constant vector of compatible size. Much theory for LCPs has been developed and used to design reliable solvers. It is known, for example, that if A is positive definite, then the LCP is equivalent to a convex quadratic program, and therefore has a unique solution and can be solved efficiently. If A is positive semidefinite, then a solution can be found efficiently also, but some unknowns may not have unique solutions. In multibody dynamics, this case arises only for frictionless systems, and for these, the body motions have unique solutions, but the contact forces do not.

If friction is included in the model, these nice properties are no longer present. Time steppers using polyhedral friction cones still yield LCPs, but A becomes indefinite, and correspondingly, the solutions are not unique; they might not exist and when they do, there can be infinitely many. In the special case when all contacts are sticking and Coulomb friction cones are approximated by four-sided pyramids, the solutions are known to exist [29]. However, the problems are more challenging to solve if the friction model is nonlinear, as is the case for quadratic Coulomb friction, where no existence results are available. Nonetheless, some solvers have success solving the nonlinear complementarity problems (NCPs) that result. We utilize the general-purpose path solver for many NCP-based simulations. Work has been done to create specialized solvers for multibody dynamics [30].

Instantaneous-Time Model.

Consider two bodies A and B as depicted in Fig. 1. When the signed gap distance ψn is positive (as in the figure), the bodies are separated. When ψn is zero or negative, they are in contact or interpenetrating. Ideally, the latter is not permitted, that is, ψn0. Generally, a contact between two bodies may be represented as a five-tuple C={Aid,Bid,pa,n̂,ψn}, where Aid and Bid are body indices or identifiers for A and B, respectively, pa is the contact point on A in world coordinates, and n̂ is the contact normal vector in the direction of A onto B by convention. The point of contact pb on B is available by pb=pa+ψnn̂. The vectors ra and rb span from the center of mass of each body to the point of contact on that body. While not required, we use the convention that the origin of a body's body-fixed reference frame is coincident with its center of mass. Vectors t̂ and ô form an orthonormal basis with n̂.

We enforce that the motion of rigid bodies satisfy the Newton–Euler equation Display Formula

(3)M(q)ν˙=G(q)λ+λext

where M(q) is a block diagonal inertia matrix where the ith block is the (6×6) inertia matrix of the ith body

Mi(q)=[miI300Ji(q)]

with mass mi and (3×3) inertia tensor Ji(q), ν˙ is the time derivative of the vector of generalized velocities ν with νi=[viTωiT]T, λext is the vector of generalized forces acting on the bodies (gravity, velocity product, and centripetal forces), λ is the vector of contact force magnitudes (one element for each contact), and G(q) is composed of the corresponding contact Jacobians where for the jth contact

Gij=[n̂ijrij×n̂ij]

where rij is the vector from the center of mass of body i to the point of application of the force due to the jth contact as in Fig. 1 above. The generalized body configurations are represented by q=[xTqT]T for position x and unit quaternion q [26].

In addition to satisfying the Newton–Euler equation, rigid bodies should not interpenetrate. The standard approach to representing nonpenetration constraints originated from the work of Signorini on problems of contact in solid mechanics [31]. For the case of rigid-body dynamics, the Signorini condition, also known as the “normal complementarity condition,” is written as Display Formula

(4)0ψn(q,t)λn0

for a contact where ψ(q,t) is the signed gap distance which is dependent on the configuration of the bodies q and time t, and λn is the magnitude of the force applied in the normal direction n̂ of the contact. Just as we will do for G and M below, we drop the functional dependence (·) to simplify notation.

The orthogonality condition of Eq. (4) enforces for each contact that if the force being applied is positive (compressive), then the gap distance must be zero, and if the gap distance is positive, then the force must be zero. The inequalities imply that the force cannot be tensile and the bodies cannot overlap. In the frictionless case, the continuous-time dynamic model of rigid bodies in contact is composed of Eqs. (3) and (4).

Discrete-Time Model.

Since no closed-form solution of the continuous-time model exists, a discrete-time model must be derived and solved numerically. One can discretize the Newton–Euler equation in several ways. For simplicity, we use a backward Euler approximation of the time derivatives over a time step h from step to +1 to obtain Display Formula

(5)Mν+1=Mν+Gnpn+1+pext

where p denotes an impulse over the time step, i.e., p=λh.

The gap function in general does not have a closed-form solution, but distances for various known configurations can be computed efficiently with collision detection algorithms. Therefore, the gap distance at the end of the next time step is represented as its Taylor series expansion, which we truncate after the linear term here (although truncation could be done after a high-order term if desired) to get Display Formula

(6)ψ+1ψ+GT(q)ν+1h+ψth

The nonpenetration constraint of Eq. (4) is then discretized to Display Formula

(7)0ψn+GnTν+1h+ψnthpn+10

Letting ρn+1=(ψn/h)+GnTν+1+(ψn/t), we can write a dynamic model in the form of a mixed LCP (MCP) as Display Formula

(8)[0ρn+1]=[MGnGnT0][ν+1pn+1]+[Mν+pextψnh+ψnt]
Display Formula
(9)0ρn+1pn+10

Equations (8) and (9) define the time-stepping subproblem of Stewart and Trinkle [16] without friction. The first row of Eq. (8) is the discrete-time Newton–Euler equation. The second row is merely a definition used to shorten the statement of the normal complementarity condition (9). Importantly, Eq. (9) defines the standard model of unilateral contact used widely in multibody simulators. The term ψn/h appearing in the standard model acts as a constraint stabilization term, which if removed, allows interpenetrations to occur and then remain indefinitely. Regardless of whether term ψn/h is removed or retained, the inequality on the left-hand side of condition (9) is the source of SM traps that we will see PEG eliminates.

State Update.

Equations (8) and (9) model all bodies that were found to be in contact during the time step. Bodies that were not in contact are assigned updated velocities according to

ν+1=ν+M1pext

At the end of each time step, body configurations are updated. Given a quaternion representation of rotation of q=a+bi+cj+dk, all dynamic bodies are then updated according to Display Formula

(10)q+1=q+[I30012(bcdadcdabcba)]ν+1h

Since Eq. (10) does not preserve length, each quaternion should be renormalized after each time step.

One approach to multibody simulation is to allow interpenetrations (and other constraint violations) to occur and then apply whatever hacks are necessary to catch poorly defined geometric configurations that inevitably result. Our approach attempts to prevent nonphysical configurations from occurring. In order to achieve this, we must detect potential contacts before interpenetration occurs.

Contact constraints have been traditionally modeled as fundamentally unilateral, however, representing features as infinite half-spaces at points of potential contact generates nonphysical behavior at corners where multiple finite features meet and terminate. Consider a vertex v near a corner as depicted in Fig. 2. When the collision detection routine identifies both C1=C(v,e1) and C2=C(v,e2), the vertex becomes erroneously trapped in the convex subspace of the free space surrounding the corner if both contacts are enforced. We refer to this as the standard-model trap (SM trap).

Figure 3 shows a simple situation in which two SM traps collude to generate a time-stepping subproblem with no solution. A common heuristic in collision detection is to choose contacts with deepest penetration. These contacts are then used to formulate a dynamics problem to move the simulation forward in time. Since the standard contact model requires that by the end of the next time step, both of the identified contacts must be non-negative, an impossible problem is generated. Poor contact choices lead to poorly defined dynamics problems, resulting in the all too familiar “popping” or “exploding” bodies in a simulation. Note that allowing bodies such as those in Fig. 3 to first interpenetrate before identifying contacts does not avoid this erroneous contact problem, as the same heuristic results in similarly poorly defined dynamics. Fundamental errors such as these are what lead to common ad hoc “corrections” to these problems, such as limiting the maximum magnitudes of contact forces.

Deriving PEG, the New Model.

Let us begin by solving the standard-model trap in 2D for a single vertex v approaching a corner composed of two edges e1 and e2. We identify two potential contacts C1=C(va,e1) and C2=C(va,e2) with corresponding gap distances ψ1 and ψ2. We wish to enforce that Display Formula

(11)max(ψ1,ψ2)0

allowing v to penetrate either half-space corresponding to e1 or e2, but not both.

Lemma 1.1. Given a,b,max(a,b)=a+max(ba,0)0.

Proof.

  • Case aba+max(ba,0)=a+0=a.

  • Case aba+max(ba,0)=a+(ba)=b.◻

Lemma 1.2.Given a,b,b=max(a,0)0bab0.

Proof.

  • Case a0 : max(a,0)=0b=00bab0.

  • Case a > 0: max(a,0)=ab=a0bab0.

  • Case ba=0,b0 : b=max(a,0)=a.

  • Case b=0,ba0 : a0max(a,0)=0=b.◻

Using Lemma 1.1, we rewrite and constrain Eq. (11) as

max(ψ1,ψ2)=ψ1+max(ψ2ψ1,0)0

Then using Lemma 1.2 and letting c=max(ψ2ψ1,0), we can constrain max(ψ1,ψ2)0 with Display Formula

(12)0c+ψ1ψ2c0c+ψ10

To give some insight into the slack variable c and its relationship with ψ1 and ψ2, consider the figure below (Fig. 4) and the value of c on each side of the bisector where ψ1=ψ2. A vertex in the free space right of the bisector has ψ1>ψ2ψ1ψ2>0, which requires c = 0 to satisfy the complementarity condition of Eq. (12). A vertex to the left of the bisector has ψ1<ψ2ψ1ψ2<0, which requires c=ψ2ψ1.

We now wish to include potential forces λ1 and λ2 while ensuring that only one is nonzero since allowing both to be simultaneously nonzero would be nonphysical.

Lemma 1.3.Given a,b,b=|min(a,0)|0b+ab0.

Proof.

  • Case a0 : |min(a,0)|=a=bb+a=00b+ab0.

  • Case a > 0: |min(a,0)|=0=b(a+b)b=00b+ab0.

  • Case b+a=0,b0 : a0|min(a,0)|=a=b.

  • Case b=0,b+a0 : a>0|min(a,0)|=0=b.◻

Lemma 1.4. Given a,b,a=0,b0(max(a,b)0)(max(a,b)+|min(a,0)|=0).

Proof. The only way to satisfy both max(a,b)0 and |min(a,0)|0 is for both to reach equality. When |min(a,0)|=0, it is possible that a0, so in order to have max(a,b)=0, we must have a=0,b0. The argument for the reverse of a and b is similar.◻

Using Lemma 1.3, we write Display Formula

(13)0d1+ψ1d100d2+ψ2d20

and using Lemma 1.4, we have Display Formula

(14)0c+ψ1+d1λ100c+ψ1+d2λ20

Essentially, we have taken a set of contacts, in this case C1 and C2, and replaced a logical AND condition of C1C2 with a NAND condition of C1¯C2. Geometrically, this corresponds to replacing line representations of a body with line segments. At this point, Eqs. (12)(14) correspond to the locally nonconvex (LNC) model proposed by Nguyen [32], which models the case for a single vertex against multiple edges. Although this model has been shown to contribute some improvement in physical fidelity over other methods [33], stopping here would correct the standard-model trap at the expense of permitting interpenetrations [34,35]. To illustrate this error, let us continue with a less trivial case.

Consider two bodies approaching as depicted in Fig. 5. In this configuration, our collision detection routine identifies four potential contacts. These are C1=C(va,eb1), C2=C(va,eb2), C3=C(vb,ea1), and C4=C(vb,ea2). We should not simply constrain all of these contacts, for this would result in a dual standard-model trap. It is possible to make a guess as to which contacts to enforce, but this inserts the likely possibility of generating contact forces where there are none. Of course, we could also ignore all four contacts and allow interpenetration, but this too is clearly a nonphysical approach and leads to instabilities.

Instead, let us write constraints including all four potential contacts in such a way as to allow for all physically feasible interactions while preventing interpenetration. We start by observing that constraints for C1 and C2 should both be able to be violated, but not simultaneously. This is similarly true for C3 and C4. We have seen how these constraints can be represented with Display Formula

(15)0c1+ψ1ψ2c100c2+ψ4ψ3c20c1+ψ10c2+ψ400d1+ψ1d100d2+ψ2d200d3+ψ3d300d4+ψ4d400c1+ψ1+d1λ100c1+ψ2+d2λ200c2+ψ3+d3λ300c2+ψ4+d4λ40

but we make the important observation that this model erroneously allows for ψ1 and ψ4 or ψ2 and ψ3 to be simultaneously violated. An example of such a consequence is depicted in Fig. 6. The model of Eq. (15) can be completed by including constraints to prevent such interpenetrations. That is, we additionally constrain max(ψ1,ψ4)0 and max(ψ2,ψ3)0 with

Display Formula

(16)0c3+ψ1ψ4c300c4+ψ2ψ3c40c3+ψ10c4+ψ20

Equations (15) and (16) together represent our first example of a geometrically accurate contact model. We will refer to the constraint type in Eq. (15) as intercontact constraint and the type in Eq. (16) as cross-contact constraint. These, along with unilateral constraint, form the set of three contact constraints necessary for our model.

In order to resolve the issues discussed in Sec. 2, we had to change the way we formulate contact constraints. In this section, we define three fundamental contact constraints in their general form and begin to describe how they can be applied.

Unilateral Contact Constraint.

The simplest of the three contact constraints is the common unilateral contact constraint (U -constraint) that was introduced by Eq. (4). We define a U -constraint as a function of a single contact Ci and write it as Display Formula

(17)U(Ci)

which corresponds to the complementarity problem Display Formula

(18)0ψiλi0

where λi is the force corresponding to contact Ci.

Intercontact Constraint.

An intercontact constraint (I -constraint) is that which was demonstrated in Eq. (15) for preventing the standard-model trap. An I -constraint deals with contacts of a single feature of one body against multiple features of another body. Here, we will define the general formulation of I -constraint for a vertex of A against multiple features of B where a subset of these potential contacts may feasibly result in a contact force at the end of the current time step.

Consider two sets of contacts: Cι which is composed of all contacts which could feasibly result in a contact force, and Cκ which contains all contacts that could not feasibly result in a contact force but are necessary to accurately represent the contact geometry. An I -constraint seeks to

  • constrain at least one gap distance ψ{ι,κ} of {Cι,Cκ} to be non-negative

  • permit a nonzero contact force λιi>0 for at most one contact in Cι

  • only permit λιi>0 when ψιi=0 and all other ψι,κ are nonpositive

We represent an I -constraint as Display Formula

(19)I(Cι,Cκ)

which corresponds for m=|Cι| and n=|Cκ| to the set of constraints Display Formula

(20)0cι2+ψι1ψι2cι200cι3+cι2+ψι1ψι3cι300cιm+cιm1++cι2+ψι1ψιmcιm00cκ1+cιm+cιm1++cι2+ψι1ψκ1cκ100cκn+...+cκ1+cιm++cι2+ψι1ψκncκn0
Display Formula
(21)cκn++cκ1+cιm++cι2+ψι10
Display Formula
(22)0dι1+ψι1dι100dιm+ψιmdιm0
Display Formula
(23)0cκn++cκ1+cιm++cι2+ψι1+dι1λι100cκn++cκ1+cιm++cι2+ψι1+dιmλιm0

Equations (20) and (21) enforce that at least one contact in {Cι,Cκ} is non-negative. Equations (22) and (23) allow at most one contact force to be nonzero, and only permit this force to be nonzero when all other gap distances are nonpositive. Solving such a set of constraints is possible as a linear program with complementarity constraints (LPCCs) [36], or by introducing an auxiliary variable to convert inequalities (21) into a linear complementarity conditions before using an LCP solver [32].

An I -constraint requires there be at least one contact in Cι. If there were not, then there would exist no nontrivial solution since these potential contact forces are necessary for preventing interpenetration. In the case that |Cι|=1, we are able to dramatically reduce the set of constraints by eliminating Eqs. (21) and (22) and replacing Eq. (23) with the single constraint Display Formula

(24)0cκn++cκ1+cιm++cι2+ψι1λι0

This greatly reduces the number of constraints and removes the inequality of Eq. (21) to leave a pure LCP. We see now the incentive for choosing our contact set with a single primary contact Cι, possibly using heuristics like those described in Sec. 4.

Cross-Contact Constraint.

In contrast to I -constraint which deals with a single feature of A against multiple features of B, a cross-contact constraint (X -constraint) deals with sets of contacts that may not depend on the same features, but could result in interpenetration if simultaneously violated, e.g., Eq. (16).

Given the sets of contacts Ca and Cb, a X -constraint attempts to constrain either

  • at least one non-negative ψa in Ca

    OR

  • at least one non-negative ψb in Cb

We represent a X -constraint as Display Formula

(25)X(Ca,Cb)

and can write the LCP in multiple forms. Conveniently, X -constraints can be written in the form of Eqs. (20) and (21). That is, X(Ca,Cb) is similar to I(Ca,Cb), but excludes Eqs. (22) and (23) which involve contact forces.

Although a X -constraint necessarily deals with contacts between the same pair of bodies, the contacts need not share features. That is, the geometric interdependencies of identified potential contacts, regarding which are to be enforced as defined by an X -constraint, may be between two otherwise seemingly independent contacts, e.g., opposite corners of the face of a box. Subsequently, it is conceivable for two potential contacts, neither of which have a potential force associated with them, i.e., contacts in two different sets of Cκ, to be X -constrained. In such cases, the dynamics problem requires sufficient U - or I -constraints so as to allow generation of forces which may contribute to achieving configurations that will satisfy all constraints. This issue, however, while an important one when employing a heuristic approach to choosing primary contacts in I -constraints, is beyond the scope of this paper.

In this section, we briefly discuss the nature of interpolytope contact and define a set of geometric tests that are useful for identifying and making heuristic choices regarding potential contacts.

We define a convex polyhedron in terms of features (Fig. 7). We may refer to a vertex either as an object v or as the vector v representing the vertex's position. Similarly, an edge object e=(vt,vh) corresponds to a vector e=vhvt. A body is composed of a set of vertices V, edges E, and faces F defined by three vertices in counterclockwise order. The outward pointing normal vector η̂ for a given face fF is perpendicular to f. We define two vectors t1 and t2 for each edge that are planar with f1 and f2, respectively, and perpendicular to e. These vectors are defined as

Display Formula

(26)t1=η̂1×et2=e×η̂2

and are useful for geometric tests when determining contacts. We additionally define a convention of edge orientation Oe as Display Formula

(27)Oe=(η̂1+η̂2)×e

for convex bodies.

As we develop these tests, we keep in mind that a contact can occur anywhere on a polytope, but the possible normal vectors are restricted by the contacted feature. Figures 812 depict the possible positions and orientations of contact normal vectors for each feature in two and three dimensions. Since a vertex does not have a uniquely defined normal, we use the normal from the edges that connected to the vertex. It is important to realize that when using a time-stepping method, we hope to identify contacts which will reflect these feasible contact normals at the end of the time step. Yet, many methods of contact identification use approaches that assume continuous or static body positions. This is why we geometrically relax the following set of tests that will prove useful for identifying and making heuristic choices regarding potential contacts.

Consider the simple example in Fig. 13 where the two vertices, va1 and va2, of body A are near an edge eb of body B. In order to achieve stability in such a configuration, it is clear that both potential contacts C1=C(va1,eb) and C2=C(va2,eb) must be considered by the dynamics formulation. If only one were included, e.g., C1, then it is possible for the other to penetrate eb given a large enough time step or velocity. We observe that the second contact normal n2 is outside of the normal cone region of va2 and does not correspond to a physically reasonable contact normal. However, this is because n2 is determined at the beginning of the time step; yet, we wish to prevent interpenetration for the end of the time step. In order to avoid interpenetrations, we must detect and consider such contacts as C2, even when they are not physically feasible at the current time step. This necessity requires that we geometrically relax the normal cone region or, equivalently, relax the definition of applicability.

Contact Applicability.

The time required to solve a time-stepping subproblem, and therefore the time required to compute a simulation, is strongly dependent on the number of contacts nc. While it may seem that the number of contacts is determined by the simulation, taking a close look at an individual discrete time step, as we do in Sec. 1.2.3, reveals that to enhance simulation accuracy and reduce instability, it is beneficial to include “potential contacts” even though they increase the problem size. We slightly modify classical applicability conditions to derive a method to choose a small set of extra contacts.

The geometric test of applicability is well known [3740] and commonly used in applications such as distance computation and collision detection for motion planning with polytopal objects. We geometrically relax the traditional motion of applicability in order to better determine which contacts to include. We have adopted notation similar to that of Latombe [38]. The basic idea is shown in Fig. 13. If the top polygon translates toward the bottom one, contact will occur between edge eb and vertex νa1, so the feature pair eb and νa1 are said to be applicable. Following Latombe's notation [38], let ω(v) be the set of all vertices connected by an edge to va. Then, we can express the applicability condition as the Boolean function Display Formula

(28)APPL{ve}(q,va,eb)=vkω(va)η̂b·vkva||vkva||0

where η̂b is the outward unit normal of eb, and the superscript ve denotes that the test is for a vertex and an edge.

Notice that the applicability test is dependent on q; more specifically, it is dependent on the relative orientation of the polygons, and not their relative translation. Later, when considering potential contacts to include in a time-stepping subproblem, we will introduce a relaxed form of inequality (28) to allow for possible rotation in a time step. In the case of our example, if the top polygon rotated enough in the clockwise direction in the time-step under consideration, contact could occur first at vertex νa2 rather than vertex νa1.

Relaxation of 2D Applicability.

Let ω(v) be the set of all vertices connected by an edge to v. We refer to va having applicability with eb when APPL{ve}(q,va,eb)ϵθ for relaxation parameter

εθ=sin(θ)

and where Display Formula

(29)APPL{ve}(q,va,eb)=minvkω(va)η̂b·vkva||vkva||

Whereas classical applicability is a Boolean function, our function returns a value as we may wish to use the value of applicability as a heuristic when comparing contacts.

Relaxation of 3D Applicability.

The 3D vertex–face case is analogous to the 2D vertex–edge case and is defined as Display Formula

(30)APPL{vf}(q,va,fb)=minvkω(va)η̂b·vkva||vkva||

Figure 14 depicts two edges, ea and eb, in contact. Just as there are vertex–face configurations that have proximity but could not result in a force, so too are there edge–edge configurations that should not be associated with potential contact forces. First, we should note that given an edge–edge contact normal n=ea×eb, we cannot immediately say if we have the correct sign for n. However, we can still determine edge–edge applicability in the following manner. Let

da1=n̂·t̂a1da2=n̂·t̂a2db1=n̂·t̂b1db2=n̂·t̂b2

We then define edge–edge applicability as Display Formula

(31)max(min(da1,da2,db1,db2),min(da1,da2,db1,db2))

Conceptually similar to Eq. (30), Eq. (31) measures the maximal extent by which the adjacent faces of an edge are pointed “in” on another body.

Contact Feasibility.

Feasibility is similar to applicability, but is dependent on feature positions as opposed to orientations. The region of feasibility is essentially the region in which we anticipate a contact between a vertex v and a facet f could occur. We define vertex–edge and vertex–face feasibility similarly. Given an angle ϕ and a small value δ greater than machine precision, a vertex with contact C(va,fb) has feasibility if

  • vaVR(f) and ψδ

    OR

  • vaVR(f) and va within region defined by ϕ below the point ϵ below the nearest point on f

    where VR(f) refers to the Voronoi region of the facet f [41]. An example of this region is depicted in Fig. 15.

We can also write an edge–edge feasibility in 3D, which checks only to see if the nearest two points pa on ea and pb on eb are within some small distance ϵ.

From the set of fundamental constraints presented in Sec. 3, which are written in terms of complementarity problems, abstract out another step and look at pairs of bodies by cases of feature pairs.

This abstraction layering is depicted in Fig. 16. At the lowest level is the complementarity problem which is utilized first as a model of unilateral contact with force, but further as a tool for generating logical NAND conditions on subsets of potential contacts. Sets of these low-level constraints can be entirely represented by abstracting to the three fundamental constraints of unilateral, intercontact, and cross-contact. In turn, the fundamental constraints are formulaically generated by the highest abstraction layer which looks only at what pairs of features between two bodies are near each other. In Sec. 5.2, we will see how precisely the same abstraction can be developed in 3D.

PEG 2D.

We consider two contact configurations in 2D: vertex–edge and vertex–vertex. We do not consider edge–edge because this case is redundantly covered by the other two configurations [38]. To be clear, edge–edge penetration in 2D is impossible if vertex–edge and vertex–vertex constraints are properly enforced, simply because edges defined and composed of vertices, and for two edges to come near each other implies that the vertices of at least one of those edges are nearing the other edge. We can formulaically generate the set of constraints for the two 2D configurations as given below.

Vertex–Edge 2D.

Given a vertex va near a single edge eb, this case is modeled by a unilateral constraint on a single contact as Display Formula

(32)U(C(va,eb))

In this case, the vertex is within a reasonably small distance ϵ of the edge, but farther than ϵ from either of the edge's vertices.

Vertex–Vertex 2D.

The 2D vertex–vertex case was depicted in Fig. 5. Using the contacts defined for this case requires the constraints Display Formula

(33)I(Cι1,Cκ1)I(Cι2,Cκ2)X({C1},{C4})X({C2},{C3})

where Cι1 contains at least one of C1 or C2, Cκ1 may contain one of C1 or C2, and similarly Cι2 contains at least one of C3 or C4, and Cκ2 may contain one of C3 or C4. As is the case for all I -constraints, in which contacts are included in Cι or Cκ these are determined using the geometric heuristics defined in Sec. 4.

PEG 3D.

There are four contact configurations in 3D: vertex–face, edge–edge, vertex–edge, and vertex–vertex. Again, we exclude certain feature pairs just as we did in 2D, namely, edge–face and face–face, because they become redundant when the other configurations are properly enforced. Edge–face is redundant since edges are composed of vertices, and for an edge to approach a face requires that either the vertices of that edge approach that face, or the edge approaches the edges at the perimeter of that face (possibly adjacent edges), reducing to either vertex–face or edge–edge, respectively, and edge–vertex in the case that the edges of the second body are adjacent. By a similar argument, face–face becomes redundant; a face, being composed of edges in turn composed of vertices, nearing another face implies that edges and vertices are approaching the other face. In all possible configurations, this results in one of the four feature combinations considered.

Vertex–Face 3D.

The vertex–face case is analogous to the 2D vertex–edge case. Given a vertex va near a face fb as depicted in Fig. 17, we include

Display Formula

(34)U(C(va,fb))

in the set of constraints.

Edge–Edge 3D.

Edge–edge contact, such as depicted in Fig. 14, is another example of a single contact with unilateral constraint and is written for an edge ea against an edge eb by Display Formula

(35)U(C(ea,eb))

Vertex–Edge 3D.

Consider Fig. 18 where a vertex va of body A approaches the edge eb of body B. There are an arbitrary number of edges (greater or equal to three) that could connect to va; however, we can limit the edge–edge contacts considered to a subset of those with relatively large magnitudes of edge orientation |ea·Oeb|. Given C1=C(va,f1) and C2=C(va,f2), the set of constraints for vertex–edge is given by

Display Formula

(36)I(Cι,Cκ)X({Cei},{C1})X({Cej},{C2})

where C1 and C2 are distributed appropriately into Cι and Cκ. Each Cei is an edge with positive edge orientation and each Cej is an edge with negative edge orientation, and both sets of edges have applicability and feasibility with eb.

Vertex–Vertex 3D.

In the case of a vertex va near another body's vertex vb as depicted in Fig. 19, we wish to avoid trapping va against the faces of B, and similarly for vb against A. This first requires the set of I -constraints

Display Formula

(37)I(Cιa,Cκa)I(Cιb,Cκb)

where Cιa and Cκa are composed of contacts between va and faces of B, and Cιb and Cκb are composed of contacts between vb and faces of A. We additionally require X -constraints for each face pair (fa, fb) in the form Display Formula

(38)X({C(va,fb),C(vb,fa)},{C(ea1,eb1),C(ea2,eb2)})

where fa has edges ea1andea2 connected to va, and fb has edges eb1 and eb2 connected to vb. These X -constraints are enforced only for edges ear and ebs, which have edge–edge applicability and feasibility. The constraint in Eq. (38) is similar to the X -constraints in Eq. (36) for the vertex–edge case, but permits a vertex to pass on either side of a face instead of either side of an edge.

Time-Stepping Formulation With PEG.

As an independent contact model, there are many methods into which we may incorporate PEG. Here however, we will use Eq. (6) to discretize the PEG contact constraints in the same manner as the Stewart–Trinkle method of Eqs. (8) and (9). Additionally, we will focus on the case where we make a heuristic choice of a single primary contact, i.e., |Cι|=1 for every set of dependent potential contacts. We have already seen the discretized unilateral constraint in Eq. (7). An I -constraint may be written in terms of velocity by substituting Eq. (6) into Eqs. (20) and (24) to get Display Formula

(39)0EI1cI+1+ψIh+GITν+1+ψItcI+10

and Display Formula

(40)0EI2cI+1+ψIιh+GIιTν+1+ψIιtpI+10

for which we also extend the Newton–Euler equations to include Display Formula

(41)Mν+1=Mν++GnIpnI+1

For mI -constraints, cI is the solution vector composed of unknown slack variables

cI+1=[cI1+1cIj+1cIm+1]

where the cardinality of each cIj+1 is k and the jth I -constraint has k + 1 potential contacts; EI1 is diagonally composed of submatrices where

EI1j=[10010111]

which is the lower triangular matrix of size k; EI2 is diagonally composed of submatrices of row vectors of ones where

EI2j=[111]

has length k, ψI we will define as

ψI=[ψI1ψIjψIm]

where each subvector ψIj is given by

ψIj=[(ψIj1h+ψIj1t)(ψIj2h+ψIj2t)(ψIj1h+ψIj1t)(ψIjkh+ψIjkt)]

ψIι is the vector of primary potential contacts

ψIι=[ψIι1h+ψIι1tψIιmh+ψIιmt]

pI+1 is a vector of unknown impulses

pI+1=[pI1+1pIj+1pIm+1]

GI is the Jacobian difference matrix diagonally composed of submatrices GIj each composed of submatrices GIjs for s from 2 to k where Display Formula

(42)GIjs=[n̂j1rj1×n̂j1][n̂jsrjs×n̂js]

and GnI is the Jacobian matrix of only primary contacts

GnIj=[n̂j1rj1×n̂j1]

A X -constraint can be written in terms of velocities and impulses and incorporated in precisely the same way as the I -constraint. The result is that the PEG model may be written (without friction) as Display Formula

(43)Mν+1=Mν+Gnpn+1+GnIpI+1+GnXpX+1+pext
Display Formula
(44)0ψnh+GnTν+1+ψntpn+10
Display Formula
(45)0EI1cI+1+ψI+GITν+1cI+10
Display Formula
(46)0EI2cI+1+ψIpI+10
Display Formula
(47)0EX1cX+1+ψX+GXTν+1cX+10

Considering that X -constraints may be represented mathematically as I -constraints with no impulses, we may simplify to

ρn+1=GnTν+1+ψnh+ψntρnI+1=GnITν+1+EI2cI+1+ψIιγI+1=GITν+1+EI1cI+1+ψI

which may then be written in the form of a mixed LCP as Display Formula

(48)[0ρn+1ρnI+1γI+1]=[MGnGnI0GnT000GnIT0EI20GIT00EI1][ν+1pn+1pnI+1cI+1]+[Mν+pextψnh+ψntψIιψI]
Display Formula
(49)0[ρn+1ρnI+1γI+1][pn+1pnI+1cI+1]0

Equations (48) and (49) constitute one possible time-stepping subproblem that incorporates PEG for accurate geometric representations of angular features of interacting rigid polytopes. As mentioned above, friction may be added easily following the discussion in Refs. [22,23].

The experiments discussed below compare simulations performed with the time-stepping methods defined by Eqs. (8) and (9) and Eqs. (48) and (49) in 2D and 3D. The only difference in the time steppers was that one used the standard contact model and the other used PEG. Both methods retained the term ψn/h for constraint stabilization, and the time-stepping subproblems were solved by path [27]. For simulations with PEG, a single primary potential contact was chosen for each set of potential contacts, which yielded the least accurate form of PEG. Despite this, the results show significant accuracy improvements when using PEG. The examples are of convex polytopes poured into nonconvex containers represented as unions of convex polytopes.

Box of Polygons.

Figure 20 shows a snapshot from the 2D experiment, which consisted of 160 simulations (ten different initial configurations for 16 different time steps). In each experiment, 20 randomly generated convex polygons were “poured” into an open box fixed in space (the same initial conditions were used for PEG as for the standard model). One polygon was dropped every 0.25 s until the simulation ended at 5.0 s. At the end of each time step, the area of overlap among all polygons was computed and used as a metric of accuracy. When simulating the standard model, the simulation frequently crashed prior to reaching 5 s of simulation time, due to SM traps that caused the solver to fail. Typically prior to failure with the standard model, configurations were reached that subsequently allowed interpenetration or nonphysically large body velocities. Pivoting solvers, such as Lemke's algorithm [24], demonstrated similar difficulties in finding physically sensible solutions in SM trap situations. By contrast, time-stepping with PEG never encountered such problems since interpenetration and overconstraining were avoided.

Figure 21 shows the median and the first and third quartiles of area of overlap error for increasing time-step size. Even for small time steps, the time stepper using PEG produced errors several orders of magnitudes smaller than the time stepper using the standard contact model. As the error with the standard model grows with the step size, the error in PEG remains imperceptibly close to machine precision.

An interesting statistic to observe is the rate of occurrence of SM traps. In other words, how frequently was PEG necessary to achieve a geometrically accurate interaction between bodies. Figure 22 depicts the number of SM traps encountered at each time step for the 2D experiment. The top curve represents the total number of contacts, and the curve below it shows the number of SM traps. It shows that PEG was necessary in approximately 25% of time steps. One should realize that as the size of the time step increases, so would the percentage of SM traps.

Box of Polyhedra.

As shown in Fig. 23, we extend the polygon pouring experiment to 3D. In the 3D simulations, the random polygons were replaced by convex polyhedra, the box became the union of five right hexahedra, and volumes of overlap replaced the error metric used in the 2D simulations. Again, there were ten sets of different initial conditions for each of 16 time-step sizes, and these sets were each used for both PEG and the standard model. A polyhedron was dropped into the box every 0.5 s until 5.0 s of simulated time was reached.

The constraint violation error defined as volume overlap is depicted in Fig. 24, which has extremely similar behavior to the 2D results. As the step size increases, the penetrations between polyhedra grow steadily larger with the standard contact model, whereas PEG virtually eliminates this error at all time steps.