Computer Methods in Applied Mechanics and Engineering North-Holland
Dynamical
equations
91 (1991) 1403-1414
in natural coordinates
Hans Brauchli and Ren6 Weber Institute of Mechanics, Swiss Federal Institute of Technology, Ziirich, Switzerland, Industries. Center for New Technologies, Ziirich, Switzerland
Received 31 October 1990 Revised manuscript received 13 February
and Oerlikon
1991
A simple set of equations for the description of a general class of multibody systems is presented. They are as well suited for simulation as for developping the dynamical theory. Applications to inverse dynamics and stability theory are given. The Delta robot of Clavel has been chosen as an example for simulation.
1. Introduction
Methods for the analysis and simulation of multibody systems are of increasing consequence in space research, robotics, vehicle dynamics and other engineering applications. Most existing methods have two major shortcomings. The resulting equations, often arrived at by computer manipulations, are useful for simulation but too complicated to be of value for analysis. For example the derivation of the equations of the first variation, necessary for stability analysis, is practically unfeasible. Furthermore, since redundant coordinates have to be introduced for the description of the system, the separation of the resulting equations into the algebraic equations for the reactive forces and the proper differential equations of motion is difficult. We present here a simple set of equations for the description of a general class of systems of rigid bodies which does not suffer from these drawbacks. Our method depends on two main features. First, the natural metric of the configuration space, i.e. the mass-metric, is respected. A consistent separation of the equations for the reactive forces and the differential equations for the constrained system is then obtained. We also introduce the concept of affine body, which allows to describe the system with holonomic velocity parameters while the mass matrix of the free system is constant. The equations obtained are quite simple. They are well suited for simulation and as matrix equations they are adapted to the needs of vector computers. They may also serve as a basis for further analysis. Two examples for this are given below: inverse dynamics and stability analysis. 2. Dynamics
in natural coordinates
2.1. Mass-orthogonal
projections
A projection method based on the classical equations of mechanics for the description of scleronomic multibody systems has been proposed previously [l-4]. The constraints of the 00457825/91/$03.50
0
1991 Elsevier Science Publishers
B.V. All rights reserved
1404
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
system are written as kinematical constrains and are used to construct two complementary projections (Yand /3 splitting the velocity space of the free system into the spaces of admissible and inadmissible velocities, respectively. They are uniquely determined by a condition of mass-orthogonality, resulting in an appropriate splitting of momentum space. Let u denote the generalized velocities. The kinematical constraints are written in matrix form as E’u = 0,
(1)
where E’ denotes an ~1,by II matrix depending on the coordinates, iz the degree of freedom of the free system and II, the number of the constraints. Let M denote the mass matrix of the free system and A its inverse. The conditions of mass orthogonality are
(Ma)’ = Ma With the fundamental
(MP)’ = Mp .
)
(2)
matrix
G=E’AE,
(3)
which is a positive definite n, by n, matrix, we solve the system of equations GX’ = E’
(4)
and find the two projections cY=I-A7T,
p=Afl
(5)
with m=n’=XE’.
(6)
The mass matrix of the constrained
system is then
M,,=dMa=M-m,
(7)
and the matrix A ao =aAa’=A-A?rA acts as its generalized a(u)
(8)
inverse. Let
= EG-‘fit
=Xl?'
.
Then the reactive forces as functions of position and velocities can be determined equation
0, = -cqu)u
+ p’niu
)
(9) from the
(10)
140.5
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
and the differential
equations
@ = 0, - MU
of motion of the system read
+ p’niu )
u=Ap,
w,
(12)
where the generalized forces $, and 0, may include correction terms for using nonholonomic velocity parameters as well as the derivative of the kinetic energy as a function of the velocities and coordinates with respect to the coordinates. 2.2. Af$ne body and natural coordinates Let us now, following [5], introduce the concept of an afine body. Consider a rigid body in space and select four material points in general position as nodal points. An arbitrary point of the body can then be represented by barycentric coordinates. If the four nodal points move independently, the body will perform an affine transformation. The corresponding mass matrix is constant. The rigid body can be obtained by imposing constraints. If this method is used, the physical interpretation of the constraint forces will have to be modified. The equations of motion, however, are considerably simplified. The correction terms in the generalized force $, disappear, and the velocities can now be identified with the time derivatives of the coordinates. The equations of motion now read
ri = Q,
- JW)u ,
Q=u=Ap,
(13) (14)
or, equivalently,
3. Inverse dynamics 3.1. Inverse kinematics We consider systems such as the Delta robot without redundant degrees of freedom. This means that the motion of the interesting part, e.g. the hand of the robot, completely determines the motion of the whole system. Then, formulating the theorem of projected velocities for two points P and Q,
(16)
v;=vl,, we obtain a linear system with coefficient matrix C for the missing velocity parameters. also formulate a theorem of projected accelerations,
We
(17) where a correction term involving rotational We have the following theorems.
velocities & orthogonal
to the line PQ appears.
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
1406
THEOREM 1. In natural coordinates the rotational velocities A orthogonal to the line PQ can be computed from the velocities of the points P and Q by
4’ = (uQ -
u,)“/I’,,
.
(18)
THEOREM 2. The coefjicient matrix C of the linear system for the missing velocities and accelerations are the same.
Therefore, tiation.
the accelerations
can easily be computed
directly without the need for differen-
3.2. Inverse kinetics We apply the principle of virtual work to the equation
of motion in the form
Mcj = Q, - fl(tj)tj
(19)
Let qo(t) be the desired motion and let 6q, denote an admissible virtual displacement. we have Theorem 3. THEOREM
3. The driving force for the desired motion can be determined
by
6W= Ql, 6q, = &,M 6q, . PROOF.
Then
(20)
From atJ2 = atxjj’
= atEG-l&’
=
0,
(21)
we find $)a’
6q, =
fjgva 6q,
=
0.
(22)
0
Theorem 3 implies that for the determination of the driving forces in natural coordinates the mass matrix M of the free system is needed.
only
4. Curvature and stability 4.1.
Geometry of dynamics
A geometrical approach to the question of stability in mechanical systems has been formulated by Synge [6]. According to the principle of Jacobi the orbits of a conservative system correspond to the geodesics of the configuration space of the mechanical system
1407
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
endowed
with a metric. For the inertial motion this metric is given by the kinetic energy, the
mass metric
ds2 = dq’ M dq .
(23)
This is the ‘kinematical line element’ undergoes a conformal change:
of Synge. If potential
forces are present,
d? = 2(h - V(q)) ds* = 2(/z - V(q)) dq’ M dq . This is the ‘dynamical line element’ of Synge. V(q) constant. A classical result of Jacobi and Levi-Civita relates This can be used to determine the local linear stability sectional curvatures. In fact, the differential equation dqk + R;&j’
i3qh = Q;h 6qh .
Here, Rfhj denotes the curvature The exprressions
Rfhj4’$ -
denotes the potential
the metric
(24)
and h the energy
the geodesic deviation to the curvature. coefficients of the mechanical system as for the geodesic deviation reads
(25)
tensor and a comma is used to indicate a partial derivative.
QTh
(26)
can be interpreted as local linear stability coefficients. The first part is due to the geometry of the configuration space with the mass metric, the second part depends on the forces applied. This method has not been widely used due to the difficulty in computing the curvature tensor of systems of interest. 4.2. Natural embedding In our setting the configuration space of the constrained system is naturally embedded in a euclidean space, i.e. the configuration space of the free system. Comparing the equation of an inertial motion in natural coordinates, Mcj+ i@j)cj=O,
with the differential
equations
(27) of a geodesic,
one might expect the term a( 4)4 to be related to the coefficients r: of the affine connection. But, in fact, the a( 4)4 correspond to the exterior curvature of the embedding. This is further illustrated by Theorem 4. THEOREM 4. The constraint reactions for an inertial motion of the constrained system as a function of position and velocity are given by
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
1408
Qp = -%i)4.
(29)
Clearly, when using natural coordinates, the interpretation of the constrain reactions is not straightforward. One might think of the nodal points as mass points connected by strings for each constraint. Then the constraint reactions are the tensions in the strings. The intrinsic curvature of an embedded manifold is related to the extrinsic curvature of the embedding. Using an extension of Gauss theorema egregium, we find Theorem 5. THEOREM 5. The curvature R of the configuration space of the constrained system due to its mass metric as a function of four typical velocities a, b, c and d is given by R(a, b, c, d) = c’K(a, b)d - d’K(a, b)c ,
(30)
K(a, b) = ati2’(a)AJ2(b)a
(31)
where
If this is introduced 6.
.
into the differential equation for the geodesic deviation, we have Theorem
Theorem 6. The first variation of the differential equation of inertial motion of a constrained mechanical system described in natural coordinates is
Stj + Aa’~‘(Gq)Ai+j)~
- Aatf2’(4)Afl(6q)tj
= 0.
(32)
Zf forces are present, it has to be replaced by
A proof of Theorems 5 and 6 will be given elsewhere. If we consider systems with geometrical constraints quadratic in the coordinates, as will be the case for the Delta-robot, 0(u)u is given as a sum of squares of rotational velocities multiplied with the columns of X to be interpreted as a dual basis for the kinematical constraints. Since the curvature merely gives the local linear stability coefficients, it may suffice to prove instability. In the case of conservative or gyroscopic systems however, stability is more difficult to establish. A general formula has been given [7] for cases where KAM-theory is applicable. For other systems of this kind algebraic criteria cannot be found. For the engineer, therefore, it makes sense to use the local linear stability coefficients as a first indication of stability or instability and then to check by numerical integration, computing e.g. Lyapunov-exponents.
5. Examples for simulation In order to simulate a specific system with our method, we have to choose the nodal points and to specify the mass matrix M of the free system, the matrices E’ and tit describing the
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
1409
constraints and their derivatives and finally the applied forces Q. Care has to be taken to choose initial data satisfying the constraints of the system. In the examples presented below, the single bodies of the system are modelled as rods. Therefore, it is sufficient to choose two nodal points in each body. A further simplification is obtained by letting the nodal points coincide with the joints.
5.1.
Inertial motion of two linked rods
As a simple test the inertial motion of two homogeneous rods linked by a frictionless joint, moving on a frictionless plane, is considered. This system is integrable. We choose the two free endpoints as node 1 and 3 and the joint as node 2. The variables are then qt
=
[x,,
y,,
x2,
Y,,
x3,
Y31.
(34)
The mass matrix is -2 0 1 0 0 020100 ml04010 M=6 0 1 0 4 0 001020 -0 0 0 1 0 where m denotes
o-
(35)
1’ 2_
the mass of one rod. The matrices for the constraints
are
(36) and
(37) In Fig. 1 a motion is shown where the first rod is initially at rest while the second rod rotates clockwise around the joint. 5.2. A simple model for the Delta robot The Delta robot [8] is a fast robot for small manipulations. The hand of the robot has three translational degrees of freedom. It is supported by three articulated arms with parallel geometry. The hand is activated by three drive moments acting from the platform to the lower rod of the arm. In the actual system the hand is a small triangle with three axles, such that the system is four times statistically indeterminate. Due to the multiple kinematical loops as well as the statically indeterminacy the modelling of the system with conventional methods is difficult. To our knowledge, such a model has not been formulated. This is why we chose it to illustrate the efficiency of our method.
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
1410
Fig. 1. Inertial
motion of two linked rods.
A simple model for the Delta robot, disregarding the statical indeterminacy, is obtained by introducing four nodal points: one for the hand and three for the joints of the three arms. Since these joints move in vertical planes, this gives a total of n = 9 degrees of freedom for the free system with II, = 6 constraints. Hence, the degree of freedom of the constrained system is n, = IZ- II, = 3. The configuration space of the Delta robot is a closed 3-manifold, but its topology is not easily determined. As coordinates we use
where x, y, z denote the position of the hand and 5, z, e.g. the position of the first joint. The coordinates 5, r], c are the projections of the lower rods on the horizontal plane. Let mH denote the mass of the hand, m, and mL the mass of the upper and lower rods, respectively. Then the mass matrix due to the three upper rods is 1 0 0 i6 M=m,
‘b -ib 0 0 0 -6
0 0 100 0100 0
t
of
;c -$c 0 0 0
0 0 0 0 io do b0
c =
sin 60” .
-&b ;c
-;!I -6c 0
0 0 f
0 0 0 0 1 I
0
0
0
6
6
0 L 6 0 0
0 0 3 0 0
0 0 0 : 0
0 0 0 0 :
1 J
0 0 0 0
(39)
where b = cos 60” ,
For the hand mH has to be added in the first three diagonal terms, for the lower rods fm, in the last six. Hence, if all three masses are equal, the total mass matrix reads
1411
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
12 0 0 1 -b -b 0 0 0 From the geometrical (x -
0 12 0 0 c -c 0 0 0
0 0 12 04 0 0 1 1 1
1 0 0
-b c 0 0 4 0 0 0 0
0 0 0 0 0
-b -c 0 0 0 4 0 0 0
0 0 0 0 111 000 000. 000 400 040 0 0
o0
(41)
4_
constraints
5)’ + y’
+ (2 -
q)’ = L:
)
(x + bq)2 + (y - ~7)~ + (z - z2)* =
L:,, (42)
(x+ b[)* + (y + CL)* + (z - q)* = L2,, ((-a)2+2;=Lt,
(7) -
a)’ + 2; = L;
(&--a)*+zt=Lt,
and
(43)
we find the matrices for the constraints
Et=
X-5 x+bq x+bl o
Y y-q y+c[ 0
Z-2, z-z2 z-z3 0
5-X 0 0 5-a
0 _ 0
0 0
0 0
0 0
0 77+ bx - cy 0 0
7-a 0
0 0 l+bx+cy 0
z1 - z 0 0
0 z* - z 0
0
Zl
0 5-a
0 0
22
0
0 0 z3-z
0 0 z3
(44)
and
j
i-i,
&.i
0
0
i, - i
x+brj
j-oj
i-i,
0
+ + bx - cj
0
0
x+bi 0
j+c,j 0
i-i, 0
0 t
0 0
0 0
0 0
0 0
0 0
i-j E’ =
;I
0
0
.2,-i
0 i,-i 0
j-+bi+cj 0
0 21
0 0
0 j_
0 0
i, 0
-
* 0 i,_ (45)
In (44), a denotes the distance of the lower axle from the origin. Notice, that the lengths L, and L, of the upper and lower rods do not enter the equations of motion. They will implicitly appear, however, in the initial data. In Fig. 2 the motion under a dead load is shown. 5.3. 5-hinged arch The 5-hinged arch is a simple plane analogon to the delta robot. Its configuration manifold is a closed 2-surface of genus 4. We may therefore expect an interesting dynamical behaviour.
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
1412
Fig. 2. Delta
robot.
With the three middle joints as nodal points, we have the variables qf
=
[x,7
Y1,
x2,
Y2, x3-
Y31*
(46)
The mass matrix is -4 0 1 0 0 040100 m104010 M=60 10 4 0 001040 -0 0 0 1 0 where m denotes
Et
4-
Yl OY’
x2-xl x2
[
and
x1
0+ a Xl
jj”=
(47)
1’
the mass of one rod. The matrices for the constraints
Xl (1"'
=
o-
[. i,;i,
Yl0 Yl j,-j2 0
0
0 .
y:
-x3
x3
px2
00
0
0
x3-a 0
0
0
4-4
j12-J4
0
ii2
j,
ii’
y3
-j3
0
.i,
T”2 x3
Y3 0 0 0 $3
1
yy2
2
-jt,
. Y3
are
*
I
(48)
(49)
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
Fig. 3. Unstable
The matrix for the generalized Q’= -mg[O
1
0
1413
motion of arch.
forces due to gravity is 1
0
l].
In Fig. 3 an unstable motion starting from an almost symmetric
(50)
rest position is shown.
Acknowledgment Presently, the Institute of Mechanics of the Swiss Federal Institute of Technology and the Center for New Technologies of Oerlikon Industries are jointly developping a general purpose software tool based on the theory presented here. We appreciate the support of the Kommission zur Forderung der wissenschaftlichen Forschung (a research fund of the Swiss Federal Department of Economic Affairs) for this project.
References
PI
H. Brauchli and R.W. Weber, Canonical approach to multibody systems using redundant coordinates, in: G. Bianchi and W. Schiehlen, eds., IUTAM-IFToMM Symposium Udine 1985 (Springer, Berlin, 1986). PI H. Brauchli, Mass-orthogonal formulation of equations of motion for multibody systems, Z. Angew. Math. Phys. 42 (1991) 169-182. zur Beschreibung ebener H. Brauchli, G. Devaquet, Chr. Schaerer and J. Rohrer, Ein Baukastensystem [31 (Springer, Berlin, 1987). Vielkorpersysteme, in: J. Halin, ed., Proc. 4. Symposiums Simulationstechnik
1414 [4] H. Brauchli,
H. Brauchli, R. Weber, Dynamical equations in natural coordinates
On the description of constrained multibody systems in natural coordinates (in preparation). [5] J. Unda, J.M. Jimenez. A. Avello and J. Garcia de Jalon, Computer analysis of multibody systems with “natural coordinates”. in: G. Bianchi and W. Schiehlen, eds.. IUTAM-IFToMM Symposium Udine 198.5 (Springer, Berlin, 1986). [6] J.L. Synge. On the geometry of dynamics, Phil. Trans. A 226 (1927) 31-106. [7] St. Kaufmann, An explicit formula for the discussion of stability in autonomous Hamiltonian systems, Z. Angew. Math. Phys. 37 (1986) 735-753. [8] R. Clavel, Delta, a fast robot with parallel geometry, Proc. Internat. Symposium on Industrial Robots (1988) 91-100.