Simulating mechanical systems on computers requires a model for describing the dynamics of the mechanical parts. Models that can describe the deformation of the mechanical parts require a high number of parameters to describe the deformation. If the core of the mechanical parts can be assumed to not deform under the considered loads, the parameters describing the state of a mechanical part can be reduced to that of a rigid body: An invariant shape with an associated mass and principal moments of inertia as well as the spatial orientation, position, linear and angular velocity of the shape. The interaction of multiple such mechanical parts must be described by another model determining the dynamics of the mechanical parts in contact, where the contact model usually allows compliance in a localized contact region. If mechanical parts collide, the contact dynamics typically occur on a time-scale that is significantly smaller than that of the motions between successive collisions. Resolving each such collision micro-dynamics in computer simulations can become computationally expensive. Alternatively, the relation between the pre- and post-collision state variables can be described by an impact modelstronge04. By condensing the impact dynamics to a single point in time, only the response of the relative contact velocities has to be specified. In order to instantaneously turn a colliding state into a non-colliding state a contact reaction impulse must be applied, so that the necessary discontinuities in the velocities can be effected. The impact model must then be combined with a non-compliant contact model in order to determine the contact reaction forces and contact reaction impulses. Non-compliant contact models alone cannot resolve collisions. For non-compliant contact models with Coulomb friction even non-colliding contact situations exist (shocks) where an impulse becomes necessary to resolve the contact stewart00. These paradoxical situations were first published by Painlevé painleve1895. The system including the non-compliant contact and impact models are often mathematically described in terms of measure differential inclusions (MDI) moreau88.
The rigid body simplification combined with the impact simplification considerably reduce the computational burden in simulations with many mechanical parts. Applications range from robotics nuseirat00; jia13, virtual reality sauer98, physics-based animation erleben04 to granular matter simulations tasora10. In particular, granular dynamics can require a very large number of particles and are insufficiently understood to date with and without an interstitial fluid phase mitarai12. Optimizing mechanical devices like powder mixers hassanpour11 or grinding mills mishra92; jayasundara11 is of economical importance. Powder mixing is important in detergent, cosmetic, food and pharmaceutical manufacturing, to name just a few applications. Understanding granular matter is also of crucial importance for safety reasons: Assessing the stability of slopes is important to prevent rock slides, land slides and snow avalanches and getting the particle distribution right in pebble-bed nuclear reactors is important to guarantee safe and performant operation tasora10.
The construction of an impact model that in combination with a non-compliant contact model leads to a consistent theory for rigid body dynamics is non-trivial. Ideally, the solution of an impact model for a collision is a limit point of a sequence of solutions of the collision based on a compliant contact model with increasing stiffness. The sequence of solutions is uniquely determined and the increasing stiffness decreases the duration of the collision towards an instantaneous event. The solution of the collision based on a compliant contact model corresponds to the integral of the contact reaction forces over the collision duration. Using such an approach Stronge constructs an energetically consistent restitution hypothesis in stronge90, Mirtich solves rigid-body dynamics with impacts for virtual reality applications, where permanent contacts are treated as sequences of collisions mirtich95; mirtich96. And lately, Jia and Wang showed in jia16 how contact reaction impulses can be computed from the limit of a contact model with linear normal stiffness and Coulomb friction for general collisions (central or eccentric, direct or oblique) in three dimensions. The authors established a condition that, if met, guarantees solution existence. Whether solutions exist unconditionally and whether the impact model in combination with a non-compliant contact model leads to a consistent theory remain open problems.
However, other impact models exist that possess solutions unconditionally. Stewart showed in stewart98 for a time-stepping scheme based on anitescu97, that it converges to a solution of an MDI as the time-step size decreases. Recently, Gavrea et al. extended the result to systems including joints in gavrea08. The MDI describes rigid-body dynamics with Coulomb friction and purely inelastic impacts (collisions and shocks), where the frictional impulses are required to directly oppose the post-impulse relative contact velocities in the tangential planes thus imitating Coulomb’s friction law in the case of impulses. A consequence of this is that Stewart proved that solutions exist for the MDI and thus resolved paradoxical configurations in rigid-body dynamics with non-compliant contacts and Coulomb friction, where apparently no solutions exist even though no collisions are present. Stewart made no attempt to show uniqueness of solutions. In fact in section LABEL:sec:non-unique we present an example demonstrating the existence of multiple solutions of a numerically constructed single-contact impact problem. The non-uniqueness is in this single-contact case directly related to the choice of the frictional impact model. In this paper we construct an alternative frictional impact model having a unique solution in the single-contact case. The friction model is based on the maximum dissipation principle stewart00 and takes into account the coupling between the normal component and the tangential components of the contact reactions. The non-uniqueness due to redundant constraints in the multi-contact case remains unaffected popa15 as well as non-uniqueness in the non-compliant contact model.
In section 2 we present integral equations describing rigid-body dynamics with impact and friction, where the Coulomb friction model on the impulsive reactions is replaced by a friction model based on the maximum dissipation principle. In section 3 the model is embedded into an impulse-velocity time-stepping scheme for numerically integrating multi-contact problems and an analytic solution of the single-contact problem is established. Subsequently, section LABEL:sec:results presents results for single-contact problems and the macro-scale behaviour of the friction model in the simulation of a large-scale granular flow problem. The paper summarizes the results and concludes in section LABEL:sec:summary.
2 Continuous System
Each particle is associated with a co-rotating body frame. The origin of the body frame corresponds to the center of mass of the particle. Let be the position function of the body frame in the inertial frame. The position function is non-smooth when impulses act. The derivative with respect to time is the discontinuous linear velocity function with left- and right-limits and
. The orientation of the body frame in the inertial frame can be represented by a unit quaternion. Instead of mixing vector and quaternion algebra a quaterniondescribing the orientation of the body frame of particle at time is represented as a vector . The orientation function is non-smooth and the left- and right limit of the derivative at time is then diebel06
where is the angular velocity of the particle beforeafter applying impulses at time . We also introduce the quaternion matrix function for abbreviating the notation.
The mass of the particle is denoted by
and is invariant with respect to time. The inertia tensor of the particle in the inertial frame changes with respect to the orientation of the body frame. It can be expressed in terms of the constant inertia tensor in the body frame. For time the inertia tensor in the inertial frame is given by
where is the rotation matrix corresponding to the orientation . The rotation matrix changes the basis from the particle’s body frame to the inertial frame. Choosing the body frame such that the axes match the principal axes of the particle, the body frame inertia tensor can be enforced to be diagonal.
In a system with particles, let (, , and ) be the vertical concatenation of all particles’ positions (orientations, linear velocities, and angular velocities) at time :
and let be the block-diagonal mass matrix containing in the upper-left quadrant and in the lower-right quadrant, where is the
identity matrix. Then given initial conditions at time, the state of the system at time is described by the integral equations
where is the set containing all points in time , where impulses are present, that is (linear) impulse or angular impulse is non-zero:
The terms and are the forces and torques acting on the particles. Note that the inverse of the mass matrix always exists, since it is symmetric positive-definite (SPD) - a property which it inherits from its diagonal blocks. The cross-product term is to be understood as the vertical concatenation of all single-particle cross-products:
The appearance of the term stems from the fact that the torque function corresponds to the time-derivative of the angular momentum function (for non-impulsive points in time), which in turn is the product of the time-varying inertia tensor and the angular velocity. Hence,
The forces, torques, linear impulses, and angular impulses at time include components from non-impulsive contact reactions and impulsive contact reactions , where is the number of contacts in the particle system. Each contact involves a pair of particles , . By convention let contact reactions act positively on the first particle and negatively on the second particle . Each contact is also associated with a contact frame. Let the first axis of the contact frame correspond to the contact normal pointing from particle towards particle
by convention, and let orthonormal vectors, and complete the contact frame. Let denote the position of the contact frame in the inertial frame. Then, subsuming all forces and torques on particle , which are not due to contact reactions, as external forces and external torques , the equations
define a wrench matrix function relating the wrenches to the contact reactions. This relation extends to impulsive reactions and linear and angular impulses:
The impulsive and non-impulsive contact reactions are then given implicitly as solutions of contact constraints. The contact constraints are usually non-linear and underdetermined depending on the specific contact model employed.
The formulation of the contact constraints requires the rigorous introduction of the contact position function , the contact normal function , the signed contact distance function and the relative contact velocity function . The latter is straightforward and for a contact given by
It can be shown that
The definition of the other three functions are difficult to state in sufficient generality. We confine ourselves here to definitions that are at least well-defined for spherical particles. Let be the set of points in the inertial frame defining the shape of particle at time , and let be the signed distance function associated with the shape . The signed distance function shall be negative in the interior of the shape. Then, let
be the contact point between the pair of particles . If the boundary of the shape is sufficiently smooth and the overlap sufficiently small, the contact position is uniquely determined and the gradient of the signed distance function exists. Then the contact normal is given by
The signed contact distance function is then simply
These specific definitions of the contact functions limit the number of contacts to the number of particle pairs . To simplify the description of the contact constraints, subscript denotes the projection of a vector to the contact normal (e.g. ) and subscript denotes the vector of projections of a vector to the contact tangential and contact orthogonal (e.g. ).
Then the contact constraints for an inelastic contact with Coulomb friction are listed in Fig. 1.
|Non-penetration constraints||Coulomb friction constraints|
The constraints can be classified into non-penetration constraints and Coulomb friction constraints. Both classes can be subdivided into impulsive and non-impulsive constraints. Impulsive constraints determine impulsive contact reactionsand non-impulsive constraints determine non-impulsive contact reactions . Some constraints are understood to be enabled only if a precondition holds. For instance the restitution hypothesis should only constrain the solution if the contact is closed (). This precondition is indicated by an arrow. The arrow originates from a constraint, which enables the constraint if the precondition becomes active. In the case of the restitution hypothesis the contact needs to close first and the arrow thus originates from the impulsive Signorini condition. The Signorini condition ensures that contact reactions are non-negative (non-adhesive) if the contact is closed and are zero if the contact is open. This relation is expressed by the complementarity condition and the corresponding inequalities. A similar chain of non-penetration constraints exists for the non-impulsive contact reactions. However, in that chain the contact reaction in the worst case can only be determined after the constraint on the acceleration level became enabled.
The impulsive and non-impulsive Coulomb friction constraints require the contact reaction to reside within a friction cone. The coefficient of friction determines the aperture of the cone. The cone is aligned along the contact normal. The friction cone condition limits the Euclidean norm of the frictional reaction by an upper bound proportional to the contact reaction in normal direction. The direction of the frictional reaction is required to oppose the relative tangential contact velocity in the case of a sliding (dynamic) contact and its Euclidean norm must be at its limit. This is expressed in the velocity-level equation. In the case of a sticking (static) contact, the velocity-level equation is universally valid. However, the zero slip enables the acceleration-level constraint. Then, the direction of the frictional reaction is required to oppose the relative tangential contact acceleration.
The work performed by the frictional contact reaction force of contact for a non-impulsive time span is
where the Coulomb friction force performs no work to the extent possible as expressed in the acceleration-level constraint. However, if sliding is inevitable, the Coulomb friction force maximizes dissipation by minimizing the integrand through the velocity-level constraint, which requires the friction force to directly oppose the relative contact velocity. The velocity-level Coulomb constraint can be formulated equivalently using the maximum dissipation principle as pointed out by Stewart in his review paper on friction and impact in rigid-body dynamics stewart00 at least if the normal reaction is considered to be given stewart11:
where the objective function corresponds to the (negated) rate of energy dissipation. Since the relative contact velocity at time is independent of the contact reaction at time , the objective function is a linear function of the frictional contact reaction at a non-impulsive point in time .
When formulating the friction constraint on the impulsive contact reactions, the situation changes subtly but drastically: The impulsive contact reactions now influence the post-impulse relative contact velocity. The drastic consequence of this is that the maximum dissipation principle and the Coulomb friction model for that matter as it is formulated for non-impulsive contact reactions in Eq. (1) cannot be transferred to impulsive contact reactions without in-depth modifications if the property of maximizing the energy dissipation is to be preserved. This statement stands in contrast to common practice jean99; stewart00; bonnefon11. In particular the term is a quadratic function of the impulsive contact reactions and it does not reflect the energy dissipated. Hence, impulsive frictional reactions directly opposing the relative contact velocity in the tangential plane also do not dissipate as much energy as allowed by the friction cone condition in general.
The system energy is the sum of the potential energy and the kinetic energy :
where the potential energy is not affected by impulses. Insertion leads to the expression for the post-impulse system energy
in terms of the pre-impulse system energy and the impulsive contact reactions. Let be the Delassus operator and let condense the terms depending linearly on the impulsive contact reactions and let condense the constant terms. A contact reaction complying with the maximum dissipation principle should minimize . Restating in terms of the -th contact reaction and assuming all other contact reactions to be constant results in
where corresponds to the -th diagonal block of the Delassus operator and where selects all columns (elements) except for column (element ). The diagonal block can be determined to be mirtich96; preclik14
where and .
Then the impulsive contact reaction complying with the maximum dissipation principle is
where Eq. (2) corresponds to the friction cone condition, Eq. (2) corresponds to the Signorini condition, Eq. (2) corresponds to the purely inelastic restitution hypothesis, and Eq. (2) is an additional constraint requiring that the contact reaction is not increasing the contact pressure. The last constraint guarantees uniqueness for a single contact. It excludes non-zero solutions if the contact opens by itself. The objective function is a quadratic function of the contact reactions and it is strictly convex since is SPD.
For open contacts () the Signorini constraint and the friction cone constraint restrict the feasible set to the reaction . The restitution hypothesis is disabled and the pressure constraint is fulfilled.
If the contact is closed () the Signorini condition reduces to and the restitution hypothesis is enabled. In order to determine whether the restitution hypothesis is active () or inactive (), the dependence of on the impulsive contact reactions must become explicit. At least for spherical particles the time-derivative of the post-impact signed distance function can be expressed in terms of the relative contact velocities:
Contacts fulfilling the property , that is contacts where a penetration is imminent if no impulsive contact reaction acts, are termed colliding in the following. Pre-impulse velocities, external impulses and impulsive reactions from other simultaneously colliding contacts determine if the contact is in a colliding state.
Contacts fulfilling the property , that is contacts where separation is imminent if no impulsive contact reaction acts, are termed separating. For separating closed contacts, the reaction fulfills all constraints and thus the restitution hypothesis is inactive. The pressure constraint ensures that no non-zero solutions exist.
For colliding closed contacts, the restitution hypothesis must be active, restricting the feasible set to the plane of maximum compression defined by Eq. (4). Combined with the friction cone condition, the feasible set forms a conic section. The normal of the plane of maximum compression is . Since is SPD and since the contact is colliding, the conic section is guaranteed to be non-empty. Since the conic sections are non-empty convex sets and since the objective function is strictly convex, the optimization problem has a unique global minimum. The pressure condition is fulfilled since it is fulfilled for any point on the plane of maximum compression in the colliding case.
The unconstrained global minimum is given by
If the contact is colliding and closed and if fulfills all constraints, then and the post-impulse relative contact velocity in the tangential plane is zero corresponding to a sticking (static) post-impulse contact state. If does not fulfill all constraints the post-impulse relative contact velocity in the tangential plane is non-zero and the post-impulse contact state thus sliding (dynamic). The solution of the contact problem then resides on the boundary of the conic section.
The friction model for contact impulses from Eq. (3) thus maximizes dissipation by minimizing the (kinetic) energy analogously to the velocity-level constraint of the Coulomb model for non-impulsive contact reactions. This is in contrast to the alleged extension of the Coulomb model to impulsive reactions, where the dissipation is not (sufficiently) maximized. At the same time the contact impulses fulfill the friction cone condition as in the Coulomb model and they fulfill the inelastic restitution hypothesis. The friction model is lazy in the sense that separating solutions (no contact reaction) are preferred over any other solutions, which is guaranteed through the pressure condition. Furthermore, static solutions (no work performed) are preferred over dynamic solutions analogously to the acceleration-level constraint of the Coulomb model.
3 Numerical Methods
3.1 Multi-Contact Problems
Integrating the equations describing the rigid-body dynamics can be approached in at least two different ways: The event-driven approach aims to predict the next impulsive point in time given an initial state at time . Then the simulation is integrated until assuming the impulsive contact reactions to be zero for . At time an impact problem given by
and Eq. (3) has to be solved. Then, the integration can be restarted having a state fulfilling all constraints at hand. The difficulty of this approach lies in the problem of predicting the next impulsive point in time, which is a priori unknown. If rigid bodies follow ballistic trajectories impact times can be predicted accurately and efficient simulation codes exist mirtich96; bannerman11. For the more general case, where e.g. impulsive points in time can not only stem from collisions but also stem from self-locking sliding frictional contacts shen11, we know of no efficient method to accurately predict the next impulsive point in time. Furthermore, simulation codes necessarily stall in situations, where the collision frequency increases unboundedly like in cases where a bouncing ball comes to rest on a plane.
The second category of approaches for integrating rigid-body dynamics are the time-stepping methods. These methods proceed in time steps independent of the impulsive points in time. The methods do not distinguish between non-impulsive and impulsive contact reactions but implicitly solve for integrals of the contact reactions . The integrals of the contact reactions are then constrained to fulfill contact conditions at selected points in time. In the following an impulse-velocity time-stepping scheme is used, which is described in detail in preclik14. It is similar to the schemes by Anitescu anitescu97, Tasora tasora11 and Stewart stewart96. The equations of motion are integrated using a discretization similar to a semi-implicit Euler method, where positions are integrated implicitly and velocities are integrated explicitly:
where primes identify state variables at time in contrast to state variables at time omitting the prime. The relative contact velocities are
The discrete non-penetration constraint for a contact is
which if the gap is closed exactly () states that the relative contact velocity at the end of the time step and in the direction of the contact normal must be non-approaching and complementary to the non-adhesive contact reaction in normal direction. The term acts as an error reduction term if penetrations are present (). In that case it can be scaled down when needed using an error reduction parameter to avoid introducing an excessive amount of energy (). If a positive gap is present (), the term ensures that contact reaction remains zero if the contact would not close within the time step. The friction cone condition
is adopted as it stands. Instead of requiring the impulsive contact reactions to fulfill the conventional Coulomb friction constraints as in
the maximum dissipation principle from Eq. (3) is used, since the discretized relative contact velocity in the tangential plane depends on the contact reactions as before. Therefore, the energy term is discretized, leading to
where . For a single contact the energy term reduces to
given all other contact reactions . The relative contact velocity in terms of and is
Hence, the contact constraint complying with the maximum dissipation principle is
Then a solution of the multi-contact time-step problem satisfies Eq. (5) for all contacts . The problem can be solved using a non-linear block Gauss-Seidel or variants thereof as demonstrated in preclik14; preclik15. The solution algorithm proposed there reduces the multi-contact problem to the problem of solving a sequence of single-contact problems. Hence, in the next section an analytic solution for the single-contact impact problem from Eq. (5) (and Eq. (3) alike) is derived.
3.2 Single-Contact Problems
The single-contact problem for impulsive contact reactions complying with the maximum dissipation principle from Eq. (3) constrains solutions to zero if contacts are open (). If contacts are closed, Eq. (3) has the same structure as the discrete single-contact time-step problem from Eq. (5):
where , is SPD and
If the contact is separating () the solution is constrained to since the pressure condition () prevents any non-zero solutions. If the contact is non-separating () the pressure condition is redundant and the constraint set is formed by the intersection of the plane of maximum compression and the friction cone. This conic section is non-empty but can take on the shape of an ellipse, parabola or hyperbola. The cases where are special in the sense that the plane of maximum compression includes the origin and thus the conic section degenerates to a point, a ray or a degenerate hyperbola. These cases will be discussed after the cases where . The cases where and are also special since the conic section degenerates to a single (non-zero) point . If and the conic section is non-degenerate. The unconstrained minimum of the objective function is . If it fulfills the friction cone condition it is contained in the constraint set and it thus solves the contact problem .
If the unconstrained minimum of the objective function does not fulfill the friction cone condition, the solution must be located on the boundary of the conic section minimizing the objective function. Eliminating the normal component using the equation for the plane of maximum compression and switching to polar coordinates leads to
where and . The friction cone condition then becomes
For angles where , the inequality poses no additional restrictions on the non-negative coordinate . For angles where , the inequality defines an upper bound on the coordinates . In both cases the component satisfies
Let be the set of angles for which holds. can be determined to be
where and . Let
describe the curve along the boundary of the conic section, where . Then the contact problem reduces to
result from eliminating the normal component. The objective function is univariate and -periodic but no longer strictly convex nor quadratic.
Fig. 4 illustrates the optimization problem for an exemplary (dynamic) contact. Fig. 4 plots the quadratic objective function for points on the plane of maximum compression and overlays the conic section. Fig. 4 plots the corresponding non-linear objective function in comparison. An iterative approach for solving this constrained minimization problem that is guaranteed to converge linearly is derived in preclik14. In the following an analytic approach is presented.
be the unit vector tangential to the curve, where
Angles minimizing must satisfy
Insertion, trigonometric identity transformations and multiplication by leads to a trigonometric equation in the form of
with constants as specified in appendix LABEL:sec:appendix. The trigonometric equation can be transformed into a quartic equation by substituting . After solving the quartic equation for , the corresponding angles have to be checked for validity. Angles are invalid as well as angles not satisfying Eq. (7) or Eq. (8) for that matter. Among all other candidates, the one with minimum objective function value amounts to the -periodic solution . Finally,
Back to the case, where and : If the conic section corresponds to a degenerate ellipse (), then the solution is exactly this point. The conic sections that correspond to a degenerate parabola () are equivalent to the ray
Degenerate hyperbolas correspond to the conical combination of two rays (along the asymptotes)