Research Papers

Simulation of Constrained Mechanical Systems — Part I: An Equation of Motion

[+] Author and Article Information
David J. Braun1

 School of Informatics, University of Edinburgh, 10 Crichton Street, Edinburgh EH8 9AB, UKdavid.braun@ed.ac.uk

Michael Goldfarb

 Department of Mechanical Engineering, Vanderbilt University, VU Station B 351592, Nashville, TN 37235,michael.goldfarb@vanderbilt.edu

In an ideal computational environment, the constraints would be exactly satisfied, which would result in, Δq=0, Δν=0, and Aqv=bq. In that case, q·=v, and Eq. 17 would reduce to q··=a+R-1Cν+(bν-Aνa)+R-1NνR-TQ0. The last two terms in this equation originate from the general solution of the linear Eq. 10. The same general solution obtained from a linear Eq. 9 would result in q··=a+M-1/2B+(bν-Aνa)+M-1/2(I-B+B)M-1/2Q0, where M=M1/2M1/2 and B=AνM-1/2. This equation is the one proposed by Udwadia and Kalaba [11].

In an ideal computational environment, the kinematic constraints and the energy law would be exactly satisfied, (i.e., Δeν=0), and thus the constrained acceleration could be exactly computed q··=a+R-1Cν+(bν-Aνa). In that case, Eq. 20 would reduce to Q0=RT[CeNν]+Ae(q··-a-R-1Cν+(bν-Aνa)), and as such it would not perform any correction, Q0=0.

The numerically exact solution is an ideal reference along which the constraints are exactly satisfied.

The constraint accuracy is reported by depicting the error on the kinematic constraints.

The energy accuracy is reported by depicting the energy error.


Corresponding author.

J. Appl. Mech 79(4), 041017 (May 16, 2012) (9 pages) doi:10.1115/1.4005572 History: Received November 11, 2010; Revised October 21, 2011; Posted January 31, 2012; Published May 16, 2012

This paper presents an equation of motion for numerical simulation of constrained mechanical systems with holonomic and nonholonomic constraints. In order to avoid the error accumulation typically experienced in such simulations, the standard equation of motion is enhanced with embedded force and impulse terms which perform continuous constraint and energy correction along the numerical solution. To avoid interference between the kinematic constraint correction and the energy correction terms, both are derived by taking the geometry of the constrained dynamics rigorously into account. In this light, enforcement of the (ideal) holonomic and nonholonomic kinematic constraints are performed using ideal forces and impulses, while the energy conservation law is considered as a nonideal nonlinear nonholonomic constraint on the simulated motion, and as such it is enforced with nonideal forces. As derived, the equation can be directly discretized and integrated with an explicit ODE solver avoiding the need for expensive implicit integration and iterative constraint stabilization. Application of the proposed equation is demonstrated on a representative example. A more elaborate discussion of practical implementation is presented in Part II of this work.

Copyright © 2012 by American Society of Mechanical Engineers
Your Session has timed out. Please sign back in to continue.



Grahic Jump Location
Figure 1

Slider-crank mechanism:  l = 1  m, m = 1  kg. The simulation is performed for t ∈ [0, 100]s started from q(0) = [2/4, 2/4, π/4, 32/4, 2/4,-π/4]T, v(0) = [0, 0, 0, 0,0, 0]T and conducted with 5 × 10 - 2 s time step. The stroboscopic view represents the motion for t ∈ [96.9, 99]s.

Grahic Jump Location
Figure 2

Horizontal motion of the slider, xs = x2 + (l/2)cos(θ2) vs time. The picture shows that the solution obtained with the proposed DAE integration method (with 5×10 - 2  s time step) is well matched with the numerically exact ODE solution.

Grahic Jump Location
Figure 3

From top to bottom the semilog pictures depict: |Φq|, |Φ·q|, |Φe|, and ɛ = dx2+dy2, where dx and dy are measured with respect to the numerically exact ODE solution. The integration here is performed with 10 - 2  s time step.

Grahic Jump Location
Figure 4

From top to bottom the log-log pictures depict the maximum constraint error, the maximum energy error and the maximum overall error. The accuracy/step-size characterizations are performed through a t ∈ [0, 100]s simulation. On the last plot, the CPU time corresponding to integration with and without energy correction is presented.

Grahic Jump Location
Figure 5

Double four-bar linkage: l = 1m, m = 1 kg for each link. Whenever the links become horizontal, y1 = y2 = y3 = 0, the linkage moves through a kinematic singularity.

Grahic Jump Location
Figure 6

Double four-bar linkage: l = 1m, m = 1 kg for each link. The simulations are performed through t ∈ [0, 1000]s, with time step 10 - 2 s and q(0) = [0, 1, 1, 1, 2, 1]Tv(0) = [1, 0, 1, 0, 1, 0]T. The ODE solutions are obtained using the unconstrained formulation of the problem: θ·· + (7g/6l)cos(θ) = 0, θ(0) = π/2, θ·(0) =  - 1 (θ denotes the orientation of the base link, see Fig. 5).

Grahic Jump Location
Figure 7

Evolution of the (least accurate) kinematic constraints, the energy conservation law and the overall accuracy ɛ3 = dx32 + dy32 where dx3, dy3 is measured between the solutions depicted in Fig. 6, and the numerically exact solutions for (x3, y3)




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