Daniel Dopico1

Department of Computer Science and

Determination of Holonomic
Department of Mechanical Engineering,
Virginia Tech,
2200 Kraft Drive,
and Nonholonomic Constraint
Blacksburg, VA 24060
e-mail: ddopico@udc.es
Reactions in an Index-3
Francisco Gonz
alez
Department of Mechanical Engineering
Augmented Lagrangian
and Centre for Intelligent Machines,
McGill University, Formulation With Velocity
817 Sherbrooke Street West,
Montreal, QC H3A 0C3, Canada
e-mail: franglez@cim.mcgill.ca
and Acceleration Projections
Javier Cuadrado Index-3 augmented Lagrangian formulations with projections of velocities and accelera-
Laboratorio de Ingeniera Mecanica,
tions represent an efficient and robust method to carry out the forward-dynamics simula-
University of La Coru~na,
tion of multibody systems modeled in dependent coordinates. Existing formalisms,
Mendizabal s/n,
however, were only established for holonomic systems, for which the expression of the
Ferrol 15403, Spain
constraints at the position-level is known. In this work, an extension of the original algo-
e-mail: javicuad@cdf.udc.es
rithms for nonholonomic systems is introduced. Moreover, projections of velocities and
accelerations have two side effects: they modify the kinetic energy of the system and they
contribute to the constraint reaction forces. Although the effects of the projections on the
Jozsef K
ovecses energy have been studied by several authors, their role in the calculation of the reaction
Department of Mechanical Engineering and
forces has not been described so far. In this work, expressions to determine the constraint
Centre for Intelligent Machines,
reactions from the Lagrange multipliers of the dynamic equations and the Lagrange mul-
McGill University,
tipliers of the velocity and acceleration projections are introduced. Simulation results
817 Sherbrooke Street West,
show that the proposed strategy can be used to expand the capabilities of index-3 aug-
Montreal, QC H3A 0C3, Canada
mented Lagrangian algorithms, making them able to deal with nonholonomic constraints
e-mail: jozsef.kovecses@mcgill.ca
and provide correct reaction efforts. [DOI: 10.1115/1.4027671]

1 Introduction improve the numerical conditioning of the penalty equations.

Both families of formulations show a dissipative behavior, and
Several different options can be used for the parameterization
none of them achieves the exact satisfaction of the constraint
of the motion of multibody systems but, in general, a set of
equations and their derivatives. The mentioned formulations were
dependent coordinates is needed to write the dynamic equations
developed for both holonomic and nonholonomic systems.
of complex systems. Classical parameterizations in multibody dy-
In Ref. [11], several new formulations for holonomic systems,
namics are reference point coordinates [1], natural coordinates
based on the idea of projecting positions, velocities, and accelera-
[2], or relative (joint) coordinates [35]. In general, they are
tions onto the constraints manifold, were introduced. Particularly
dependent and related by means of a set of constraint equations.
interesting is the index-3 formulation based on velocity and accel-
The equations resulting from these redundant or dependent param-
eration projections that exactly satisfies constraints at position,
eterizations usually constitute a system of index-3 differential
velocity, and acceleration levels. The mentioned formulation was
algebraic equations (DAE) whose solution entails some numerical
subsequently employed in Refs. [12] and [13], but alternative
difficulties [6] that point toward the need to reformulate the prob-
projections were proposed to reduce the computational cost of the
lem. Different techniques to transform the equations of motion
original ones.
have been proposed, including one or several of the following
In Refs. [1416], an index-3 penalty formulation for holonomic
options: index reduction, index reduction with stabilization [7],
systems, with position penalty terms only, was combined with an
state space reformulation [8,9], penalty and augmented Lagran-
energy-momentum integrator. The resulting formulation does not
gian techniques [10], and projections onto the constraint mani-
need the velocity and acceleration penalty terms or the projections
folds [11]. The present work focuses on augmented Lagrangian
to stabilize the equations of motion. The resulting formulation
techniques and projections onto the constraints manifolds.
exactly conserves energy in conservative systems but it does not
In Ref. [10], the penalty formulation in the context of multi-
fulfill exactly either the constraints or their derivatives, the latter
body dynamics was introduced, based on the idea of completely
not being enforced at all. An improved version of the formulation
replacing the constraint equations of the original DAE system
using the augmented Lagrangian scheme was presented in
with some fictitious terms, which avoid the violation of the con-
Ref. [17]. The resulting formulation fully satisfies the constraints
straints and convert the original system into a set of ordinary dif-
but again their derivatives are not enforced.
ferential equations (ODEs) of smaller size than the original
The velocity and acceleration projections affect velocities and
problem. The solution of the penalty equations converges to the
accelerations, thus modifying energy and constraint forces. The
original DAE system solution for large penalty factors, although
effect of the velocity projections on energy conservation was stud-
at the expense of poor numerical conditioning. The augmented
ied in Refs. [1820], providing more general expressions for the
Lagrangian scheme was also proposed in Ref. [10] in order to
projections developed in Refs. [11] and [12] and proving under
which conditions the projections unconditionally dissipate energy.
1
Corresponding author.
In this paper, two topics not studied before are discussed: the
Manuscript received July 20, 2013; final manuscript received May 9, 2014; extension to nonholonomic systems of the index-3 augmented
published online July 11, 2014. Assoc. Editor: Dan Negrut. Lagrangian formulation with projections (ALI3-P) [11,12] and the

Journal of Computational and Nonlinear Dynamics OCTOBER 2014, Vol. 9 / 041006-1

C 2014 by ASME
Copyright V

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 effect of the velocity and acceleration projections on the con- transforms the equations of the augmented Lagrangian formula-
straint reactions and their evaluation. tion presented in Ref. [11] into the following expressions:
h i
2 Index-3 DAE Formulation for Holonomic Systems Mqdm UTq ki1 UTq aU Qdf (7a)
df
Let us assume that the configuration of a multibody system is i1 i i1
given by a set of nc coordinates q 2 Rnc , related by a set of m kn1 kn1 aUn1 ; i>0 (7b)
holonomic constraint equations.
The equations of motion constitute an index-3 system of nc m where
DAE, with the following expression:
i
  dm 1  dm   n1 dm   n (8a)
Mq
q UTq k Qt; q; q_ (1a)   i
  df 1  df   n1 df   n (8b)
Uq; t 0 (1b)  
i1 i1
  df 1  df   n1 df   n (8c)
where M Mq 2 Rnc nc is the mass matrix, Q Qt; q; q_
2 Rnc contains the generalized forces and may also include the and where i 0,1,2,..., a is a diagonal matrix that contains the
Coriolis and centrifugal effects (if the formulation needs them), penalty factors associated with the constraints, df and dm are scalar
U Uq; t 2 Rm is the array of constraints acting on the de- parameters of the generalized-a method, n is the time-step index,
pendent coordinates q, Uq 2 Rmn is the Jacobian matrix of the and i is the iteration index of the approximate Lagrange multi-
constraints, and k 2 Rm contains the Lagrange multipliers associ- pliers kn1 . These converge for i ! 1 to kn1 , which are the ones
ated to the constraints. resulting from the solution of the original index-3 DAE system.
The remaining quantities have been described in the nomencla-
3 ALI3-P formulation for Nonholonomic Systems ture. In Eq. (7a), the inertia, external, and constraint forces are
The ALI3-P introduced here is based on the index-3 weighted between time-steps n and n 1, following the expres-
formulations described in Refs. [11] and [12] and employs the sions (8a)(8c). For initializing the Lagrange multipliers at each
1
generalized-a integrator described in Ref. [21]. The mentioned time step (when i 0), the authors suggest kn1 kn , which
formulations deal with holonomic constraints, both scleronomic usually offers better results than the options proposed in
and rheonomic, but they do not cover the nonholonomic case. The Refs. [1113].
aim of this section is to extend the mentioned formulations to the It is very important to note that, since the scheme is an index-3
case in which nonholonomic constraints are present. DAE, Eq. (7a) only include penalty terms for the constraints at
Let us consider a multibody system modeled with q 2 Rnc position-level (2). Since the nonholonomic constraints (3) are non-
coordinates related by m holonomic constraints (some of them integrable by definition, a configuration-level expression does not
rheonomic) and m  nonholonomic constraints. Only nonholonomic exist for them, and therefore, they do not appear in Eqs. (7a) and
constraints whose expression is a linear function of the general- (7b). The projections of velocities and accelerations described in
ized velocities will be considered in this paper; nevertheless, the Sec. 3.2 become the only stage at which the nonholonomic con-
extension to more general kinds of constraints is easy straints are imposed. Moreover, special attention must be paid to
the Jacobian matrix Uq , which only includes the partial deriva-
Uq; t 0 holonomic (2) tives of the holonomic constraints but not the nonholonomic con-
 q; q; straints Jacobian A. The inclusion of the latter would result in an
U _ t Aq; tq_ bq; t 0 nonholonomic (3) incorrect scheme and convergence issues in the dynamics, among
other problems.
Equation (2) is given at the position or configuration-level, while The time-stepping equations for the method are the Newmark
Eq. (3) is given at the velocity-level, and no equivalent expression expressions [23],
at position-level exists for them.
The velocities of the system have to fulfill Eq. (3) and the first c ^_
q_ n1 q q (9)
derivative of Eq. (2), given by bh n1 n

_ 1 ^
Uq; _ t Uq q_ Ut 0
q; (4) n1
q q q (10)
bh2 n1 n

where Uq 2 Rmnc is the Jacobian matrix of the constraints, as

where
described in the Nomenclature.
Similarly, the accelerations of the system have to fulfill the sec-      
ond derivative of Eq. (2) and the first derivative of Eq. (3) given ^_  c q c  1 q_ c  1 h
q qn (11)
n
by expressions bh n b n
2b

q; q;
U _ q U
; t U q q _ q q_ U
_t0 (5)    
_ q; q; ^  1 q 1 q_ 1  1 q
q n (12)
U _ q; t Aq A_ q_ b_ 0 (6) n
bh2 n
bh n
2b

The formulation developed in this section exactly satisfies all where h is the time step and b and c are scalar parameters. These
the kinematic relations (2) and (3) and their derivatives, (4)(6). can be chosen together with the generalized-a parameters df and
This will be accomplished in two stages: first, the equations of dm, according to Ref. [22], to obtain different integrators: trape-
motion stage, in which only Eq. (2) along with the dynamic equi- zoidal rule, Newmark, and HHT (HilberHughesTaylor) are par-
librium equations are enforced; second, the velocity and accelera- ticular cases of the expressions presented here.
tion projections, in which Eqs. (3)(6) are imposed. Substituting Eqs. (9) and (10) in Eq. (7a) and scaling by a factor
of bh2 =1  dm , a nonlinear system of equations for qn1 is
3.1 ALI3-P Formulation: Equations of Motion. As obtained, which can be solved by means of the following
described in Ref. [22], the use of the generalized-a integrator NewtonRaphson iteration:

041006-2 / Vol. 9, OCTOBER 2014 Transactions of the ASME

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

   i  An option to deal with the constrained problem is to use an aug-
@f qn1 i1 i
 i
qn1  qn1  f qn1 (13) mented Lagrangian formulation to transform the constrained mini-
@qn1 mization problem (17) into an equivalent unconstrained one
i1 i i1
kn1 kn1 aUn1 ; i>0 (14) 1 1 _T _
min V  q_  q_  T Pq_  q_  1U _ Tr
aU U
q_ 2 2
where i 0,1,2,...,
1 T   Tr
1U aU U  (20)
0 h2 1
2
qn1 qn hq_ n q ; kn1 kn
2 n
where r and r  are the Lagrange multipliers for the minimization
and problem, associated with the holonomic and nonholonomic con-
straints, respectively, and a is the penalty matrix for the nonholo-
 h  i nomic constraints. The necessary condition for the minimum can

 i bh2
f qn1 M qdm UTq aU ki1  Q be obtained as
1  dm df

(15) @V  _ T_ aU _ T_ r 1U
Pq_  q_  1U _ U  T_ aU U
 T_ r
q
  i @ q_ q q q
@f qn1 i 1  df
Mn1 chCn1 Pq_  q_  1UTq aU _ UT r 1AT aU  AT r
@qn1 1  dm 
q

  ii  T T
Pq_  q_ 1Uq a Uq q_ Ut Uq r
bh2 UTqq aU ki1 UTq aUq K
n1 1AT aAq_ b AT r
0 (21)
(16)
i 1  df h  ii Expression (21) represents a nonlinear system of n equations
Mn1 chC bh2 UTq aUq K with n m m  unknowns, q; _ r, and r; therefore, the equations
1  dm n1
have to be complemented with m m  velocity-level constraint
 equations (18) and (19). For the solution of the nonlinear system,
where K  Qq, C  Qq_ , and UTqq UTq is the third order the following fixed-point iteration can be used:
q
tensor formed according with the rules provided in the Nomencla- 
ture section. For the tensorvector products involving UTqq , the P 1UTq aUq 1AT aA q_ i1
rules explained in the Nomenclature also apply.
Pq_   1UTq aUt  UTq ri1  1AT ab  AT r
i1 (22a)
From Eqs. (13) and (14), it can be noticed that, for efficiency
reasons, the same iteration, i, is used to update the generalized ri1 ri 1aU
_ (22b)
coordinates q and the Lagrange multipliers k . From a theoretical
i1 i 
point of view, the augmented Lagrangian scheme consists of solv- r
 r
 1aU (22c)
ing Eqs. (7a) and (8b) for a constant set of Lagrange multipliers
and, once the convergence is attained, the update formula (7b) is As in the augmented Lagrangian schemes described in Sec. 3.1,
applied and the dynamics equations solved again for the new mul- the update equations (22b) and (22c) are intended to make the
tipliers. This process has to be repeated until the Lagrange multi- Lagrange multipliers r and r  converge together with the fixed-
pliers converge. Nevertheless, since the proposed integrator is point iteration (22a).
implicit, advantage can be taken from the corrector iteration of If a penalty formulation is used for the projection instead of an
the integrator to update the multipliers, thus obtaining the conver- augmented Lagrangian one, a computationally less demanding
gence of both quantities simultaneously. noniterative scheme can be obtained

P 1UTq aUq 1AT aA q_ Pq_   1UTq aUt  1AT ab (23)
3.2 ALI3-P Formulation: Velocity and Acceleration Pro-
jections for Nonholonomic Systems. The aim of this section is
to properly extend the projections presented in Refs. [11] and [12] 3.2.2 Acceleration Projection. Similar to the velocity projec-
to nonholonomic systems in order to force the dependent veloc- tion, the problem to solve for the accelerations can be stated as
ities and accelerations coming from the integration stage of
1
Sec. 3.1 to fulfill conditions (3)(6). min V q q  T Pq
q  (24)
2
3.2.1 Velocity Projection. The problem to solve for the veloc-  
s:t : q; q;
1U _ q; t 1 Uq qU _ qq U _t 0 (25)
ities can be formulated as
:  
 q; q;
1 U _ q; t 1 Aq A_ q_ b_ 0 (26)
1
min V q_  q_  T Pq_  q_  (17)
2
_
  where q are the acceleration estimates coming from the solution
s:t : 1Uq; q;
_ t 1 Uq q_ Ut 0 (18) of the equations of motion after the NewtonRaphson iteration
 q; q;
1U _ t 1Aq_ b 0 (19) (13) and (14) and q are the accelerations resulting from the projec-
tion. Again, in order to simplify the notation, the subindex n 1 is
avoided in this section.
where q_  are the velocity estimates coming from the solution of
The problem in Eqs. (24)(26) can be transformed into the fol-
the equations of motion after the NewtonRaphson iteration (13)
lowing unconstrained one:
and (14), q_ are the velocities resulting from the projection, P is
the weight matrix (or projection matrix), and 1; 1 are scalar con- 1 1 T
stants for the weighting of the constraints in the minimization min V  q
q  T Pq
q  1U Tj
aU U

q 2 2
problem. Observe that, in order to simplify the notation, the subin-
dex n 1 indicating the time-step in which all the quantities are 1 _ T _ T
_ j
1U aU U  (27)
evaluated is avoided in this section. 2

Journal of Computational and Nonlinear Dynamics OCTOBER 2014, Vol. 9 / 041006-3

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 where j and j  are the Lagrange multipliers for the minimization 1  df 2
1 1 bh (35)
problem, associated with the holonomic and nonholonomic con- 1  dm
straints, respectively.
Proceeding as in the case of the velocity projection the follow- The expressions obtained in this case represent an evolution of
ing iterative solution is obtained: Cuadrados projections [12] and will be called evolved Cuadra-
 dos projections in this work.
P 1UTq aUq 1AT aA q i1 Considering that the weight of the constraints Eq. (33) is much
  smaller than the weight Eq. (31), it is apparent that Bayos projec-
P q  1UTq a U _ q q_ U _ t  UT ji1
q tions will satisfy much better the constraints than Cuadrados for
  the same number of iterations. This is something to take into
 1AT a A_ q_ b_  AT j i1 (28a)
account, especially if a noniterative penalty scheme is chosen for
ji1 ji 1aU
(28b) the projections.


j i1

ji _
1aU (28c)
3.3 ALI3-P Formulation: Implementation. The present sec-
tion summarizes the steps required to implement the formulation
Again, an inexpensive noniterative scheme is possible, less
accurate but more efficient than Eq. (28) (1) (1) n 0; t t0 ; q0 ; q_ 0 compatible with the constraints.
 0 (e.g., by means of Eq. (7a) with k 0).
(2) Calculate q
P 1UTq aUq 1AT aA q
    (3) Time-step loop: n n 1; t t h.
 T _ q q_ U_ t  1AT a A_ q_ b_
q  1Uq a U
P (29)
0 0 0
(4) Predictor: qn1 ; q_ n1 ! Eq. (9); q
n1 ! Eq. (10).
3.2.3 Selection of the Projection Matrix. There is a number of
possibilities to select the projection matrix P and constants (5) Iteration loop: i 0
1 and 1. This selection strongly affects the behavior of the projec- i1
tions. In the case of the velocity projection, the choice has a useful (6) If i > 0, kn1 ! Eq. (14).
physical meaning in terms of energy dissipation, as described in
Ref. [18]. There are two interesting options: (7) Residual f i ! Eq. (15) and tangent matrix @f i =@q
! Eq. (16).
i1 i1
(1) The mass-orthogonal projections of Bayo and Ledesma (8) Corrector: qn1 ! Eq. (13); q_ n1 ! Eq. (9); q i1
n1
[11]: they are easy to extend to the nonholonomic case by ! Eq. (10).
choosing i
(9) Until convergence:
  i 1, go
 
to 6. As convergence
 i1 i 
PM (30) criteria, qn1  qn1  < e or f i  < e.
(10) Velocity Eq. (22) [or Eq. (23)] and acceleration projec-
1 1 1 (31) tions Eq. (28) [or Eq. (29)].
(11) Until the end of the simulation: go to 3.
The resulting expressions will be called Bayos projections
in this work. 4 Calculation of the Constraint Reactions
(2) The mass-stiffness-damping-orthogonal projections of Cua- for the ALI3-P Formulation
drado et al. [12]: the aim of these projections was to obtain
The velocity and acceleration projections of Sec. 3.2 modify
the same coefficient matrix for the projections and for the
the velocities and accelerations coming from the dynamic equilib-
tangent matrix of the equations of motion, in order to take
rium stage of Sec. 3.1. As a consequence, the reaction forces
advantage of the previous factorization of this matrix to
obtained at the end of the NewtonRaphson iteration cannot be
carry out the projections in a very efficient manner. In Ref.
considered accurate since the modification of velocities and accel-
[13], the expressions were developed for holonomic systems
erations performed in the projection step must be considered for
with the trapezoidal rule as integrator. The extension to non-
the correct evaluation of the constraint reactions. This is espe-
holonomic systems with the generalized-a integrator is
cially true if nonholonomic constraints exist, since they are not
straightforward, and the following relations provide equiva-
imposed at any level by the formulation of the equations of
lent expressions for the more general case presented here:
motion (7a).
1  df   In this section, the calculation of the constraint reactions for the
PM chCn1 bh2 Kn1  1AT aA (32) ALI3-P formulation of Sec. 3 is presented.
1  dm
1  df 2
1 1 bh (33) 4.1 Constraint Reactions Coming From the Dynamic
1  dm
Equilibrium Stage. The constraint forces resulting from the
dynamic equilibrium stage of the ALI3-P formulation can be esti-
This choice makes the coefficient matrix of systems (22a) and mated by the following expression obtained from Eq. (7a):
(28a) equal to the approximate expression of the tangent matrix h i
(16), and therefore, the previous factorization of the tangent ma- Qdyn T 
U  Uq k aU (36)
trix, carried out to solve the system of equations (13), can be used. df
If the exact expression of the tangent matrix is used instead of
Eq. (16), the following selection provides the same advantage: where the subscript df means weighted evaluation between time-
steps n and n 1 like in Eq. (8b). The Lagrange multipliers kn1
are the Lagrange multipliers after convergence of Eq. (7). There-
1  df fore, the use of the iteration subscript i was avoided here for
PM chCn1 simplicity.
1  dm
h i  Note that the nonholonomic constraints do not contribute at all
bh2 UTqq aU k K  1AT aA (34) to the reaction forces of equation (36). Moreover, the expression
n1 is also valid if the constraint equations are fulfilled only up to

041006-4 / Vol. 9, OCTOBER 2014 Transactions of the ASME

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 some extent and there remains a violation of the constraints at the the reaction forces estimated from the penalty noniterative projec-
position-level, which is a common fact due to the reduced number tions Eq. (29).
of iterations of equation (7b). Observe that if the constraint equa- Making use of Eq. (41), a more convenient expression can be
tions (2) are completely satisfied, Eq. (36) would become obtained for the constraint forces coming from the projection of
h i accelerations
Qdyn
U  Uq k
T 
(37)  
df
Qacc
U Mqq  UTq 1aU
j
 
which is the classical result of expressing the reaction forces as a  j
 AT 1aU   P  Mq
q (42)
linear combination of the directions provided by the transpose of
the Jacobian matrix, UTq with the multipliers k as linear combina-
tion coefficients, playing the role of generalized magnitudes of It is important to remark that Eq. (42) is still the contribution of
these reactions. all the constraints to the generalized constraint forces vector,
Expression (36) represents the generalized constraint forces nevertheless from the two first terms of the right hand side, one
resulting from all the constraints. Nevertheless, the contribution of can identify separate contributions of individual constraints.
all the constraints together is not usually required, but the contri- Unfortunately, due to the last term, it is not possible to exactly
bution of each one of them separately. The contribution of con- extract separately the contribution of each constraint for any
straint r to the generalized constraint forces is the following: choice of the projection matrix P, as will be shown in Secs. 4.2.1
and 4.2.2 in which Eq. (42) is particularized for Bayos and
h i
T  evolved Cuadrados projections of Sec. 3.2.3.
Qdyn
Ur  Ur q kr ar Ur (38)
df
4.2.1 Constraint Reactions for Bayos Projections. Using
Eqs. (30) and (31) in Eq. (42), the constraint reactions can be
To illustrate the interest of using expression (38) instead of Eq. directly obtained as
(36) consider the 1 DOF system consisting of three masses joined
by two constant distance constraints and constrained to move only    :
Qacc U T
aU j  AT a U
 j
 (43)
in one direction. The equations of motion (1a) of such a system U q
are as follows:
2 32 3 2 3 2 3 If the projection iteration has fully converged and the satisfaction
m1 0 0 x1 1 0   F1 of the constraint equations (5) and (6) is good, then the terms
4 0 m2 0 54 x2 5 4 1 1 5 k1 4 F2 5 (39) coming from the constraints can be neglected and a simplified
k2
0 0 m3 x3 0 1 F3 expression is obtained as

If one is interested in the reaction forces Qacc T T

U Uq j  A j (44)
 acting on the second
T
mass, the result offered by Eq. (36) is Uq k 2 k1 k2
Comparing expressions (44) and (37), it is clear that the physi-
which is the sum of the forces exerted by both rods. Equation
cal interpretation of the Lagrange multipliers of the projections j
(38), instead, allows for distinguishing
 between the forces of the and that of the Lagrange multipliers of the dynamics k is exactly
first and the second rod U1 Tq k1 2 k1 and the same for Bayos projections.
 From the expressions of the reactions of Bayos projections
U2 Tq k2 2 k2 . (43), the contribution of each constraint can be easily identified

T 
Qacc
Ur Ur q ar Ur jr (45)
4.2 Constraint Reactions Coming From the Projection
Stage. The projections cause changes in the velocities and accel- 
erations of the system in order to make them compatible with the Qacc T
 r Ar a
_ r j
r U r (46)
U
constraints. Considering the linearity of the equations of motion
Eq. (7a) with respect to the accelerations, any change in the accel- where Eq. (45) is for individual contributions of holonomic
erations to fulfill the constraint Eqs. (5) and (6) can be interpreted constraints and Eq. (46) is for individual contributions of nonholo-
in terms of missing constraint reaction force terms. Thus, the nomic constraints.
missing reactions coming from the acceleration projection (28a)
can be expressed as
4.2.2 Constraint Reactions for the Evolved Cuadrados
Qacc
U Mq 
q (40) Projections. Using Eqs. (32) and (33) in Eq. (42), the constraint
reactions can be directly obtained as
Equation (40) gives us the total contribution of all the constraint
reactions to all the generalized coordinates. In this sense, it is  
T 1  df 2
equivalent to Eq. (36), and therefore, it is not useful to distinguish Qacc
U U q bh a U j
the contribution of an individual constraint to the reactions, as 1  dm
 
explained before. The following derivation intends to obtain 1  df 2 :
expressions equivalent to Eq. (38), which make it easier to iden-  AT 
bh a U j
1  dm
tify the contribution of each constraint separately. From Eq. (28a)  
1  df  
 :  chCn1 bh2 Kn1  1AT aA q  (47)
q
  1  dm
Pq  UTq 1aU
q j  AT 1a U
 j
 (41)

Note that in Eq. (41), it is important to keep both terms, the multi- Again, in the case in which convergence of the projection itera-
pliers and the acceleration constraints, so that it is not necessary tion has been achieved and the satisfaction of the constraint equa-
to achieve the full convergence of the projection process. Even tions (5) and (6) is good, the following simplified expression is
the use of Lagrange multipliers can be completely avoided and obtained:

Journal of Computational and Nonlinear Dynamics OCTOBER 2014, Vol. 9 / 041006-5

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Qacc T
U Uq j  A j
T
 reactions in the case in which only velocity projections are used.
  The exact relations are those presented before based on the accel-
1  df  
 chCn1 bh2 Kn1  1AT aA q 
q eration projection.
1  dm One small limitation of the calculation proposed, which is
(48) really a limitation of the formulation itself, is revealed by the fact
that the velocity projection does not play any role in the reaction
Unlike Bayos projections, for evolved Cuadrados projections, forces. For the case of velocity-dependent generalized forces
the physical interpretation of the Lagrange multipliers of the pro- Qq_ , it is clear that the velocity projection should have an effect
jections j and the Lagrange multipliers of the dynamics k is not on the dynamics of the system and thus on the reactions. Unfortu-
the same. j are simply multipliers that enforce the fulfillment of nately, the formulation developed never recalculates the general-
the constraints imposed, but no physical interpretation in terms of ized forces vector after projecting the velocities so the mentioned
reactions can be extracted from them. effect, even if small, is not handled neither by the formulation,
In contrast to the results obtained for Bayos projections, from nor by the formulae here proposed to calculate the reactions.
the expression of the reactions obtained using the evolved Cuadra-
dos projections (47), the contribution of each constraint cannot
5 Numerical Results
be extracted separately without assuming some additional approx-
imations: considering that a is a diagonal matrix with very large The method to evaluate the constraint reaction forces described
penalty factors in the diagonal and h a small time-step, the terms in Sec. 4 was assessed in the forward-dynamics simulation of
chCn1 and bh2 Kn1 can be neglected compared with AT aA, and three test examples.
therefore, the contribution of each constraint can be approximated
by the following formulae: 5.1 Two Point-Masses System With Nonholonomic Con-
  straint. The first example consists of two point masses moving on
T 1  df 2 the xy plane (Fig. 1). The motion of this system can be described
Qacc
Ur   U r q bh ar U r j r (49)
1  dm with four generalized coordinates xA, yA, xB, and yB. In this case,
  the velocity of the mass at point A is constrained to be pointing to-
T 1  df 2  _ j
Qacc

Ur
 A r bh 
ar U r r ward B at every instant by the velocity-level constraint
1  dm (50) yB  yA x_ A  xB  xA y_A 0. It can be shown that this con-
 T 
1Ar ar Ar qq  straint is nonholonomic, since no equivalent position-level expres-
sion can be found for it. Both masses are initially moving with
where Eq. (49) represents individual contributions of holonomic constant velocity x_A x_ B 1 m=s, while y_A y_B 0. A force
constraints and Eq. (50) individual contributions of nonholonomic fB 2 N is acting in the vertical direction on point mass B. A 5 s
constraints. long forward-dynamics simulation was carried out, using two dif-
ferent methods to solve the dynamic equations. The first one was
the direct solution of the system resulting from expressing (1a)
4.3 Total Constraint Reactions. Finally, the total general-
and (1b) with a Lagrangian formulation, establishing Qc AT k
ized constraint forces, considering the contribution from both the
equations of motion and the projections and including all the     
M AT q Q
constraints, are (54)
AT 0 k A_ q_  2nxU

QU Qdyn acc
U QU (51)
The trapezoidal rule was used as integrator, with an integration
time-step h 103 s.
where Qdyn
U is given by Eq. (36); Qacc
U in general is given by
Second, the simulation was repeated with the same parameters
Eq. (42), being the particular expression for Bayos and evolved
using the index-3 augmented Lagrangian formulation and Bayos
Cuadrados projections (43) and (47), respectively.
projections with P M and 1 1. The reaction forces associated
As explained before, Eq. (51) represents the contribution of all
with the nonholonomic constraint were evaluated using the results
the constraints to the generalized constraint forces. It is generally
of the projection process as described in Sec. 4. Results are shown
more useful to identify contributions of individual constraints.
in Fig. 2. The direct solution method and the augmented Lagran-
Thus, the total generalized reaction force associated with a single
gian formulation yielded the same constraint reaction forces.
constraint r, is
Moreover, the obtained results with the augmented Lagrangian
formulation were the same, regardless of whether the reaction
QUr Qdyn acc
U r QU r (52) forces were obtained from the projection of accelerations or from
QU r Qacc
r (53) the projection of velocities.
U

where Eq. (52) is the expression for holonomic constraints and 5.2 Three-Dimensional Rolling Disk. The second test prob-
Eq. (53) the one for nonholonomic constraints; Qdyn lem is a 3D rigid disk rolling on the horizontal xy plane (Fig. 3)
Ur is given by
Eq. (38); Qacc
Ur and QU
acc
 r are given by Eqs. (45) and (46), respec-
tively, for Bayos projections and they are given by Eqs. (49) and
(50), respectively, for evolved Cuadrados projections. It is impor-
tant to note that, from Eq. (53), only the projections contribute to
the generalized constraint forces of nonholonomic constraints.
One important fact to note about the generalized reactions
described before is that they are obtained from the dynamic equi-
librium stage and the acceleration projection stage exclusively.
The velocity projection stage does not play any role in the expres-
sions obtained. In the Appendix, new approximate formulae to
determine the constraint reactions from the velocity projections
are presented as additional material for the reader, but they are Fig. 1 Nonholonomic constraint restraining the velocity of
just approximate relations that allow to have an estimation of the mass A to be collinear with segment AB

041006-6 / Vol. 9, OCTOBER 2014 Transactions of the ASME

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 Fig. 2 Reaction forces acting on mass A in the x and y direc-
tions due to the nonholonomic constraint in the first example,
computed with direct solution (54) and with the ALI3-P
formulation Fig. 3 3D motion of a disk rolling on a plane

under gravity effects. The disk has mass m 5 kg and radius

q 0.5 m. This system was modeled using a set of 12 natural
coordinates [2], related by a set of constraint equations, defined
both at position and velocity-levels. The set of natural coordinates
comprised the x, y and z coordinates of the center of the disk (P)
and three orthogonal unit vectors u1 (normal to the plane of the
disk), u2 and u3 (both contained in the plane of the disk). These
vectors and point P form the local reference frame of the disk. At
time t 0 s, vectors u1 ; u2 , and u3 are aligned with the inertial ref-
erence frame x, y, and z. Six kinematic constraints were intro-
duced to ensure that the set of coordinates behaves as a rigid
body. The norm of vectors u1 ; u2 , and u3 was constrained to be 1,
and the angles between them were made constant during motion.
Moreover, the z coordinate of the contact point location C was
constrained to be 0 with Fig.4 Reaction fy at contact point C during 2D motion of the
disk, obtained with the projection of velocities (fy vel) and
q accelerations (fy acc)
zC zP  q 1  u21z 0 (55)
Second, the initial velocity of point P was set to
where zC and zP are the z coordinates of points C and P and u1z is x_P 0:05 m=s; y_P 1:3 m=s, and z_P 0:0 m=s. The reaction
the z component of vector u1 . The rolling condition can be forces fx and fy obtained using Eq. (51) are shown in Fig. 5. They
expressed as are compared to m xP and myP for verification. Results show that
the obtained reaction forces verify fx mxP and fy m yP , which
vC vP  qx  d 0 (56) are the values required by the motion developed by the system.
where x is the angular velocity of the disk and
u1  n  u1 5.3 Marine Crane M95.20. The last test problem is the ma-
d (57) rine crane M95.20 from Industrias Guerra shown in Fig. 6.
k u1  n k
If the rotation of the crane is neglected, the system has three
where n is the unit vector normal to plane xy. The first two equa- DOF with hydraulic actuators between pairs of points 34, 58,
tions in Eq. (56) are nonholonomic, and they were used to impose and 1011. The mass properties of the system were provided by
x_C 0 and y_C 0. The third one is the velocity-level equivalent the manufacturer. The hydraulic actuators were modeled by three
to Eq. (55). The imposition of rolling introduces a reaction force rheonomic constraints with a prescribed motion given by time-
at the contact point with components fx, fy, and fz, which cannot be dependent functions:
correctly evaluated with the previously existing index-3
algorithms. s1 t 0:990 0:25t2 m (58)
Two 5 s forward-dynamics simulations of the motion of the 2
s2 t 1:013 0:25t m (59)
disk were carried out with an integration time step h 5  104 s.
2
The augmented Lagrangian formulation of index-3 and the s3 t 0:076 0:20t m (60)
evolved Cuadrados projections with dm 0, df 0, c 1/2, and
b 1/4, described in Secs. 3 and 4, were used to integrate the The dynamics of the system under the prescribed motion was
dynamic equations of the system. First, a constant force f 3 N solved for 1 s of simulation time, using the ALI3-P formulation
was applied in the global y direction to point P, between t 1 with Bayos and evolved Cuadrados projections. As it was com-
and t 3 s. Under these conditions, the motion is two-dimensional mented before, the fulfillment of the constraints in Bayos projec-
in the plane yz, and the reaction force at point C can be deter- tions is much better than in evolved Cuadrados for the same
mined using the free body diagram of the wheel to be fy 1 N. number of iterations. Thus, in this problem for Bayos projections,
Figure 4 shows that the forces obtained from the projections a noniterative penalty scheme is sufficient to obtain accurate reac-
match the expected result. tions while for evolved Cuadrados an iterative augmented

Journal of Computational and Nonlinear Dynamics OCTOBER 2014, Vol. 9 / 041006-7

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

 this case is an iterative augmented Lagrangian scheme to guaran-
tee accurate results.
Results show that the classical calculation from the Lagrange
multipliers of the dynamics without taking into account the pro-
jections can offer completely inaccurate results even for holo-
nomic systems. This limitation is overcome if the formulation
described in this paper is used.

6 Conclusions
Existing ALI3-P formalisms were not designed to tackle nonho-
lonomic constraints. The need for velocity and acceleration pro-
jections, which are essential for the stability and accuracy of the
formulation, introduces constraint reaction forces required to ful-
fill the constraint derivatives, which are not imposed by the equa-
tions of motion alone.
Fig. 5 Reaction forces fx and fy at contact point C during In this work, the extension of the ALI3-P formulation to nonho-
motion of the disk, compared to mxP and myP lonomic systems was accomplished, and a simple and efficient
technique to obtain the reaction forces coming from both the
equations of motion and the projections was presented. Results
obtained from the forward-dynamics simulation of three examples
showed that the method is able to provide correct values of the
generalized constraint forces associated with holonomic and non-
holonomic constraints.

Acknowledgment
The first author acknowledges the support by Award No. NSF
CMMI-1130667 and by the Computational Science Laboratory at
Virginia Tech.
The second and fourth authors acknowledge the support by the
Natural Sciences and Engineering Research Council of Canada
(NSERC) and CMLabs Simulations, Inc.
Fig. 6 Three DOF crane

Nomenclature
 
Dq @D @D @D

@q1 @qj @qnc
2 Rsrnc
third order tensor of derivatives of matrix
D 2 Rsr with respect to vector q
Dq v Dq v
 
@D @D @D
v v v 2 Rsnc
@q1 @qj @qnc
tensor-vector product, where v 2 Rr is a
vector
Mq 2 Rnc nc mass matrix
q 2 Rnc vector of coordinates of the system
_ t 2 Rnc
Qq; q; vector of generalized forces, including Corio-
lis and centrifugal effects if present
t time
(.) d d2 @ @
; 2 ; t ; q
dt dt @t @q

Fig. 7 Actuation forces: (a) Bayos (top) and (b) evolved Cua-
drados (bottom)
Appendix: Approximate Value of the Constraint Reac-
tions Obtained From the Velocity Projection Formulae
Lagrangian scheme is mandatory to obtain acceptable results for Since velocities and accelerations are related by the integrator
the reactions. equations and both velocity and acceleration projections impose
The actuation forces corresponding to the first rheonomic con- exactly the same constraints, it is possible to approximately obtain
straint Eq. (58) using Bayos projections are shown in Fig. 7(a), the value of the constraint forces from the velocity projection.
where the contribution of the different terms and the total force From Eq. (21)
were represented. The theoretical value of the actuation force is    
also represented in Fig. 7(a), calculated by means of the virtual Mq_  q_  UT 1aU
q
_ r  AT 1aU  r
work principle, in order to compare with the expressions devel-  P  Mq_  q_  (A1)
oped in this work. The scheme chosen for these projections is the
cheapest, noniterative penalty scheme. Velocities and accelerations coming from the equations of
The same forces solving with evolved Cuadrados projections motion stage and from the projections stage are, respectively,
are shown in Fig. 7(b). The scheme chosen for the projections in related by means of the following Newmark formulae:

041006-8 / Vol. 9, OCTOBER 2014 Transactions of the ASME

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use

   
q_  q_ n h1  c q 
qn c (A2) 1 1  df 2 
Qvel
r  
U
ATr bh ar Ur r
r
q_  q_ n h1  c
qn c
q (A3) ch 1  dm
  
 1ATr ar Ar q_  q_  (A11)
where q_ n and q
n are velocities and accelerations evaluated in the
previous time step. Observe that Eq. (A2) is exact, because it
relates magnitudes after projecting in the previous time step with
magnitudes before projecting in the current time-step, which is where Eq. (A10) is for individual contributions of holonomic con-
directly Eq. (9). On the other hand, Eq. (A3) is approximate straints and Eq. (A11) is for individual contributions of nonholo-
because it relates magnitudes after projecting in the previous time nomic constraints.
step with magnitudes obtained after projecting in the current time
step and they are not directly related by the integrator equations.
Substituting Eqs. (A2) and (A3) into Eq. (A1). References
[1] Haug, E. J., 1989, Computer Aided Kinematics and Dynamics of Mechanical
Systems, Allyn and Bacon, Inc., Boston, MA.
Mfq_ n h1  c
qn c
q  q_ n  h1  c q g
qn c on, J., and Bayo, E., 1994, Kinematic and Dynamic Simulation of
[2] Garca de Jal

   
ch Mq q T
 U 1aU _ r  A 1aU
T  r  Multibody Systems. The Real-Time Challenge, Springer, New York.
q [3] Paul, B., and Krajcinovic, D., 1970, Computer Analysis of Machines With Pla-
 P  Mq_  q_  (A4) nar Motion, ASME J. Appl. Mech., 37(3), pp. 697712.
[4] Sheth, P., and Uicker, J., 1972, IMP (Integrated Mechanism Program): A
Computer-Aided Design Analysis System for Mechanisms and Linkages,
Then, the constraint forces obtained from the velocity projec- ASME J. Eng. Ind., 94(2), pp. 454464.
tion become [5] Smith, D., Chace, M., and Rubens, A., 1973, The Automatic Generation of a
Mathematical Model for Machinery Systems, ASME J. Eng. Ind., 95(2), pp.
629635.
1 n T _    [6] Brenan, K., Campbell, S., and Petzold, L., 1989, Numerical Solution of Initial-
Qvel
U Mq   
q  r
Uq 1aU r AT 1aU  Value Problems in Differential-Algebraic Equations, North-Holland, New York.
ch
o [7] Baumgarte, J., 1972, Stabilization of Constraints and Integrals of Motion in
P  Mq_  q_  (A5) Dynamical Systems, Comput. Methods Appl. Mech. Eng., 1(1), pp. 116.
[8] Serna, M. A., Aviles, R., and Garca de Jalon, J., 1982, Dynamic Analysis of
Plane Mechanisms With Lower Pairs in Basic Coordinates, Mech. Mach.
Theor., 17(6), pp. 397403.
Equation (A5) is an approximation to Eq. (42). Particularizing [9] Wehage, R., and Haug, E., 1982, Generalized Coordinate Partitioning for
the previous expression for Bayos projections Dimension Reduction in Analysis of Constrained Mechanical Systems, ASME
1 n T _  
 r
o J. Mech. Des., 104, pp. 247255.
Qvel
U   Uq aU r AT aU  (A6) [10] Bayo, E., Garca de Jalon, J., and Serna, M. A., 1988, A Modified Lagrangian
ch Formulation for the Dynamic Analysis of Constrained Mechanical Systems,
Comput. Methods Appl. Mech. Eng., 71(2), pp. 183195.
[11] Bayo, E., and Ledesma, R., 1996, Augmented Lagrangian and Mass-
From the expressions of the reactions of Bayos projections Orthogonal Projection Methods for Constrained Multibody Dynamics, Nonlin-
ear Dyn., 9(12), pp. 113130.
(A6), the contribution of each constraint can be easily identified [12] Cuadrado, J., Cardenal, J., Morer, P., and Bayo, E., 2000, Intelligent Simula-
tion of Multibody Dynamics: Space-State and Descriptor Methods in Sequential
1n  o and Parallel Computing Environments, Multibody Syst. Dyn., 4(1), pp. 5573.
Qvel
Ur   Ur Tq ar U_ r rr (A7) [13] Cuadrado, J., Gutierrez, R., Naya, M., and Morer, P., 2001, A Comparison in
ch
Terms of Accuracy and Efficiency Between a MBS Dynamic Formulation With
1 n T  ^
o Stress Analysis and a Non-Linear FEA Code, Int. J. Numer. Methods Eng.,
Qvel
r
U
 Ar ar Ur rr (A8) 51(9), pp. 10331052.
ch [14] Garca Orden, J., and Goicolea, J., 2000, Conserving Properties in Constrained
Dynamics of Flexible Multibody Systems, Multibody Syst. Dyn., 4(23), pp.
where Eq. (A7) is for individual contributions of holonomic con- 225244.
[15] Goicolea, J., and Garca Orden, J., 2000, Dynamic Analysis of Rigid and De-
straints and Eq. (A8) is for individual contributions of nonholo- formable Multibody Systems With Penalty Methods and Energy-Momentum
nomic constraints. Schemes, Comput. Methods Appl. Mech. Eng., 188(4), pp. 789804.
The particularization of Eq. (A5) for evolved Cuadrados pro- [16] Goicolea, J., and Garca Orden, J., 2002, Quadratic and Higher-Order Con-
jections gives the following expressions: straints in Energy-Conserving Formulations of Flexible Multibody Systems,
Multibody Syst. Dyn., 7(1), pp. 329.
   [17] Garca Orden, J., and Ortega, R., 2006, A Conservative Augmented Lagran-
vel 1 T 1  df 2 _ gian Algorithm for the Dynamics of Constrained Mechanical Systems, Mech.
QU   Uq bh aU r
ch 1  dm Based Des. Struct. Mach., 34(4), pp. 449468.
  [18] Garca Orden, J. C., and Dopico, D., 2007, Multibody Dynamics. Computa-
1  df 2  tional Methods and Applications, On the Stabilizing Properties of Energy-
AT bh aU r
 Momentum Integrators and Coordinate Projections for Constrained Mechani-
1  dm
  cal Systems, Springer, Dordrecht, The Netherlands, pp. 4967.
1  df   [19] Garca Orden, J., 2010, Energy Considerations for the Stabilization of Con-
chCn1 bh2 Kn1  1AT aA q_  q_  strained Mechanical Systems With Velocity Projection, Nonlinear Dyn.,
1  dm 60(12), pp. 4962.
(A9) [20] Garca Orden, J., and Conde, S., 2012, Controllable Velocity Projection for
Constraint Stabilization in Multibody Dynamics, Nonlinear Dyn., 68(12), pp.
245257.
Assuming the same approximations mentioned in Sec. 4.2.2, [21] Cuadrado, J., Dopico, D., Naya, M. A., and Gonzalez, M., 2004, Penalty,
the contribution of each constraint can be approximated by the Semi-Recursive and Hybrid Methods for MBS Real-Time Dynamics in the
following formulae: Context of Structural Integrators, Multibody Syst. Dyn., 12(2), pp. 117132.
[22] Geradin, M., and Cardona, A., 2001, Flexible Multibody DynamicsA Finite
  Element Approach, 2nd ed., Wiley, Chichester, England.
1 T 1  df _
Qvel
Ur   U
r q bh2
a U
r r r r (A10) [23] Newmark, N. M., 1959, A Method of Computation for Structural Dynamics,
ch 1  dm ASCE J. Eng. Mech. Div., 85(EM3), pp. 6794.

Journal of Computational and Nonlinear Dynamics OCTOBER 2014, Vol. 9 / 041006-9

Downloaded From: http://computationalnonlinear.asmedigitalcollection.asme.org/ on 01/27/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use