See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/235770791

Developing a Numerical Simulation Software for 3D Multibody Systems Based on

a Unified Computational Modeling Technique

Article · August 2009

DOI: 10.1115/DETC2009-86550

CITATION READS

1 165

2 authors, including:

Shahram Shokouhfar
Concordia University Montreal
5 PUBLICATIONS 13 CITATIONS

SEE PROFILE

All content following this page was uploaded by Shahram Shokouhfar on 20 May 2014.

The user has requested enhancement of the downloaded file.

 Proceedings of the ASME 2009 International Design Engineering Technical Conferences &
Computers and Information in Engineering Conference
IDETC/CIE 2009
August 30 - September 2, 2009, San Diego, California, USA

DETC2009-86550

DEVELOPING A NUMERICAL SIMULATION SOFTWARE FOR 3D MULTIBODY

SYSTEMS BASED ON A UNIFIED COMPUTATIONAL MODELING TECHNIQUE

1 2
Shahram Shokouhfar Sayyid Mahdi Khorsandijou
Mechanical Engineering Group Mechatronics Group
Islamic Azad University-Parand Branch, Iran Islamic Azad University-South Tehran Branch, Iran

ABSTRACT
This article represents the features and capabilities of a INTRODUCTION
newly developed application namely MASS (Mechanisms It is quite possible to derive equations governing the
Analysis and Simulation Software) and the formulation and kinematics and dynamics of a particular 3D multibody system
techniques therein. MASS is a general C++ application program with the aid of the minimum possible number of the generalized
whose main task is to construct and solve the governing coordinates. The derived equations are useful only for a
algebraic differential motion equations of 3D multibody particular type of systems and are numerically solved with
systems automatically in matrix forms complying with the minimum amount of calculations. To have the general
computational algorithms required for numerical simulation. governing equations of 3D multibody systems sacrifices the
Newton-Raphson and SVD methods have been used for advantage of the minimum possible number of the generalized
kinematical assembling and producing consistent initial coordinates. Moreover, some kinematical constraints among the
conditions. Adaptive time-step Runge-Kutta-Fehlberg numerical generalized coordinates emerge and increase the amount of
integration methods might be used for forward dynamics calculations required for numerical simulation. The kinematical
problems. The governing equations perfectly describe the constraints represent the restrictions caused by the kinematical
kinematics and dynamics of multibody systems within which 3D joints of the system.
kinematical joints and collisions between rigid bodies might be In order to develop a general computational tool for
taken into consideration. The unified computational technique kinematics and dynamics analysis of 3D multibody systems, e.g.
for mathematical modeling of kinematical joints is the most 3D mechanisms, coordinate viewpoint has been taken into
important concept on top of which MASS has been consideration [1]. The coordinate viewpoint has led to the
implemented. It has occurred due to the existence of thirteen creation of thirteen basic kinematical constraints [1] for unified
basic kinematical constraint equations. Each kinematical joint dynamic modeling of the systems. It should be noted that the
might be defined by a set of algebraic equations being selected unification of the dynamic models has occurred due to the
from the mentioned basic equations. The unified dynamic existence of the thirteen basic equations. Point, line and vector
models for collisions and impulsive loads have been also viewpoints have been considered in [2-6] and four basic
achieved using the mentioned technique. Simulation results kinematical constraints have been used for dynamic modeling of
obtained from MASS have been compared with that of the general mechanical systems. Numerical analyses required for
corresponding software of Working Model ver. 6 and a the kinematics and dynamics simulation of the systems having
discussion about the coincidences and differences has been multi rigid or flexible bodies have been perfectly studied in [4].
exposed. The present article exposes the technique of thirteen basic
kinematical constraints required for unified dynamic modeling
Keywords: Unified modeling technique; Basic constraints; of 3D multibody systems in a C++ application, namely MASS.
Kinematical joints; Collision; Multibody systems; Dynamics; Based on the thirteen constraints, some outstanding novelties
Numerical simulation; Software development; MASS. have arisen and proposed in this article. Using SVD (Singular
Value Decomposition) method [7], MASS can receive

1
Instructor of Mechanical Engineering, Corresponding Author, Email: [email protected]
2
Assistant Professor of Mechanical Engineering

1 Copyright © 2009 by ASME

consistent initial conditions with independent and adequate n h : Number of collisions detected at t C
components, and can avoid joint violation and cumulative
integration errors. MASS might solve the forward dynamics n k1 , n k 2 : Normal unit vectors of contact points of the k1th and
problems with the aid of the adaptive time-step Runge-Kutta- k th
2 bodies being collided
Fehlberg 5-6 or 7-8 numerical integration methods as well as P : System vector of Euler parameters
other simpler methods. MASS shows the real DoF of the th T
systems and the kinetic, potential and total energies therein as a Pn : Euler parameters of the n body; ëéa n bn cn d n ûù
verifying instrument. This paper has compared simulation T
Q : System vector of generalized coordinates; éR T PT ù
results of MASS with that of Working Model ver. 6, and has ë û
declared the reasons of the existing coincidences and T
R : System vector of positions; éR1T R T2 L R TnB ù
differences. ë û
th
R i : Mass center of the i body expressed in GCS
NOMENCLATURE Rot : Rotational
A i : The ith body Rotation matrix projecting from CCS on GCS ri : Position of J i expressed in CCS of the ith body
A Ji : The ith body Rotation matrix projecting from JCS on CCS
T Free : Free angle of rotational spring
A Ji J j : Rotation matrix projecting from J j JCS on J i JCS Tr : Translational
t C : Exact instant of collision
C i : Mass center of the ith body expressed in GCS
t k1 , t k 2 : Unit vectors along relative velocities of contact
C L : Damping coefficient of reciprocating damper
C T : Damping coefficient of rotational damper points of the k1th and k th
2 bodies being collided
D mn : Direction cosines; Dependent rotational joint variables V : System vector of velocities
DoF: Degree(s) of freedom Vn : Velocity of the nth body
f Actuator : Force applied by a reciprocating actuator V p : System vector of quaternion derivatives
f A : Active forces externally applied to bodies Vnp : Quaternion derivatives of the nth body
f C : Constraint forces applied to bodies X , Y , Z : Cartesian coordinates; Translational joint variables
é -b n an dn - cn ù
th ê ú e : Infinitesimal time for impulse or collision
G n : G of the n body; ê -c n -d n an bn ú z i : Point of action for linear impulse applied to the ith body
êë -dn cn -bn a n úû
hk , hk1 , hk 2 : Contact point in the kth detected collision
G : G 3n of System; diag (G1, G2 , L, GnB )
B ´4 n B Q , F , Y : Euler angles; Independent rotational joint variables
I B : System matrix of rotary inertia with respect to each CCS
l C : Lagrange multipliers due to joints constraints
ILi : Linear impulse applied to the i body th
l p : Lagrange multipliers due to quaternions constraints
IiT : Angular impulse applied to the i body th
mk : The kth collision friction coefficient at t C
J i : Origin of JCS attached to the ith body sk : The kth collision restitution coefficient at t C
K L : Reciprocating spring constant jC = 0 : Joints constraints
K T : Rotational spring constant j K = 0 : Basic kinematical constraints
L : Rod or jack joint variable
LFree : Free length of reciprocating spring jp = 0 : Quaternions constraints
M : System mass matrix wn : Angular velocity of the nth body
m A : Active moments externally applied to the system
m Actuator : Moment applied by rotational actuator UNIFIED MODELING TECHNIQUE
m C : Constraint moments applied to bodies A unified technique for modeling of 3D joints is planned on
the basis of coordinate viewpoint. It represents the restrictions
N : System normal reaction impulses due to collisions caused by each type of the joints. The restrictions lie in a set of
N k : Normal reaction force due to the kth collision algebraic equations being selected from a comprehensive set of
n B : Number of bodies in a mechanical multibody system thirteen basic constraint equations. The unified dynamic models
n C : Number of constraint equations for collisions and impulsive loads are also achieved using the
mentioned technique.
n L : Number of linear impulsive loads
n T : Number of angular impulsive loads

2 Copyright © 2009 by ASME

Coordinate Systems éX ù é D11 D12 D13 ù
Generally, four types of coordinate systems are required. uuur ê ú ê ú
J i J j = ê Y ú , A Ji J j = ê D 21 D 22 D 23 ú (6)
Global Coordinate System (GCS) is an inertial coordinate
reference frame anchored to the ground. Body Coordinate ëê Z ûú êë D31 D32 D33 ûú
System (BCS), Central Coordinate System (CCS) and Joint Exp.'s (6) show the relative translational and dependent
Coordinate System (JCS) are body-fixed frames. BCS describes rotational variables of a joint. The relative independent
the geometric details of a body. CCS is parallel to BCS having an rotational joint variables are shown by the local xyz or global
origin coincided with mass center. Kinematical joints, springs, zyx Euler angles of Exp. (7). The length of the relative
dampers and actuators can only be defined on the basis of two translational joint variables, i.e. L , has been also considered as
JCS's separately possessed by two bodies. a rod or jack joint variable.
Position of the nth body can be described by the origin of its é1 0 0 ù é cos F 0 sin F ù écos Y -sin Y 0 ù
CCS expressed in GCS, as shown by Eq. (1). A Ji J j = êê0 cos Q -sin Q úú êê 0 1 0 úú êê sin Y cos Y 0 úú (7)
R n = [ xn zn ]
T
yn (1) êë0 sin Q cos Q úû êë-sin F 0 cos F úû êë 0 0 1 úû

The orientation of the nth body can be described by four Two jointed rigid box-shaped bodies, their JCS's and the
Euler parameters, namely quaternions, as shown by Eq. (2). joint assigned between them are shown in Figure 2 within
MASS environment.
Pn = [ a n dn ]
T
bn cn (2)

Relation between Euler parameters

Euler parameters representing the orientation of the nth
body are dependent due to the quaternions constraint Eq. (3).
The first and second time-derivatives of Eq. (3) are given by
Eq.'s (4).
jpn = PnT Pn - 1 = 0 (3)
¶jpn p ¶jpn & p
Vn = npn , Vn = g pn (4) Figure 2: A time-dependent joint in MASS environment
¶Pn ¶Pn
The Jacobian and the right hand sides of Eq.'s (4) are given To produce various types of kinematical joints, the joint
by Eq.'s (5). variables might be assigned as free, fixed, or a function of time
¶jpn T in MASS, as shown by Figure 1. If a joint variable is left free,
= 2PnT , npn = 0 , g pn = -2 éë Vnp ùû Vnp (5)
¶Pn its corresponding equations are ignored. The joint variables are
implicitly expressed in terms of positions and orientations of the
jointed bodies in the thirteen basic kinematical constraints,
Basic Kinematical Constraints
namely Eq. (8) or Eq.'s (9), (11), and (12). The agent variables
A kinematical joint is defined between two bodies in order
are defined by Exp.'s (10) and (13).
to restrict their relative motions. Restriction equations should be
ùû = [ 0]13´1
T
expressed through a unified modeling technique that simplifies jK = éë j1K jK2 L j13 K
(8)
the programming stage. When the joints are defined by user, the The first three equations of Eq. (8), i.e. Eq. (9), coincides
system restriction equations are automatically formed with the the origins of the two JCS's with each other. They are
aid of selection from a set of thirteen comprehensive basic independent and might be used for creating a spherical joint for
constraints equations and substitution of the constants therein. instance.
éj1K ù éX ù
ê Kú ê ú
êj2 ú = A Ji éëA i D - ri ùû - ê Y ú = [ 0 ]3´1
T T
(9)
êj3K ú ëê Z ûú
ë û
D = R j - R i + A j rj (10)
The 4 up to 12th equations of Eq. (8), i.e. Eq. (11), assigns
th

Figure 1: Joint coordinate systems and independent variables the relative orientations of the two JCS's and eventually that of
the two jointed bodies. They are dependent and might be used
A general joint and eventually the basic constraints are for creating a prismatic joint on space for instance.
readily established by means of two JCS's J i and J j as shown é jK4 j5K jK6 ù é D11 D12 D13 ù
ê K Kú
j9 ú = A Ji A i A j A J j - êêD 21 D 23 úú = [0]3´3
th th
by Figure 1. They are respectively attached to the i and j ê j7 jK
8
T T
D 22 (11)
bodies. ê j10
K
jK K ú
j12 êë D31 D 32 D 33 úû
ë 11 û

3 Copyright © 2009 by ASME

 Only three equations out of the nine equations of Eq. (11) and Table 2. They form the joints constraints of the system.
are independent, but to cope with some numerical difficulties, According to Eq. (3), each body has one quaternions constraint,
all nine dependent equations have been considered in MASS. It so Eq. (14) is the system constraints composed of body
should be noted that the redundant equations have been quaternions constraints along with the joints constraints.
automatically skipped in a pivoting process. é jC ( R , P ) ù
The 13th equation of Eq. (8), i.e. Eq. (12), assigns the j=ê p ú = [ 0] n + n ´1 (14)
êë j ( P ) úû ( C B)
length of the origins of the two JCS's.
j13
K
= S T S - L2 = 0 (12)
S = D - A i ri (13) Table 4: Combining basic Tr. and Rot. restrictions
on ¬ Spherical Revolute Rigid
Eq. (12) might be used for creating a rigid rod joint or a
¯ Rot 3 Rot 1 Rot 0
reciprocating time-dependent kinematical actuator depending
upon the assignment of L .
Point
Tr 0
Creating kinematical Joints Using Basic constraints
The present paper idiomatically uses the words spherical,
revolute and rigid, and point, slot, plane and space to Slot
respectively indicate that three, one and zero relative rotational Tr 1
DoF, and zero, one, two and three relative translational DoF are
assigned to be free rather than a constant value or a time- Plane
dependent function. Constraints restricting relative translational Tr 2
and rotational movements are given by Table 1 and Table 2. As
arisen in Table 3, the words rod and Jack indicate a constant
and a time-dependent distance, respectively. Space
Tr 3

Table 1: Basic constraints restricting translational DoF

Translational Situation Basic Constraints Real DoF of Systems
Point Determined translation j1K = 0, jK2 = 0, jK3 = 0 As a matter of fact, the DoF suggested by the famous
Free on slot X jK2 = 0, jK3 = 0 Gruebler's mobility criterion is mathematically the number of
Slot j1K = 0, jK3 = 0 the generalized coordinates minus the number of the constraints
Free on slot Y
defining each joint without regard to the fact that the whole
Free on slot Z j1K = 0, jK2 = 0 joints constraints must be independent. MASS numerically
Free on plane YZ j1K = 0 obtains the number of the independent joints constraints of the
Plane Free on plane ZX jK2 = 0 systems through computing the rank of Jacobian matrix, i.e.
é¶jC ¶R ¶jC ¶P ù , of the whole joints constraints. As a result,
Free on plane XY j3K = 0 ëê ú û

Space Free translation - MASS numerically computes the real DoF of the systems.

Table 2: Basic constraints restricting rotational DoF Time-dependent Joints

Rotational Situation Basic Constraints Time-dependent joints are kinematical actuators in fact.
Spherical Free rotation - Considering Figure 1, position and orientation of JCS J j
Free about X axis jK4 = 0, j5K = 0, jK6 = 0 relative to JCS J i are determined by the independent joint
Revolute Free about Y axis j = 0, j = 0, j = 0
K
7
K
8
K
9
variables, i.e. X (t ) , Y (t ) , Z (t ) , Q(t ) , F (t) and Y (t ) . When
Free about Z axis j = 0, j = 0, j = 0
K K K they are assigned as functions of time, time-dependent joints are
10 11 12
created. In inverse dynamics problems the system whole
Rigid Determined rotation jK4 = 0, jK5 = 0, ..., j12
K
=0
motions are governed by time-dependent joints and it is
required to find the constraint forces and moments.
Table 3: Basic constraint restricting distance between points
Distance Situation Basic Constraint The First Time Derivatives of Constraints
Rod / Jack Assigned distance j13
K
=0 The first derivative of the basic kinematical constraint Eq.
(8) with respect to time is given by Eq. (15).
As shown by Table 4, various types of kinematical joints
might be constructed in MASS by the combination of the basic
translational and rotational constraints illustrated by Table 1

4 Copyright © 2009 by ASME

 é Vi ù é ¶jC ¶jC ù
ê ú ê ú
¶P ú é V ù é n ù
C
é é ¶j ù é ¶j ù é ¶j ù é ¶j ù ù ê Vj ú ê ¶R
K K K K
êê ú ê ú ê ú ê úú ê ú = n
K
(15) ê ú = ê ú (22)
êë ë ¶R i û w ê ¶jp ú êë V úû êë n p úû
p
ëê ¶R j úû ë ¶P i û ¶P
ëê j ûú úû ê i ú ê éë0 nB ´3n B ùû ú
êë wj úû ë ¶P û
Eq. (15) has been reformed into Eq. (16) by replacing wi
The Second Time Derivatives of Constraints
¶jK ¶jK
and with 2G i Vip and 2 G i respectively. Considering the fact that the angular accelerations w& i is
¶Pi ¶P i
equal to 2G i V& ip , the second derivative of the basic kinematical
é Vi ù
constraint Eq. (8) with respect to time is derived and shown by
é é ¶j ù é ¶j ù é ¶j ù é ¶j ù ù ê Vj ú
K K K K
Eq. (23).
êê ê ú K
ú ê ú ê ú ê úú p = n (16)
êë ë ¶R i û ëê ¶R j ûú ë ¶Pi û ëê ¶Pj ûú úû êê Vi úú éV& ù
i
p ê& ú
êë Vj úû é é ¶j ù é ¶j ù é ¶j ù é ¶j ù ù ê Vj ú
K K K K
êê ú ê ú ê ú ê úú ê & p ú = g
K
(23)
The Jacobian of the thirteen basic constraints in Eq. (16), is ¶R ¶ R ¶P ¶P
ëê ë i û ëê j ûú ë i û ëê j ûú ûú ê i ú V
given by Eq. (17) within which E1 , E 2 and E3 are the columns & pú
êë V j û
of matrix A J j [1].
The right hand side of Eq. (23) is shown by Exp. (24).
é é ¶jK ù é ¶jK ù é ¶jK ù é ¶jK ù ù g K = éë g1K g 2K L g13 K
ùû
T
(24)
êê ú ê ú ê ú ê úú =
êë ë ¶R i û êë ¶R j úû ë ¶Pi û êë ¶Pj úû úû The components of Exp. (24) are given by Eq.'s (25)-(27).
² é g1 ù é X
K && ù
é A T é é -A T ù éA T ù é 2A T ù é -2A Ti A j r%jG j ù ù ù
ê Ji ê ë
ë i û ë i û ë i D Gi û ë û ûú ú ê K ú ê && ú
ê ú êg 2 ú = ê Y ú -
ê A TJi é[ 03´3 ] [ 03´3 ] é 2A ²T
AEG ù °G ù ù ú
é -2A iT A j E && ú
ê g 3K ú ê Z (25)
êë ë i j 1 iû ë 1 j ûú
ûú (17) ë û ë û
ê
ê T é é ²
ê A Ji êë[ 03´3 ] [ 03´3 ] ë 2A i A jE 2G i û
T ù ±
é -2A i A j E 2G j ù ù ú
ë
T
û úû
ú A TJi é-2w
ë i i j (
± AT V - V + A w
i
± ±± T ) ±± ù
j j rj + wi wi A i D + A i A j wj wj rj
T
û
ê ú
ê A T é[ 0 ] [ 0 ] é 2A ² T
AEG ù ±G ù ù ú
é -2A iT A j E
ê J i ê 3´ 3
ë 3´3
ë i j 3 iû ë 3 j ûú
ûú é g K4 g K5g K6 ù é D && &&
D D&& ù
11 12 13
ê ú ê K Kú ê && && && ú -
êë 2 S éë[ -I 3´3 ] [ I 3´3 ] éë 2A i r%i G i ùû é - 2A j r%jG j ù ù ú ê g7 g 8Kg 9 ú = êD 21 D D
T
23 ú
ë ûû û 22
ê g10
K
g11
K
g12K ú êD&& &&
D D&& ú (26)
The right hand side of Eq. (16) is given by Exp. (18). ë û ë 31 32 33 û

nK = éën1K n K2 L n13
T
± ± ± ± ±w
A Ji éëwi wi A i A j - 2 wi A i A j wj + A iT A j w
T T T ±ù
ùû j j û AJj
K
(18)
The components of Exp. (18) are given by Eq.'s (19)-(20).
én1K ù é X& ù é n K4 n K5 n K6 ù é D &
11
&
D 12
& ù
D 13 g13
K
( ) ±±
&& + L& 2 + 2 S T A w
= 2 LL ( ±±
i i wi ri - A j w j w j rj ) -
ê Kú ê&ú ê K ú ê & & & ú (27)
n = Y , n n K
n K
= D D D
ê ú ê ú ê ú ê 23 ú (19)
( ± r +A w±r ) ( -V + V - A w± r + A w± r )
T
2 7 8 9 21 22
2 - Vi + Vj - A i w
ên 3K ú ê Z& ú ên10 K ú ê & & & ú i i j j j i j i i i j j j

ë û ë û ë
K
n K
11 n12 û D
ë 31 D 32 D 33 û
Eq. (28) is required to be substituted in Eq. (26). In Eq.
n13 = 2 LL&
K
(20) & J is the angular acceleration of JCS J j relative to JCS J i .
(28), w
Eq. (21) is required to be substituted in Eq.'s (19). In Eq. && && && ù é D
éD D D D12 D13 ù
(21), wJ is the angular velocity of JCS J j relative to JCS J i . ê &&
11 12 13
ú
11
± ±J w
±J ù
D
ê 21
&&
D D 23 ú = êêD 21
&& D 22 D 23 úú é w
&J +w (28)
&
éD & & ù éD
22
ë û
D D D12 D13 ù êD&& &&
D && ú ê D
D D 32 D 33 ûú
ê&
11 12 13
& ú = êD
11
±J ë 31 33 û ë 31
& D 23 úú w
32
D
ê 21 D D 23 ú ê 21 D (21)
22 22
Considering Eq.'s (4)-(5), acceleration of the whole
&
êD D& D& ú êD D D ú
ë 31 32 33 û ë 31 32 33 û constraints of a system is derived from differentiation of Eq.
Considering Eq.'s (4)-(5), velocity of the whole constraints (22) with respect to time and is proposed by Eq. (29).
of a system is derived from differentiation of Eq. (14) with é ¶jC ¶jC ù
ê ú &
respect to time. It is given by Eq. (22). ¶P ú é V ù é g ù
C
ê ¶R ê ú =ê ú (29)
ê & p ú ê gp ú
¶jp ú êë V
é ù
ê ë 0 n B ´3 n B û ú û ë û
ë ¶P û

5 Copyright © 2009 by ASME

SPRINGS, DAMPERS AND ACTUATORS
Similar to the joints, reciprocating and rotational elements Ji ë (
ìM Spring = A é 0 0 +K T atan [d , d ] - T Free ù T
ï i 21 11 )û
í (33)
of spring-damper-actuator (SDA) might be created in MASS.
( )
T

As an example, a reciprocating and a rotational element are

ïM Spring
î j
= A Jj ë
é 0 0 - K T
atan [ d 21 , d 11 ] - T Free
ù
û
illustrated by Figure 3 within MASS environment. ìM i Damper
= C éëA i A j wj - wi ùû
T T

In Eq. (30), the vector ¡ is expressed in GCS and ï

í Damper (34)
represents the position of J j relative to J i . ïîM j = C T éëA Tj A i wi - wj ùû
¡ = R j - R i + A j rj - A i ri (30) ìM Actuator = m Actuator A Ji [ 0 0 -1]
T
ï i
í Actuator (35)
= m Actuator A J j [ 0 0 +1]
T
ïîM j
Similar to Eq. (32), the total moment applied to each body
can be readily obtained by adding up the moments in Eq.'s (33)-
(35).

DYNAMICS
The Newton-Euler motion equations of a constrained
multibody system are represented by Eq.'s (36). It is obvious
Figure 3: Reciprocating and rotational SDA that the force and moment Eq.'s (36) are respectively expressed
in GCS and CCS of each body.
Reciprocating SDA MV & = fA + f C , I w& +w
% I w = mA + mC (36)
As a matter of fact, a reciprocating spring-damper-actuator
The Eq.'s (37)-(38) are obviously derived from Eq.'s (36)
is a rod joint that exerts two opposite equimagnitude forces on
within which the constraint forces and moments are expressed
the matched bodies. The forces are parallel with ¡ shown by
in terms of Jacobian and Lagrange multipliers corresponding to
Eq. (30).
the constraints. It should be noted that the Jacobian and
The reciprocating spring, damper, and actuator forces
Lagrange multipliers corresponding to the quaternions
exerted on point J i of the ith body are given by Eq.'s (31). They constraints have been also considered. Eq. (38) has been
are expressed in GCS. expressed in terms of quaternions.
ì Spring é LFree ù T
& + éê ¶j ùú l C = f A
C
ïFi = K L ê1 - ú¡ MV (37)
ï êë ¡ úû êë ¶R úû
ï
CL
ï Damper
(
= T é Vj - Vi + A j w ±r - A w )
± T ¡ù ¡ T T
é pù
& p + éê ¶j ùú l C + ê ¶j ú l p = 2G T ém A - w
C
íFi i i ri (31) °Iwù
ï ( ¡ ¡) ë ê j j úû é 4G T IG ù V
ë û
ëê ¶P ûú êë ¶P ûú ë û (38)
ï Actuator
ïFActuator = -f ¡
The method of solution is determined by the problem type.
ï i
¡ In inverse dynamics problems, algebraic equations and in
î
forward dynamics problems ordinary differential equations are
Resultant force exerted on points J i and J j of the ith and jth
to be solved. It is due to the fact that kinematics is decoupled
bodies are given by Eq. (32). from kinetics in inverse dynamics problems, on the contrary of
ìïFiSDA = FiSpring + FiDamper + FiActuator forward dynamics problems.
í SDA (32) An advanced numerical technique has been used in MASS
ïîFj = -FiSDA for solving holonomic algebraic differential Eq. (39) composing
The above forces in Eq. (32) produce two moments about of Eq.'s (29), (37)-(38). The technique iteratively modifies the
the mass centers of the matched bodies. Their projections onto results of explicit numerical integration in order for the
their CCS's are ri ´ A iT FiSDA and rj ´ A Tj FjSDA . constraints Eq.'s (14) and (22) to be complied with. The
modification can be interpreted as kinematical assembling
during integration. It is optional and might be activated. It is
Rotational SDA
essential in inverse dynamics problems, whereas it is sometimes
As a matter of fact, a rotational spring-damper-actuator is a
constructive and sometimes destructive in forward dynamics
revolute joint that exerts two opposite equimagnitude moments
problems. In inverse dynamics problems the number of
on the matched bodies. The moments are about the axis of the
constraints and coordinates are the same. In forward dynamics
revolute joint.
problems, however, the number of constraints is less than that of
The rotational spring, damper, and actuator moments
the coordinates, so MASS uses the optimization of SVD
exerted on the ith and jth bodies are given by Eq.'s (33)-(35).
They are expressed in their CCS's.

6 Copyright © 2009 by ASME

method after each step of integration to satisfy constraints and numerical integration step in forward dynamics problems. Eq.
their velocities to prevent cumulative errors violating joints. (39) must be written in state space in order to be numerically
é é ¶jC ù
T ù integrated. MASS provides users to adjust simulation
ê M 0 ê ú 0 ú parameters of numerical integration methods.
ê ëê ¶R ûú ú
ê ú
& é fA ù
é ¶jp ù ú é V ù ê
T T
ê é ¶jC ù Handling Impulsive Loads
ê 0 4G T IB G ú ê ú ú
ê ú ê ú & T é A ° ù Initial impulsive loads might be considered in MASS via
êë ¶P úû ú ê V ú = êê 2G ëm - wIB wû úú
p
ê êë ¶P úû
ê úê Cú Eq. (40). Using the initial linear and angular impulsive loads, i.e
ê é ¶jC ù é ¶jC ù ú êl ú ê gC ú
êê ú ê ú 0 0 úê p ú ê ú ILi and IiT , Eq. (40) changes the initial velocities and
ê ëê ¶R ûú ëê ¶P ûú ú ë l û êë
p
g úû
simultaneously gives the initial impulses of the Lagrange's
ê ú
ê 0 é ¶jp ù multipliers.
ê ú 0 0 ú
ê êë ¶P úû ú Eq. (40) has been obtained by integrating Eq. (39) with
ë û
respect to time. It computes velocities immediately after
(39) uuuuur
applying impulsive loads. In Eq. (40), C i z i is expressed in CCS
Kinematical Assembling of Systems of the ith body.
Before simulation starts, the system must be assembled, é T ù
é ¶jC ù
that is, the initial conditions, i.e. positions and velocities, of the ê M 0 ê ú 0 ú
bodies must be modified in order to satisfy the constraints and ê ëê ¶R ûú ú
ê ú
p ù ú é V (e) ù
their velocities. MASS employs Newton-Raphson method to ê T T
solve the nonlinear system of constraint Eq. (14) and utilizes é ¶j C ù é ¶j ê ú
ê 0 4G T I B G ê ú ê ú ú ê V p ( e )ú
Gaussian elimination process to solve the linear system of ê ëê ¶P ûú ëê ¶P ûú êú
velocity Eq. (22). A complete pivoting technique is used to skip ê ú C
ú= (40)
redundant equations. The initial guess for positions can be ê é ¶j ù é ¶j ù
C C
úê l ú
êê ú ê ú 0 0 úê p ú
prescribed by user in order to make the process converge. If the ê êë ¶R úû êë ¶P úû ú êë l úû
number of equations is less than the number of unknowns, user ê ú
ê 0 é ¶jp ù
has option to mark and determine some positions and velocities ê ú 0 0 ú
in order to reach a desirable initial configuration in terms of ê ê ¶P ú ú
ë ë û û
positions and velocities of the bodies. If user does not mark
é é 03( i -1)´1 ù ù
enough items to compensate the shortage of equations, MASS ê nL ê ú
ú
employs SVD method to calculate an optimum configuration
that fulfills the constraint equations.
ê
ê å ê ILi
ê
i =1 0
ú + MV ( 0 )
ú
ú
ú
ê ê 3( n B - i )´1 ú ú
ê ë û ú
Inverse Dynamics Problems ê ú
é é 03( j -1)´1 ù é 03( i -1)´1 ù ù
In inverse dynamics problems, motion is completely ê ê nT ê ú
ú n L ê uuuuur úú
prescribed and eventually positions, velocities and accelerations
are found through solving some nonlinear and linear algebraic ê å
ê 2G T ê ê I T
ê j =1 ê j
ú+ å
ú i =1 ê
êC i z i ´ A i Ii ú + 4G I B GV ( 0 ) ú
T L ú
úú
T p
ú
equations of kinematics without regard to kinetics. Motion Eq.'s ê ê ê03( n - j )´1 ú ê 03( n -i )´1 ú ú ú
ê ë ë B û ë B
ûû ú
(36) are algebraic and linear with respect to the constraint loads ê C C ú
and the applied loads required for the prescribed motion. In the ê ¶j ¶j ú
V ( 0) + Vp ( 0)
last stage of solving an inverse dynamics problem the constraint ê ¶R ¶P ú
loads are achieved via solving Eq.'s (36). ê p ú
ê ¶j ú
êë V ( 0)
p
úû
¶P
Forward Dynamics Problems
In forward dynamics problems, kinematics and kinetics are
coupled and eventually the positions and velocities are achieved Handling Rigid Collisions
through solving some ODE's (Ordinary Differential Equations). MASS is capable of modeling collisions among sphere-
As well as the present article, where the generalized coordinates shaped rigid bodies. At each time-step of numerical integration,
are more than DoF the mentioned ODE's are coupled to some detection of probable collision happenings might be activated.
linear algebraic equations of constraints acceleration. It is Considering Figure 4, if a collision has occurred, the current
obvious that the algebraic differential motion Eq. (39) is time-step is ignored and the current numerical process is
algebraic and linear with respect to the accelerations, repeated with a new time-step being equal to half of the current
Lagrange's multipliers and eventually the constraint loads. time-step. When this cycle is stopped, the exact instant of
Accelerations and Lagrange's multipliers are simultaneously collisions has been captured and the maximum amount of
found through solving Eq. (39) in the first stage of each

7 Copyright © 2009 by ASME

interference due to all collisions at the instant has become less é
T
é01´3( n -k ) ù ùú
T T
Ckf = ê éë 01´3(k1 -1) ùû éCkf 1 ù é01´3( k -k -1) ù éCkf ù
ë û ë 2û ë 2 û
(46)
than a prescribed tolerance. ë ë û 2 1 B
û
C kf 1 = mk t k1 - n k1 , C kf 2 = mk t k 2 - n k 2 (47)

C m = éC1m Cm L Cm ù (48)
ë 2 nh û
T
éé ù é mù
T
é01´4( k -k -1) ù éC km ù
T
é01´4( n -k ) ù ùú
Cm
k = ê ë 01´ 4(k1 -1) û ë C k1 û ë û ë 2û ë B 2 û
(49)
ë û
2 1

uuuuuuur
ì C m = 2G T é C h ´ A T C f ù
Figure 4: Contact points on body surfaces ï k1 k 1 ë k1 k1 k1 k1 û
í m uuuuuuuur (50)
ïC k 2 = 2GkT2 éë Ck 2 hk 2 ´ A kT2 C kf 2 ùû
User should specify the friction and restitution coefficients î
between colliding bodies in MASS. When a number of
collisions are identified and captured by MASS, velocities FEATURES AND CAPABILITIES OF MASS
immediately after collision are found via the set of linear Eq.'s In this section, features and capabilities of MASS are
(41). Simultaneously, the impulses of Lagrange's multipliers and demonstrated via explaining the required steps for modeling
the impulses of normal colliding forces and eventually the and simulation of a multibody system in MASS.
impulses of the corresponding friction forces are computed.
é T ù Modeling a Four-Bar Linkage
é ¶jC ù Planar systems are three-dimensionally modeled and
ê M 0 ê ú 0 -C úf
ê êë ¶R úû ú numerically simulated by MASS. A 2D/3D linkage can be
ê ú é V ( t + e) ù
ê T T úê C
ú defined in MASS through the following steps:
ê 0 é ¶jC ù é ¶jp ù m ú ê Vp ( t + e )ú
T
4G IB G ê ú ê ú -C ú C Step 1: Creating simple shapes or importing complicated
ê ê ú
ê ëê ¶P ûú êë ¶P úû úê shapes from AutoCAD.
lC ú
ê C úê ú Step 2: Attaching JCS's to the bodies and locating them
ê ¶j ¶jC
0 0 0 úê lp ú properly so that they can be used in defining joints.
ê ¶R ¶P úê ú
ê úê N úû
ê 0 ¶jp ë
0 0 0 ú
ê ¶P ú (41)
ê ú
ë a 0 0 0 0 û
é MV ( t C ) ù
ê ú
ê 4G IB GV ( t C )
T p
ú
ê C C ú
ê ¶j ¶j ú
= ê ¶R V ( t C ) + ¶P V ( t C ) ú
p

ê ú
ê ¶jp p ú
ê V ( C)
t ú
ê ¶ P ú
ëê b úû
For the kth collision, the collided bodies, i.e. the k1th and Figure 5: Connecting bodies by joints
k th
2 bodies shown by Figure 4, have normal and tangential unit
uuuuuuuur Step 3: Creating joints by selecting the attached JCS's and
vectors n k1 , n k 2 , t k1 and t k 2 . The vectors C k1 hk1 and specifying the type of joints via determining each joint variable
uuuuuuuur as a free, fixed or a function of time.
C k 2 hk 2 are expressed in CCS's of the collided bodies. Exp.'s
Step 4: In order to reach the desired initial configuration
(42)-(50) are required for Eq. (41). and velocities, one should emphasize and mark some IC's
T
a = é a1T a T2 L a Tnh ù (42) (initial conditions) as far as the mechanism DoF permits. In the
ë û example shown by Figure 5, the revolute joint between ground
ak = éê éë01´3(k1 -1) ùû é n Tk1 ù éë01´3(k 2 - k1 -1) ùû é n kT2 ù éë01´3( n B - k 2 ) ùû ùú (43) and the 1st link as shown by Figure 6, and its angle and angular
ë ë û ë û û velocity are emphasized in the dialogue box shown by Figure 7.
T
b = éb1 b2 L bnh ù , bk = -sk a k V ( t C ) (44) This can be treated as one of the novelties of MASS in
ë û comparison with other application tools.
C f = é C1f C f2 L C fn h ù (45) The positions and velocities of the other bodies must be
ë û calculated via assembling process for having a set of initial
conditions consistent with constraints and their velocities.

8 Copyright © 2009 by ASME

 Figure 9: Kinematical assembling of the mechanism

Step 7: Simulating the model and observing the motion of

the mechanism during solving process.
Step 8: As shown by Figure 10, observing motion tracks,
vectors and diagrams of positions, velocities, accelerations and
constraint loads with respect to time.
Figure 6: Defining joints Figure 7: Emphasizing IC's

Step 5: In the dialogue box shown by Figure 8, one might

adjust the simulation settings such as time-step, integration
method and tolerances. MASS is provided with constant and
adaptive time-step explicit numerical integration methods.
User can optionally mark choices to satisfy constraints and
their velocities after each step of integration to prevent
cumulative errors violating joints.

Figure 10: Simulating a forward dynamics problem

It will be a forward dynamics problem if all four joints are

defined as revolute; but if one of the joint variables, e.g. F , in
the joint between ground and the 1st link is defined as a function
of time, it will be an inverse dynamics one.

VERIFICATIONS
In this section, simulation results of two selected problems
are verified via comparing with the corresponding results of
Working Model 4D ver. 6. Models built in the two applications
are completely similar in all aspects.
Figure 8: Simulation settings
Case I: Spatial Pendulum
Step 6: As shown by Figure 9, kinematical assembling of The spatial pendulum shown in the Figure 11 consists of
the mechanism in order to find positions and velocities three revolute joints about three different axes defined between
satisfying constraints and their velocities. The Newton-Raphson three similar boxes and ground. Time history diagrams of
process might be graphically traced for educational purposes. position and orientation of Box 3 as well as constraint forces
The user-defined conditions are used as initial guess for starting and moments applied to Box 1 due to its revolute joint with the
the process. The emphasized conditions are assumed as known ground are chosen to be compared in Figure 12 (a-h). The
values. results of MASS in black and thin can be readily distinguished
from the results of Working Model in green and thick. All

9 Copyright © 2009 by ASME

diagrams are coincident and the deviations in orientation neglected. Initial velocities are such that the mass center of the
diagrams are due to the arc conventions. Positions, angles, body remains stationary. Initial pose (position and orientation)
forces and moments are respectively in terms of m, deg, N and locates the edges of the body coincident with the axes of the
N-m. Simulation time is 10 seconds. MASS and Working GCS. It is obvious that the principal moments of inertia about
Model solving times are respectively 59 and 12 seconds. the mass center of the body might be calculated as follows:
ìï Imax = 0.002378318028 Kg.m 2
í 2
(51)
ïî Imin = Iint = 0.001725404482 Kg.m
The principal direction corresponding to I max might be
calculated as follows:
1
u max = [1 1 1]T (52)
3
It should be noted that, any axis perpendicular to u max
passing through the mass center is a principal axis. So, u t
might be considered as a principal direction.
1
ut = [ 0 -1 1]T (53)
2
Figure 11: A 3D pendulum with three revolute joints Motions of the body due to the initial angular velocities
T T
éë0.5 0.5 0.5ùû , éë0 -0.5 0.5ùû and éë0 0 0.5ùû T rad/s are
respectively shown by Figure 13 (a-c) derived from MASS, and
by Figure 14 (a-c) derived from Working Model. Some JCS's
have been attached to the body and their tracks have been
plotted to illustrate motion. Comparing Figure 13 (a-c)
respectively with Figure 14 (a-c) bears witness that the results
(a) Box 3 – Position – X (d) Box 3 – Orientation – q
of the two applications are not the same for a unique problem.

(b) Box 3 – Position – Y (e) Box 3 – Orientation – f

Simulation Time: 2.4 s Simulation Time: 8.9 s Simulation Time: 23.37 s

Solving Time: 0.64 s Solving Time: 2.4 s Solving Time: 6.2 s
(c) Box 3 – Position – Z (f) Box 3 – Orientation – y
(a) (b) (c)
Figure 13: Results from MASS

r r
(g) Constraint Force F (h) Constraint Moment M
Figure 12: Comparison between MASS and WM 4D results

Case II: Torque Free Motion

Torque free motion of a rigid body with homogeneous mass
distribution having a density of 1.0 Kg/m3 has been simulated. Simulation Time: 2.4 s Simulation Time: 8.9 s Simulation Time: 23.37 s
The shape of the body, as shown by the figures, is an octant of a Solving Time: 0.6 s Solving Time: 2.3 s Solving Time: 5.8 s
sphere having a radius of 0.5 meters. The gravity and the (a) (b) (c)
applied torques about the mass center of the body have been Figure 14: Incorrect results of Working Model 4D

10 Copyright © 2009 by ASME

 Unlike the initial angular velocity ëé0 0 0.5ûù T , the initial ► MASS, as a C++ application, has been established on
T T
top of the mentioned dynamic model. It has successfully
angular velocities éë0.5 0.5 0.5ùû and éë0 -0.5 0.5ùû are simulated several systems, some of which has been exposed in
respectively about u max and u t that are the principal the section of "Verifications" within which the results have been
directions of the rotary inertia tensor of the body. Principal axes compared with that of Working Model. The results of the
of the rotary inertia tensor about the mass center of a body are comparison have been also discussed in the mentioned section.
the stable or unstable equilibrium axes for the body spatial ► MASS might solve forward dynamics problems using
rotation. Therefore, unlike Figure 14 (a-c) derived from the adaptive time-step Runge-Kutta-Fehlberg 5-6 or 7-8
Working Model, the Figure 13 (a-c) derived from MASS is numerical integration methods as well as some simpler
acceptable. procedures. Computing consistent initial conditions for a
Figure 14 (a-c) has been derived from Working Model simulation manually is as difficult as simulating the problem by
within which the body has a uniform density as shown by Figure hand without computer. SVD method has been found to be
16 that involves the unexpected inertia tensor with respect to automatically capable of providing users with consistent initial
the principal frame. In Working Model, as shown by Figure 17, conditions having independent and adequate components. To
one might select custom density and might manually calculate assemble the system during integration has been observed to be
and enter the expected components of the inertia tensor with successful in eradicating cumulative integration errors that
respect to the CCS that originates from the mass center. In this violates the definitions of the joints.
case, the acceptable results coinciding with Figure 13 (a-c) are ► The quantity of kinetic, potential and total mechanical
achieved and are depicted by Figure 15 (a-c). energies, constraints' position and velocity errors and DoF are
shown in the status bar below the main graphic view of MASS,
as shown by Figure 10. These quantities are a few necessary,
but not sufficient, check points that convince users about the
validity of the numerical results. One might monitor these
quantities and be sure about the fact that the numerical error is
not high enough to break the energy conservation law, for
instance, in a conservative or non-conservative system being
simulated by MASS. Real DoF of a system numerically
evaluated by MASS is a demanded engineering parameter for
Simulation Time: 2.4 s Simulation Time: 8.9 s Simulation Time: 23.37 s some mechanical systems such as spatial or planar mechanisms.
Solving Time: 0.6 s Solving Time: 2.3 s Solving Time: 5.8 s The evaluations are too long or actually impractical to be
(a) (b) (c) calculated by hand without a computer program. The
Figure 15: Correct Results from Working Model 4D importance of this DoF might be found when the Gruebler's
mobility criterion offers a number lower than the real DoF of a
linkage, whereas one might easily imagine the real DoF thereof.
It might be verified that the real DoF of inverse and forward
dynamics problems are respectively zero and a positive number.

REFERENCES
[1] Khorsandijou, S.M., and Shokouhfar, S., Computer Aided
Dynamic Modeling of Spatial Mechanisms (in Farsi),
Figure 16: Uniform density Figure 17: Custom density Ayeh Cultural Publication Institute, Tehran, Iran, 2004,
ISBN: 964-5741-03-3, Along with the CD of MASS.
CONCLUSIONS [2] Nikravesh, P., Computer Aided Analysis of Mechanical
► Unified dynamic model of 3D multibody systems has Systems, Prentice Hall, Englewood Cliffs, 1988.
been achieved based on coordinate viewpoint and thirteen basic [3] Haug, E.J., Computer Aided Kinematics and Dynamics of
constraints Eq.'s (9), (11)-(12). It includes the algebraic Mechanical Systems, Allyn and Bacon, Boston, 1989.
differential motion Eq. (39), impulse and collision governing [4] Shabana, A.A., Dynamics of Multibody Systems,
Eq.'s (40)-(41). Cambridge University Press, Cambridge, 2005.
► To observe these equations in detail, reveals that [5] Shabana, A.A., Computational Dynamics, John Wiley &
Jacobian and g in the constraints' velocity and acceleration Sons, Inc., New York, 2001.
Eq.'s (17), (25)-(27) have been exposed in quite fresh forms. [6] Haddadpour, H., Analysis of Kinematics and Dynamics of
They have the highest possible generality and quality of being 3D Mechanisms Using Constraint Equations Method,
closed-form. M.Sc. Thesis, Tehran University, 1995.
[7] Press, W.H., Numerical Recipes in C (The Art of Scientific
Computing), Cambridge University Press, 1992.

11 Copyright © 2009 by ASME

View publication stats