0
Research Papers

Complementarity Techniques for Minimal Coordinate Contact Dynamics

[+] Author and Article Information
Harshavardhan Mylapilli

Department of Aerospace and
Mechanical Engineering,
University of Southern California,
Los Angeles, CA 90089
e-mail: mylapill@usc.edu

Abhinandan Jain

Jet Propulsion Laboratory,
California Institute of Technology,
Pasadena, CA 91109
e-mail: Abhinandan.Jain@jpl.nasa.gov

Contributed by the Design Engineering Division of ASME for publication in the JOURNAL OF COMPUTATIONAL AND NONLINEAR DYNAMICS. Manuscript received October 3, 2015; final manuscript received April 20, 2016; published online December 2, 2016. Assoc. Editor: Andreas Mueller.The United States Government retains, and by accepting the article for publication, the publisher acknowledges that the United States Government retains, a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for United States government purposes.

J. Comput. Nonlinear Dynam 12(2), 021004 (Dec 02, 2016) (18 pages) Paper No: CND-15-1324; doi: 10.1115/1.4033520 History: Received October 03, 2015; Revised April 20, 2016

In this paper, nonsmooth contact dynamics of articulated rigid multibody systems is formulated as a complementarity problem. Minimal coordinate (MC) formulation is used to derive the dynamic equations of motion as it provides significant computational cost benefits and leads to a smaller-sized complementarity problem when compared with the frequently used redundant coordinate (RC) formulation. Additionally, an operational space (OS) formulation is employed to take advantage of the low-order structure-based recursive algorithms that do not require mass matrix inversion, leading to a further reduction in these computational costs. Based on the accuracy with which Coulomb's friction cone is modeled, the complementarity problem can be posed either as a linear complementarity problem (LCP), where the friction cone is approximated using a polygon, or as a nonlinear complementarity problem (NCP), where the friction cone is modeled exactly. Both formulations are studied in this paper. These complementarity problems are further recast as nonsmooth unconstrained optimization problems, which are solved by employing a class of Levenberg–Marquardt (LM) algorithms. The necessary theory detailing these techniques is discussed and five solvers are implemented to solve contact dynamics problems. A simple test case of a sphere moving on a plane surface is used to validate these solvers for a single contact, whereas a 12-link complex pendulum example is chosen to compare the accuracy of the solvers for the case of multiple simultaneous contacts. The simulation results validate the MC-based NCP formulations developed in this paper. Moreover, we observe that the LCP solvers deliver accuracy comparable to that of the NCP solvers when the friction cone direction vectors in the contact tangent plane are aligned with the sliding contact velocity at each time step. The theory and simulation results show that the NCP approach can be seamlessly recast into an MC OS formulation, thus allowing for accurate modeling of frictional contacts, while at the same time reducing overall computational costs associated with contact and collision dynamics problems in articulated rigid body systems.

Copyright © 2017 by ASME
Your Session has timed out. Please sign back in to continue.

References

Figures

Grahic Jump Location
Fig. 1

Overview of the five contact dynamic solvers studied in this paper. Acronyms listed in the figure are defined in the Nomenclature section. The figure details the design choices that have been made while developing each of the five solvers. Contact dynamics problems are approached using the MC OS formulation and cast as a complementarity problem. Both linear and nonlinear complementarity formulations are studied. Among the many nonlinear complementarity formulations available in the literature, we consider Todorov's approach in the present study for its simplicity and ease of implementation, whereas the prox formulation, which is also widely used in the literature, is tabled for a future course of study.

Grahic Jump Location
Fig. 2

RC formulation. Absolute coordinates are used to characterize each link in the system. Each interlink hinge is modeled explicitly as a bilateral constraint as illustrated. Figure adapted from Ref. [3].

Grahic Jump Location
Fig. 3

MC formulation. A minimal set of coordinates are used to characterize the dynamics of the system, which is obtained by eliminating all the interlink bilateral constraints. Only loop-closure constraints remain in the formulation. Figure adapted from Ref. [3].

Grahic Jump Location
Fig. 4

CE formulation. This formulation is an extension to the MC formulation, where not only the interlink hinge constraints but also the loop-closure constraints are eliminated. Bodies affected by eliminating the closure constraints are aggregated and modeled as compound bodies. Figure adapted from Ref.[3].

Grahic Jump Location
Fig. 5

Approximating the friction cone by a friction polyhedron using a finite number of direction vectors. The polygon is spanned by the direction vectors as illustrated in the figure. For each direction vector in the set, we include the opposite direction vector as well, so that the entire friction cone is spanned by the direction vectors (as opposed to a sector of the cone). Figure adapted from Ref. [3].

Grahic Jump Location
Fig. 9

Uniform sphere moving on a horizontal plane surface

Grahic Jump Location
Fig. 10

Polyhedral approximation of the friction cone for the MLCP set of solvers: (a) case 1 (validation case)—one of the friction cone direction vectors is aligned precisely with the tangential friction impulse direction and (b) case 2—friction cone direction vectors are misaligned with tangential friction impulse direction. The circle represents the circular friction limit set, the solid arrows denote the direction vectors, and the dotted lines denote the polygonal approximation of the circular friction limit set. The motion of the sphere is along the positive x-direction (and so is the tangential contact velocity vector), whereas the tangential friction impulse vector is along the negative x-direction. Four direction vectors that are at right angles to each other are chosen to span the friction cone.

Grahic Jump Location
Fig. 13

Initial setup of the 12-link pendulum

Grahic Jump Location
Fig. 11

Motion data for the sphere example (case 1), where the friction cone direction vectors are aligned with the tangential contact velocity: (a) time history of the generalized velocities and (b) time history of friction forces in contact coordinates

Grahic Jump Location
Fig. 12

Time history of generalized velocities for the sphere example (case 2), where the friction cone direction vectors are not aligned with the tangential contact velocity

Grahic Jump Location
Fig. 14

Snapshots of the motion of a 12-link swinging pendulum from t=0  to  t=20 s. Through the course of the simulation, the pendulum collides with itself and its surrounding environment.

Grahic Jump Location
Fig. 15

A time-history plot of the height of the last link of the pendulum computed using the five different solvers. The plots of the five solvers are superposed on top of one another with minor differences showing up at the end of the 20 s simulation. The height of the last link is measured from the origin, which is located at the midsection of the floor (of thickness 0.2 m).

Grahic Jump Location
Fig. 16

A time-history plot of the z-component of the linear velocity (vz) of the last link of the pendulum computed using the five different solvers. The plots of the five solvers are superposed on top of one another with minor differences showing up at the end of the 20 s simulation.

Grahic Jump Location
Fig. 17

A zoomed-in plot depicting the height of the last link of the pendulum at the end of the 20 s simulation. The solvers show a grouping behavior, with the MLCP set of solvers showing a similar solution and the MNCP set of solvers show a slightly different solution. The height of the last link is measured from the origin, which is located at the midsection of the floor (of thickness 0.2 m).

Grahic Jump Location
Fig. 18

A zoomed-in plot depicting vz of the last link of the pendulum at the end of the 20s simulation. The solvers show a grouping behavior, with the MLCP set of solvers showing a similar solution and MNCP set of solvers show a slightly different solution.

Tables

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In