413
Mathematics and Computers in Simulation 29 (1987) 413-420 North-Holland
A METHOD FOR THE NUMERICAL INTEGRATION OF MECHANICAL SYSTEMS WITH UNILATERAL CONSTRAINTS: STUDY OF IMPACT IN MULTIBODY SYSTEMS * Vojin DRENOVAC Mathematical
Institute,
Belgrade,
Yugoslavia
A&tract: Studyof mechanical systanswith unilateral constraints is associated with forming two systemsof equations, a systemof differential equations 4 a systemof algebraical equations. Differential equations are used to describethe lrotion untilthe nmentof impact, i.e.untilactivation of unilateral constraints. Algebraic equations are used to describethe -ct. Duringnlrmerical integration, transition franone systemto anotheroccursat the pints of *ct. hren in simpleproblans, formingalgebraic equations represents a cunplex task. Thispaperpresentsa method,the so-called Reduction Methcd,whichprovidesfor the analysisof thesesystemswithoutformingthe algebraic equations. They are substituted by a new systemwhichis easilyderivedfranequations of na>tion. Comparedto ITH~XYIS basedon the classical impacttheory,usingReduction Methcd,velocities afterthe impact are easily canpltedrwardlessof the degreesof freedan. 1. Introduction Unilateral constraints have local effect on the motion of mechanical systems, restricted to the study of impact. These constraints are active in link elements such as joints, guides etc., and can be modelled with massless rods. As illustration, let us observe motion of a ball connected to a fixed point by a massless string. Effect of the string, i.e. the unilateral constraint is nonexistent, until the string is extended. Tight string corresponds to the impact and leads to the change of the ball velocity. The described system can be represented by the model displayed in Fig.1, where the string has been replaced by two massless rods.
Fig.1:Replacing the stringwith two mSSleSS
rods
* Thispaperwas writtenwhilethe authorwas at the Institute B for Mechanics of the Stuttgart university, (DEAD). suppXtedbytheGer?Mn Exchangeoffice
037%4754/87/$3.50
0 1987, IMACS/Elsevier
Science Publishers B.V. (North-Holland)
V. Drenovac / Impact in multibody systems
414
By applying this modelling approach to mechanical systems with unilateral constraints, equations of motion with a semidefinite mass matrix are obtained. Regardless of the numerical method used, numerical integration is interrupted in singular points. These singular points correspond to the points of impact.
2. Eguations of motion Equations of motion for a mechanical system with f degrees of freedom are M(ytt)j;+ k(y,j,t) = q(y,j,t)
(2.1)
where y is fxl vector comprising all generalized coordinates, and M(y,t) is a fxf mass matrix. The fxl vector k describes all Coriolis and centrifugal forces, and the fxl vector q describes all generalized forces, 161. The mass matrix is always symmetric and positive definite in real mechanical systems. In case the unilateral constraints are taken into account when forming the equations of motion, this matrix will be semidefinite. Equations (2.1) transformed into a form suitable for numerical integration are ; = Y(x,t)
(2.2)
where x is a 2fxl vector (x=[yTiTIT), and Y(x,t) is 2fxl vector function r
Y(x,t)
=
Y
I
M-'(q-k)
I.
3. Classification and the Reduction Method Numerical integration of system (2.2) cannot be performed in case of a semidefinite mass matrix. The deter?inant of the mass matrix M is zero at the singular points and therefore the inverse matrix M is not defined. From the mathematical point of view, equations (2.2) can be observed in the neighbourhood of singular points as a system of equations with small parameters at derivatives. The determinant of matrix M is in this case very small and can be treated as a small parameter. In the study of such systems, the problem is to find a solution when all small parameters converge to zero. This purely mathematical problem has been studied at length in literature. Let us observe a system of differential equations with small parameters at derivatives
dy -=
p =diag{ul,..,ur},
g(y,z,t),
rSf, s=f-r,
dt
(3.1)
dz v - = F(y,z,t) dt
for r
p = diag{$,..,ur,l&), S
where y=(yl,..,yf) and z=(yl 1..rY ‘f) are vectors and F=(Fl,..,F') and.g=(gl,..,gf) are vector functions. Without loss of generality. it can be assumed here that ul=u (i=l(l)r), where u is also a small parameter. The-solution‘ of the system (3.1) r(t,v),
dt,v)
(3.2)
is defined by the initial conditions Y(L,P)
= yo
,
z(to,lJ)
= zo
.
V. Drenovac
For u=O, system (3.1)
/ Impact
in multibody
415
systems
yields a degenerative system
dy - = g(ytz,t)t dt
(3.3)
z = @(y,t). z=a(y,t) is the solution of the sistem F(y,z,t)=O. The question is when and under what conditions will the solution (3.2) converge to the solution of the degenerative system, when II+0 , and under the same initial conditions. The answer is provided by Tichonow's theory [3],[4],[5]. Equations (2.2) have special structure compared to equations (3.1). First f equations are trivial equations of the form y=y. Remaining equations have small parameters at derivatives in the neighbourhood of the singular point. Right hand side of these equations has the form adjM(q-k). Respective to the number of small parameters, there are two cases: a) the number of small parameters is equal to the number of generalized coordinates, r=f, b) the number of small parameters is smaller than the number of generalized coordinates, r
< E
9
z
f
tiy,t)
equation F(y,z,t)=O is not true. 2. Transition from system (3.1) to system (3.3) corresponds to the associated system dz - = F(ytz,t), d-r
(T is a parameter).
(3.4)
Characteristic of system (3.4) is that the isolated solution z=@(y,t) is here a stable characteristic solution, with y and t as parameters. 3. Only stable characteristic solutions in the sence of Liapunov are of interest. When numerically integrating equations (2.2), there is a neighbourhood of the singular point where these equations do not hold numerically. Numerical values of particular coordinates exceed the range of the computer causing an interruption of integration. Equations (2.2), in case b, can be represented in the form *k
k
Y =z, *j ==j Y 9 -k
nz
(3.5)
= Fk(yk,yj,zk,zj,t),
z*j = fj(yk,yj,zk,zj,t). Set of variables y i.e. z is divided in two subsets: one including all variables having small parameters at derivatives; the other subset is the complement of the first.
416
V. Drenovac / Impact in multibody systems For v=O, we have Fk(yk,yj,zk,zj,t)=O.
(3.6)
Only the stable characteristic solutions of the above equation are of interest. Thus, from equations (3.4), it follows that .a"=~:is a constant. Equations (3.5) now become .k
k
Y
=Z
-j
=
Y
c
zj 9
’
zj = fj(yk,y',zZ,z',t). that is j;'= f'(yk,yj,+j,t) where y
k
(3.7)
are linear functions.
Method: By eliminating respective columns and rows in equation (2.1) we obtain the same equations (3.7). We shall demonstrate that condit.ion (3.6) can be introduced into the system of equations (2.1). System (2.1) can be written in this form
Reduction
a,B = l(l)f
a,S j;@= S,
(3.8)
where aclSare elements of the mass matrix. Coordinates of the q-k are included in So. Without loss of generality, we can assume that only one coordinate, i.e., yk is singular. From the condition that yk=zE, it follows that yk=O. Column in aclSthat corresponds to this coordinate can be eliminated, so that we obtain akj j;j= Sk, .. a_ij y’ = Si,
(3.9) (i,j=l(l),..,k-l,k+l,..,f).
(3.10)
Coefficients aij represent the remaining elements of matrix M. From equations (3.9) and (3.10) it follows
aji Si = Sk. akj -..
(3.11)
Using the relation .ki ji=
akj
__
9
akk we obtain
lakS
S B = 0,
1 *a@ = a@, P
u = laoSI.
(3.12)
From equation (3.5) in the form Y
'k
= zk
,j
= zj
9
,,;k = ;j we
=
,
lakB
(3.13) s
8'
P =
la,6
It
ajBSB.
also obtain equations (3.12). Coefficients
lakS
represent only elements of the k-th row in
417
V. Drenovac / Impact in multibody systems
the adj[aoH]. They do not contain small parameter u. Matrix elements aj@ also do not contain the small parameters because they are contained in the elements of matrix adj[aoS]. According to the condition (3.6), equation (3.12) follows. From equations (3.13), we obtain the equation .. yl
= ajB
(3.14)
S 13
which corresponds to equation (3.10). Conditions (3.6) are introduced into these equations,[2].
4. Numerical integration Numerical integration of equations (2.2) with the semidefinite mass matrix necessitates a few explanations: i) When does the interruption of integration occur ? ii) How large is the neighbourhood of the singular point ? iii) How are the new initial conditions for integration defined after interruption ? The answers are the following: i) Figure 2 illustrates the usual dependency of the mass matrix determinant on time. The condition det M=m(t)bd, (d>O, d-O), excludes the points of integration interruption (m(t)=O) and must be examined after each integration step.
Fig.2:Interruption mall
1.
2.
3.
5.
5.
7.
8.
9.
18. T
of integration using
parameter
T
ii) Computing the neighbourhood of the singular point At is related to the feature of conservative systems with respect to the volume integral in the phase space, Liouville's theorem, ]ll* I 2f = ... dy 1...dyfdpl...dpf = const. (4.1) I i Generalized impulses pi correspond to generalized coordinates yi (i=l(l)f). In case of one t;g;;e of freedom system, after integrating in the interval {[-y,o],[o,p]), it follows from . y.j,= const.
(4.2)
where impuls p is substituted with y. In conservative systems with two degrees of freedom, by developing equations of motion up to the terms of second order , we obtain again equation corresponding to equation (4.2), [7]. In systems with many degrees of freedom where the singularity is caused by one coordinate only, i.e. yl, equation (4.1) yields again equation (4.2). Surface elements dyidpi (i=2(l)f) have no effect. From (4.2), we obtain for the interval At
* (4.3) where y * and jr*are values at the interruption point of integration. Let us assume that the
V. Drenovac / Impact in multibody systems
418
singularity is at the center of interval At. iii) Singular coordinates in the interval At are linear functions. We shall assume that their velocities correspond to the mean velocity $, Figure 3, 2Y* f=--._=_ At
2y*jr* = 2y* .
(4.4)
y*
As before, we shall distinguish between two cases, a and b. In the first case, values of singular coordinates at the end of interval At are defined through the mean velocity 7. In case b, integration of the reduced system (3.10) produces values of nonsingular coordinates and velocities. Singular coordinates and velocities are,defined as in a.
t
Y
At
5.
Fig.3: Meminig
of the
mean velocity
2Y*
Example
Mass point with mass m is linked with massless rods AM and AB whose length is respectively a and r, to the slider B with mass ml, Figure 4.
Fig.4: Impact
pnduhm
Equations of motion read as follows
vcos2@ + sin'@
!
nsin@ sino
nsin$ sina
n2
(1-v)sin@ co&
;d . 11
b2+ nsin@
+ a 1
nsincrcosb ;9'
where the following symbols are used v = m,/m ,
q = a/r ,
u2= g/a
(g .s the gravity constant).
Mass matrix becomes semidefinite for ml=O. The singular position is defined with b=O, and the interval At according to (4.3) is
419
V. Drenovac / Impact in multibody systems
at=
@ II -
.
(5.2)
6
Using the Reduction Method we obtain the angle value a.and the jump value &; in case b, 1 a = - sino cosr$$'- w2sinu rl
(5.3)
Here, @ is a linear function, and & is a constant. By integrating this equation, we obtain
,I,
=
cl’ + -
2
sincesin@ 6. rl Variables a' and &" are velocity values before and after the singular point (point of impact). Second term in (5.3) is small compared to the first one and thus can be disregarded. The same results are obtained using the classical impact theory, [B]. The results are illustrated in Fig.5, for small values of parameter u=O.lE-7 and v=O (o=@=n/4). Other values are: n=9.81, d=O.OOl, r=0.5, g=9.81. m .a .6 .4
--sM3ll
parameter
.I 8. B
reduction
-. * -. 4 -. 5 -. B B.B I. 5
ix
.5
1.B
1. 5
2. B
.5
3.5
4. B
!-_
r’
snlallparameter
F---k
5. B
.5
1.5
2. kl
Fig.5: Canparative results of integmting
‘I”“““““‘,““‘.“.
2. 5
i
-1-=t
reduction
B.B
1. 5
T
l.B-
----
3.8
2. 5
3. 8
3. 5
1. B
L. 5
5. B
With a small parameter (V=O.lE-7)
and the Reduction M&ho3
(V=O)
6. Conclusion This paper is a contribution to the study of mechanical systems with unilateral constraints. In such systems, the mass matrix is semidefinite as the result of modelling real mechanical systems with massless rods. Continuous integration of these systems is not possible regardless of the method used. Interruption occurs at singular points which correspond to the points of impact. Transition through each of these points, according to the classical impact
420
V. Drenovac / Impact in multibody systems
theory, necessitates forming an algebraic system of equations. The proposed Reduction Method relies, instead of on the algebraic system, on the so-called reduced system which is easily obtained from equations of motion, by eliminating certain rows and columns. Integration of this system produces discontinuities (jumps) of the velocity function for each point of impact.
References [I] Arnold, I. V., Mathematical Methods of Classical Mechanics (Springer-Verlag, New-York/.., [Z] kziz?;ac, V A Method for the Numerical Integration of Equations of Motion of Singular Multibody S;&ems (in German), Ph.D. Thesis, University of Stuttgart, (Feb. 1985). [3] Tichonow, N. A., System of differential equations containing small parameters at derivatives, Mat. sb. Vo1.31(73), No.3(1952) 575-586. [4] Tichonow, N. A., On the dependance of the solution of differential equations on small parameter, Mat. sb. Vo1.22(64), No.2(1948) 193-204. [5] Tichonow. N. A., On Systems of Differential Equations Containing Parameters, Mat. sb. Vo1.27(69), No.l(l950) 147-156. [6] Schiehlen, W. and Kreuzer. E., Generation of Equations of Motion of Multibody Systems Supported by a Computer (in German), Ingenieur-Archiv 46(1977) 185-194. [7] Schiehlen. W., Nonlinear Oscilations in Multibody Systems, Proc. of the IXth ICNO, Kiev (1981). [8] Wittenburg. J., Dynamics of Systems of Rigid Bodies (B. G. Teubner, Stuttgart, 1977).
Address: Dr. Vojin Drenovac Generala Zdanova 1 11000 Belgrade Yugoslavia
I have a great debt to Miss Vesna Taw&
for her help, criticism and translation.