- Email: [email protected]
On the application of the modal integration method to flexible multibody systems
Pergamon
Computers & Srruc~ures Vol. 59, No. 3, pp. 553-559. 1996 Copyright Q 1996 Elsevier Science Ltd Printed in Great Britain. All rights rexrved 00457949/96 $15.00 + 0.00
00457949(95)00268-5
ON THE APPLICATION OF THE MODAL METHOD TO FLEXIBLE MULTIBODY
INTEGRATION SYSTEMS
Jung-Hwan Lee Research
and Development
Department, Kumho and Co. Inc., 555 Sochon-dong, Kwangju, 506-040, Republic of South Korea (Received
23 September
Kwangsan-ku,
1994)
to flexible multibody systems. Even dynamics, this technique has been rarely used in multibody mechanical systems containing a mix of rigid and flexible bodies. However, there are several advantages to the modal integration technique when it is applied to so-called stiff differential equations. Due to the orthogonality of the modal matrix, we can obtain the uncoupled equations of motion. This allows exact solutions for the linearized equations of motion in modal coordinates. Therefore, the time step taken is only restricted by the range over which the linearization holds. This drastically increases the size of the time steps which can be used. Since modal integration is performed on linearized differential equations, linearization schemes are discussed in the first part of this paper. The modal integration method described here is applied to two representative mechanisms with flexible bodies, and the results are included.
Abstract-In
this paper,
the modal
integration
method
is applied
though the concept of this integration method is popular in structural
INTRODUCTION
Dynamic modeling of multibody mechanical system articulated by kinematic joints often ends up with systems of differential-algebraic equations, abbreviated DAE [l, 21. That is, equations of motion
W(r)l{?J+ {V, i, t)} = if) and constraint
(1)
equations; {C(r, r)> = (0).
(2)
In the literature, there are in principle two different ways to deal with differential equations with algebraic constraints. One way is to directly solve the simultaneous differential algebraic equations [l]. There are two problems pertaining to these methods. First of all, flexible multibody mechanical systems usually turn out to be stiff systems, i.e. characterized by an extremely wide range of natural frequencies. Another problem related to this method is that solutions often do not satisfy the constraint equations after several time steps are taken. These errors introduce large forces during the time integration. To correct this trouble, several stabilization methods have been proposed [3]. Considering the problems mentioned above, direct numerical integration of DAE is an ongoing active area of research. Due to these problems, methods based on transformation of variables are still widely used [4]. These methods first transform the differential equations into a reduced set of ordinary differential equations by use of constraint equations. Then the reduced system of 553
ordinary differential equations (ODE) are integrated. The literature on numerical integration of ODE is vast [5, 61. However, it is known that no algorithm is good for all problems, and choosing an optimum algorithm is problem-dependent [7]. In general, numerical integration methods for finding the dynamic response of a system are either direct integration methods or modal methods. The direct integration methods are further characterized as explicit or implicit methods. While numerical integration is performed directly on the given dynamic equations in the direct integration methods, the dynamic equations are transformed into the modal coordinates in the modal method. The modal method is often effective due to linear independence and orthogonality of the eigenvectors. These properties change coefficient matrices of the dynamic equations into a simple, uncoupled form with modal coordinates. Even though the concept of this integration method is popular in structural dynamics, this technique has been rarely used in flexible multibody system dynamics. This is because eigensystem calculation costs are high. However, there are several advantages to the modal integration technique. The motive for use of the modal technique is that one can have exact closed form solutions for the uncoupled equations in terms of the modal coordinates. Due to this fact, we can avoid the problems caused by inherently stiff systems. This also means that the time step taken is only restricted by the range over which the linearization holds. Therefore, the number of iterations in the total simulation is reduced drastically. Another fact that makes the modal integration technique attractive is
Jung-Hwan
554
that coordinate reduction is possible in the time integration step. An accurate time history analysis can often be achieved using only a small range of low frequency information. If this is the situation, coordinate reduction can be done through modal transformation before each time step is taken. Eigenvalues of the previous time step can also be used for initial estimates of eigenvalues in subsequent steps. In the literature, it has been suggested that the inverse iteration method is suitable for tracking the history of an eigenvalue [8]. The modal technique presented here has originally been proposed by Cipra [9] and implemented into the IMP program [IO]. In order to use the modal integration method, one needs to linearize the nonlinear differential equations shown in eqn (I). By differentiating constraint equations, we can get constraint equations in velocity and acceleration levels, in addition to position level constraint. If we assume that the linearization has been performed, then we can have a linearized equation together with constraint equations in the following form:
[A]{si}
+[C1{6i}
+[R]{6r}
= {.f)
Lee LINEARIZATION OF NONLINEAR DIFFERENTIAL EQUATIONS
Since we need to linearize a set of nonlinear differential equations for modal integration, linearization schemes will be discussed briefly in this section. For comparison of the linearization schemes described below, an example of a simple pendulum attached to a moving cart is chosen, as shown in Fig. I. Nonlinear equations of motion for this system are 2m
[ ml cos 0
ml case ml2 +
.i: Ii e I -mii’sin
0 + kx
mgl sin
= (0).
e
(8)
Scheme I Probably the most general linearization scheme is the following form which results from the perturbation of a nonlinear function g(r, i, P) = 0 about an equilibrium position or about some nominal trajectory.
(3)
? W.1 = WI W,“, j
W 1=
WI isi,,, 1
iw = Wl{~~“,}+ iq,
(4)
(5) (6)
where {Y} comprises the quadratic velocity terms resulting from differentiation of the constraint equations in eqn (2). Putting eqns (4)-(6) into eqn (3) and then premultiplying by [H]’ gives
WI {K”d1 + [Cl{bid >+ WI{bY,“d i = {.I”1.
(7)
where
Equation (7) represents the linearized equations of motion for the total system with minimal coordinates. It is worthwhile emphasizing here that the linearized equations are only valid during a small time step around the current operating position at time T. These equations are to be integrated to find the next operating position at T + dt from the information of the current operating position at time T. Then eqn (7) is updated to represent the linearized equation around the new operating position.
g=gu+6g=g~+~6i+~6i+~6r=O.
(9)
Such linearized equations are usually used for stability analysis of the motion, and for the study of sensitivity analysis to small disturbances and parameter variations [l I]. In order to use this linearization scheme, we need to know the nominal trajectory beforehand. However, we usually do not know the nominal trajectory except for special cases such as, for example, railway vehicles traveling at a known speed along a given track shape. Even though we may assume nominal trajectories, this could result in erroneous results. Therefore, it is not adequate to use the linearization of eqn (9) for the general case of finding the next operating position at T + dt from information of the current operating position at time T. If we assume that the motion of the cart and pendulum in Fig. 1 is small enough, then we can find the motion of the system around an equilibrium position by integrating these linearized equations. By using eqn (9), the linearized equations of motion
m= 0’2 lb
s&/in
1=19.3 In
Fig. I. Simple pendulum on a moving cart.
Modal integration method applied to flexible multibody around an equilibrium as follows:
ml
cos eO
0
0
[ k
06
se
ml2
- 2mI sin e,&
+0[
+0
SR 1-iI I(se I
mgl cos eO- ml sin
6x 0 x{68 = 0' 111
^ v
04
E ;
0.2
E E
6i
-mlsinB,&-mlcos0,6~
555
08
(0, = 0) are obtained
m~ cos eO
2m
[
position
systems
0
=-02 E: G -0.4
e,z 1
-0.6 -0
6 0 20
0 00
(10)
0 40 time
0 60
0 60
)
(set
I 00
Fig. 3. Angular displacement of the pendulum of Fig. I.
Scheme II There is another scheme for which we do not need to find the equilibrium position or nominal trajectory. In this scheme, perturbation about a current operating point at time T is represented as t=T+&
on the right hand side are assumed to be constant during the time step. The constant values for 6r, 6f and 6i: are found at 6t = 0, i.e. Gr(St = 0) = 0, &(6t = 0) = 6i(T = t) and 6f(& = 0) = iV(T = t). By this approach, the linearized equations of motion for the system shown in Fig. 1 are found as
[ ml cos e,
i =6i +
where, since linearization takes place about the operating point r* which is assumed to be constant, the differential of r * is zero. Using perturbation, Taylor’s expansion of function g(r, f, f) = 0 can then be written as g(r * + 6r, f, f) = g(r *, 63,6f) +~g(r,&,SP)(*br. a4
(11)
The higher-order terms resulting from this scheme are moved to the right side of the equations. Then, the variables 6r, 6i and 6i in the higher order terms
=
ml2
k
0
0
mgI ~0~8,
-ml
c5.f I{ 1
ml cos eO
2m r =r*+Sr
sine,@
se
- kx,
mgl sin eO
.. 1
(12)
For purposes of comparison and validation, the nonlinear equation in eqn (10) and the linearized equations of eqn (12) have been solved by the fourthorder Runge-Kutta method. The coefficients of the nonlinear equations, as well as of the linearized equations, were evaluated at every time step. In Figs 1 and 2, displacement x and 0 are obtained for both cases, and both cases show the same results. It has also been observed that the modal integration method described in the next section showed identical results for this example problem. MODAL TIME INTEGRATION METHOD It is assumed that the equations of motion have been obtained at the current time (T) and linearized about that operation point, as shown in eqn (7). The resulting linearized equations are represented as
0.6 04 d E
0.2
E
0
wfl(~ } + Pqa I+ Kl{x I= IF 17
H E-0.2 z -0.4 -0.6 ’ 0.00
0.20
’ 0 40 time
’
’ 0 60 (set)
’ 060
’
1
(13)
where matrices [Ml, [Cl, [K] and {F} are assumed to be constant over the coming (short) time interval. The above equations can now be transformed into firstorder form as
1 00
Fig. 2. Displacement of the moving cart of Fig. I.
[;
“C’](;}+[-;
;I{;}={;}.
(14)
Jung-Hwan Lee
556
By introducing obtain
A=[;
new
E]
Therefore, equations are
the
variables
and
{v} = [iixl’,
B=[-y
solutions
of
i].
the
we
the
initial
{Y(O)) = Pl{v)
(16)
PI {rl) = i Y(O)) - [Bl-’
The matrix {q} are obtained term [Q][[B] as follows:
Pl’P]P]l(rll
then the homogeneous generalized eigenvalue
(18)
form of eqn (15) becomes problem
(26)
+ P-’
; . 11
(27)
by premultiplying
by the
- [@]‘[B][B]-’ From the orthogonality
; 11
. (28)
of the eigenvectors,
the
(]9)
The solutions of eqn (19) give 2n eigenvalues (1) and 2n corresponding eigenvectors {rp}. With the eigenvalues and eigenvectors, the homogeneous solution becomes
the matrix [B] diagonalized
where
{II> = [~~,l-‘Pl’[4{Y(O))
Also rearranging
(21)
eL,’
by
t h e modal matrix [@I. Then eqn (27) becomes
-[\KJ’Pl’
[\B,l-‘[q’([Bl{l.(o)}
PI = [cpl(P2.. ~4ohl
above
= PlVl~Y(O)~
where [%I represents
G-LA1+ Wl){cp) = (0).
the
(17)
where { yC } and { y,} represent the homogeneous and particular solutions, respectively. If we assume the homogeneous solutions are of the form {JG> = icp} e”
conditions,
(15)
differential
{Y} = {VC) + {Y,f?
If we apply equation yields
- {;I).
(30)
eqn (29) yields
[B]-’ = [aq’[\8 \ I-‘[@].
(31)
0 ehl [\e*‘,]
=
.. 0 r
(22)
* e;.*r i
@I>= [rll 92. . ’ fl2J.
(23)
The column matrix fv} represents unknown modal amplitudes at time t = T which will be determined from the initial conditions for each time step. Since {F} is assumed constant over the short time interval, we have the particular solutions of eqn (15) as
Therefore
the general
solution
From eqns (30) and (31) the general solution eqn (24) becomes
for
1.~1 = Pl~\eA’,lPB\l-‘Pl’
(32) We can somewhat simplify the above equation by adding zero: -{y(O)} + {y(O)}. Since [El-‘[B] = I,
becomes (25)
- [W’[~l~YW~
+ {Y(O)).
(33)
557
Modal integration method applied to flexible multibody systems
Then, from eqn (31)
I
I -----
(~1 = ~~l~~~^‘,l~~B\l-‘~~l’PI{Y(O)I( (8
y
01
dt 0.0001 I
-
dt = 0.00
-
dt z 0.002
I
S z f a
- [~l[\B\l-‘[~l’[Bl{y(O)} + {Y(W~ (34)
0
5 .G -0. I
which reduces to {y} =[~][\e*‘-I,][\B,]-‘[~I
0.02
0.00
time
x([BllY(o)}-{~})+{Y(o)J.
(35)
If we explicitly express matrices [B] and {y(O)}, then
(rec.)
0.04
0 06
Fig. 5. Transverse displacement u2 of the coupler. the case. Denoting the modal vectors and eigenvalues corresponding to the reduced set of low frequencies by subscript 1, we have the closed form solutions as
{Yl = WWe:1- I,lN3,1-‘[@,I - {;} Since the initial conditions (35) can be simplified as
(36)
= {i:}.
x,, = x(8t = 0) = 0, eqn
X where [@,I=
{y} = [@][\e^’ - I,]@,]-‘[@I’ {
cp,l, W?,l
(37)
Equation (37) represents the exact solutions of the dynamic equations in eqn (13). By differentiating the above equation with respect to time, we can obtain the accelerations of x in the following equations:
= [Q] [A,] [\e*‘,][\B\]-‘[a]’
-yi” i
(38) 1
As mentioned in the second section, the modal coordinate reductions are possible for the time integration step. For design purposes, one may wish to see the wide range of the frequencies. However, only the lower range of frequencies are usually required for time history analysis with acceptable accuracy. One of the advantages in the final form of the solution in eqn (37) is that this equation exactly fits
J$+++ XC, a/ ’ wu, in its initial
and m is the number of the lower frequency modes used; m < n. We can recognize that [@,I, [\ey,], and [\&I are simply submatrices of [a], [\eAt,], and [G,], respectively. EXAMPLE
PROBLEMS
To demonstrate the validity of the simulation procedure described here, two representative examples of mechanisms-a slider-crank mechanism and a four-bar linkage-are selected. The slidercrank mechanism considered has been studied by several researchers [12] and is therefore chosen as a benchmark for comparison purposes. For the modeling of the mechanisms, the procedure shown in eqns (I)-(7) has been used. That is, each link has been separated from the mechanism and then modeled. Next, the dynamic equations for each link are linearized. At the same time, constraint equations are formulated. Finally, linearized system dynamic equations with minimum coordinates are obtained as shown in eqn (7). Slider -crank mechanism with massless slider
“I
Fig. 4. Slider-crank mechanism
.
1 + {Y(O)).
u.
[cp, v2..
-ri”
“.I configuration.
Figure 4 shows the slider-crank mechanism which consists of a rigid crank, a flexible connecting rod and a massless slider. The crank and the connecting rod are both assumed to have uniform circular cross-
Jung-Hwan
558
Lee
---dt:OOOl 01
--o-
dt T 0 002
k z 0 7
9 K 5
t_
0
i
-0
I
4% Fig. 8. Four-bar
-0 2 0 00
0 02
)
time (set Fig. 6. Displacement
0 04
u2 of the coupler angle is 90”.
0 06
when initial
crank
sectional area. The flexible connecting rod is modeled by two plane frame elements [13] which can deform both axially and in bending. Therefore, each node initially has three degrees of freedom as shown in Fig. 4. In Fig. 4, the body reference coordinate system, as well as the inertial coordinate system (X0 Y,,), are shown. The body reference coordinate systems are attached to the middle of each link in its initial configuration. These are attached to the bodies in such a way that x axes are always parallel to the line between the two pin joint axes during the motion. It is also assumed that the nodal displacements of the connecting rod are initially zero. Mechanical properties of the chosen mechanism are the same as those in the literature [12]. For the given mechanism with the initial positions shown in Fig. 4, the following cases are tested: dt =O.OOOl; dt = 0.001; and dt = 0.002. Figure 5 shows the transverse displacement of the node located at the middle of the connecting rod. The displacement is measured with respect to the reference coordinate system. The results with dt = 0.0001 show very good numerical agreement with the results in the literature. As we increase the step size, a noticeable time delay appears in the time history. This time delay is due to the initial configuration of the
02m
I __
full
tiechanism
in its initial position
a so-called “dead-center position”. mechanism, Figure 6 shows the results for the same node when the initial position of the crank is 90”. In this case, increase of the step size does not display time delay. Figure 7 is the result when frequency reduction at the time integration stage is performed. Eight modes of high frequencies are deleted from the first-order differential equations, as shown in eqn (39). This is equivalent to deleting four modes from the secondorder differential equations. As we can see in Fig. 7, this does not affect the displacement r2. Four-bar mechanism In this example, the connecting rod and output link (rocker) shown in Fig. 8 are considered to be flexible. Again, the crank is considered rigid. The data for the mechanism are: length of crank is 4in; length of coupler is 12 in; length of rocker is 6 in; and the distance between ground joints is 11 in. Figure 8 also shows the body reference coordinate systems attached to the middle of each body at its initial configuration. The flexible connecting rod and rocker are modeled with two plane frame elements the same as in the previous example. Therefore, each elastic link has six elastic degrees of freedom. All the links are assumed to be uniform circular rods. It is also assumed that the initial displacements of the nodes are zero. The properties of the four bar mechanism are: Young’s modulus of the rods is 30 x lOh lb in ‘; diameter of the rods is 0.25 in; mass density of the rods is 0.0007331 lb.~~~inP; and crank speed is 125 rad SC’ (constant).
I
modes
041i
.
I ---m-dt.00001 1 w(--dt:OOOl -02-
0 00
001
Fig. 7. Displacement
0 02
0 03 time (set
0 04
)
v2 of the connecting modes deleted.
0 05
0 06
-0
4
-
0 00 rod with eight Fig. 9. Transverse
dt : 0 002
001
0 02
time
displacement
1 0 03 (set
1
004
at the center
0 05
of coupler.
Modal integration method applied to flexible multibody
systems
559
a reasonably accurate time history with only 25 time steps for 360” rotation of crank has been obtained. It is also observed that we can obtain good approximate solutions with only a few low frequency modes. These results show that modal integration method in multibody dynamic simulations can be a useful tool.
REFERENCES
+ -04
0 00
16 modes
I 001
’ 0 02 time
Fig.
10. Displacement
deleted
’
(set
0 03 )
’
0 04
ur of the connecting modes deleted.
’
0 05
rod with
16
Transverse displacement of the node located at the middle of the connecting rod are shown in Fig. 9. These are tested with two finite elements on the coupler and output link. The following cases are tested: dt = 0.0001; dt = 0.001; and dt = 0.002. The results show good numerical agreement. Figure 10 shows the results when frequency reduction is performed at the time integration stage. It is observed that we can get a good approximate solution on displacement at the center of the connecting rod with 16 modes deleted in eqn (39). CONCLUSION
In this paper, the modal integration method has been presented for the dynamic analysis of multibody mechanical systems. The motive for use of the modal method was that we can avoid the problems caused by inherently stiff systems. Even though eigensystem calculations can be slower, there are several advantages as far as efficiency is concerned. The time step is only restricted by the range over which the linearization holds. In addition, an accurate time history analysis can usually be achieved using only a small range of low frequency information. The procedure proposed here was applied to the two representative mechanisms. The results obtained are comparable to the results in the literature. It has been observed that
1. E. J. Haug (ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, NATO ASI Series F, Vol. 9. Springer, Berlin (1984). 2. A. Shabana, Dynamics of Multibody Systems. Wiley, New York (1989). 3. G.-P. Ostermeyer, On Baumgarte stabilization for differential algebraic equations. In: Real-Time fntegration Methods for Mechanical System Simulation (Edited by E. J. Haug and R. C. Deyo), NATO ASI Series F, Vol. 69, pp. 193-207. Springer, Berlin (1991). 4. W. Blajer. D. Bestle and W. Schiehlen, An orthogonal complement matrix formulation for constrained multibody systems. J. Mech. Des. Trans. ASME 116,423428 (1994). 5. C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, NJ (1971). 6. E. Isaacson and H. B. Keller, Analysis of Numerical Methods. Wiley, New York (1966). 7. K. J. Bathe and E. L. Wilson, Numerical Solution in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, NJ (1971). 8. G. Peter and J. H. Wilkinson, Ax = IBx and generalized eigenproblem. SIAM J. numer. Anal. 7, 479492 (1970). 9. R. J. Cipra and J. J. Uicker Jr, On the dynamic simulation of large nonlinear mechanical system: part I -an overview of the simulation technique. Substructuring and frequency domain consideration. J. Mech. Des. Trans. ASME 103, 849-856 (1981). 10. P. N. Sheth and J. J. Uicker Jr, IMP (integrated mechanism program), a computer-aided design analysis system for mechanism and linkage. J. Engng Ind. Trans. ASME 94, 454464 (1972). Principles of Dynamics, 2nd edn. Il. D. T. Greenwood, Prentice-Hall, Englewood Cliffs, NJ (1988). and R. A. Wehage, Variable degree of 12. A. Shabana freedom component mode analysis of inertia variant flexible mechanical systems. J. Mech. Transmis. Automat. Des. 105, 371-378 (1983). 13. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd edn. Wiley, New York (1989).
Computers & Srruc~ures Vol. 59, No. 3, pp. 553-559. 1996 Copyright Q 1996 Elsevier Science Ltd Printed in Great Britain. All rights rexrved 00457949/96 $15.00 + 0.00
00457949(95)00268-5
ON THE APPLICATION OF THE MODAL METHOD TO FLEXIBLE MULTIBODY
INTEGRATION SYSTEMS
Jung-Hwan Lee Research
and Development
Department, Kumho and Co. Inc., 555 Sochon-dong, Kwangju, 506-040, Republic of South Korea (Received
23 September
Kwangsan-ku,
1994)
to flexible multibody systems. Even dynamics, this technique has been rarely used in multibody mechanical systems containing a mix of rigid and flexible bodies. However, there are several advantages to the modal integration technique when it is applied to so-called stiff differential equations. Due to the orthogonality of the modal matrix, we can obtain the uncoupled equations of motion. This allows exact solutions for the linearized equations of motion in modal coordinates. Therefore, the time step taken is only restricted by the range over which the linearization holds. This drastically increases the size of the time steps which can be used. Since modal integration is performed on linearized differential equations, linearization schemes are discussed in the first part of this paper. The modal integration method described here is applied to two representative mechanisms with flexible bodies, and the results are included.
Abstract-In
this paper,
the modal
integration
method
is applied
though the concept of this integration method is popular in structural
INTRODUCTION
Dynamic modeling of multibody mechanical system articulated by kinematic joints often ends up with systems of differential-algebraic equations, abbreviated DAE [l, 21. That is, equations of motion
W(r)l{?J+ {V, i, t)} = if) and constraint
(1)
equations; {C(r, r)> = (0).
(2)
In the literature, there are in principle two different ways to deal with differential equations with algebraic constraints. One way is to directly solve the simultaneous differential algebraic equations [l]. There are two problems pertaining to these methods. First of all, flexible multibody mechanical systems usually turn out to be stiff systems, i.e. characterized by an extremely wide range of natural frequencies. Another problem related to this method is that solutions often do not satisfy the constraint equations after several time steps are taken. These errors introduce large forces during the time integration. To correct this trouble, several stabilization methods have been proposed [3]. Considering the problems mentioned above, direct numerical integration of DAE is an ongoing active area of research. Due to these problems, methods based on transformation of variables are still widely used [4]. These methods first transform the differential equations into a reduced set of ordinary differential equations by use of constraint equations. Then the reduced system of 553
ordinary differential equations (ODE) are integrated. The literature on numerical integration of ODE is vast [5, 61. However, it is known that no algorithm is good for all problems, and choosing an optimum algorithm is problem-dependent [7]. In general, numerical integration methods for finding the dynamic response of a system are either direct integration methods or modal methods. The direct integration methods are further characterized as explicit or implicit methods. While numerical integration is performed directly on the given dynamic equations in the direct integration methods, the dynamic equations are transformed into the modal coordinates in the modal method. The modal method is often effective due to linear independence and orthogonality of the eigenvectors. These properties change coefficient matrices of the dynamic equations into a simple, uncoupled form with modal coordinates. Even though the concept of this integration method is popular in structural dynamics, this technique has been rarely used in flexible multibody system dynamics. This is because eigensystem calculation costs are high. However, there are several advantages to the modal integration technique. The motive for use of the modal technique is that one can have exact closed form solutions for the uncoupled equations in terms of the modal coordinates. Due to this fact, we can avoid the problems caused by inherently stiff systems. This also means that the time step taken is only restricted by the range over which the linearization holds. Therefore, the number of iterations in the total simulation is reduced drastically. Another fact that makes the modal integration technique attractive is
Jung-Hwan
554
that coordinate reduction is possible in the time integration step. An accurate time history analysis can often be achieved using only a small range of low frequency information. If this is the situation, coordinate reduction can be done through modal transformation before each time step is taken. Eigenvalues of the previous time step can also be used for initial estimates of eigenvalues in subsequent steps. In the literature, it has been suggested that the inverse iteration method is suitable for tracking the history of an eigenvalue [8]. The modal technique presented here has originally been proposed by Cipra [9] and implemented into the IMP program [IO]. In order to use the modal integration method, one needs to linearize the nonlinear differential equations shown in eqn (I). By differentiating constraint equations, we can get constraint equations in velocity and acceleration levels, in addition to position level constraint. If we assume that the linearization has been performed, then we can have a linearized equation together with constraint equations in the following form:
[A]{si}
+[C1{6i}
+[R]{6r}
= {.f)
Lee LINEARIZATION OF NONLINEAR DIFFERENTIAL EQUATIONS
Since we need to linearize a set of nonlinear differential equations for modal integration, linearization schemes will be discussed briefly in this section. For comparison of the linearization schemes described below, an example of a simple pendulum attached to a moving cart is chosen, as shown in Fig. I. Nonlinear equations of motion for this system are 2m
[ ml cos 0
ml case ml2 +
.i: Ii e I -mii’sin
0 + kx
mgl sin
= (0).
e
(8)
Scheme I Probably the most general linearization scheme is the following form which results from the perturbation of a nonlinear function g(r, i, P) = 0 about an equilibrium position or about some nominal trajectory.
(3)
? W.1 = WI W,“, j
W 1=
WI isi,,, 1
iw = Wl{~~“,}+ iq,
(4)
(5) (6)
where {Y} comprises the quadratic velocity terms resulting from differentiation of the constraint equations in eqn (2). Putting eqns (4)-(6) into eqn (3) and then premultiplying by [H]’ gives
WI {K”d1 + [Cl{bid >+ WI{bY,“d i = {.I”1.
(7)
where
Equation (7) represents the linearized equations of motion for the total system with minimal coordinates. It is worthwhile emphasizing here that the linearized equations are only valid during a small time step around the current operating position at time T. These equations are to be integrated to find the next operating position at T + dt from the information of the current operating position at time T. Then eqn (7) is updated to represent the linearized equation around the new operating position.
g=gu+6g=g~+~6i+~6i+~6r=O.
(9)
Such linearized equations are usually used for stability analysis of the motion, and for the study of sensitivity analysis to small disturbances and parameter variations [l I]. In order to use this linearization scheme, we need to know the nominal trajectory beforehand. However, we usually do not know the nominal trajectory except for special cases such as, for example, railway vehicles traveling at a known speed along a given track shape. Even though we may assume nominal trajectories, this could result in erroneous results. Therefore, it is not adequate to use the linearization of eqn (9) for the general case of finding the next operating position at T + dt from information of the current operating position at time T. If we assume that the motion of the cart and pendulum in Fig. 1 is small enough, then we can find the motion of the system around an equilibrium position by integrating these linearized equations. By using eqn (9), the linearized equations of motion
m= 0’2 lb
s&/in
1=19.3 In
Fig. I. Simple pendulum on a moving cart.
Modal integration method applied to flexible multibody around an equilibrium as follows:
ml
cos eO
0
0
[ k
06
se
ml2
- 2mI sin e,&
+0[
+0
SR 1-iI I(se I
mgl cos eO- ml sin
6x 0 x{68 = 0' 111
^ v
04
E ;
0.2
E E
6i
-mlsinB,&-mlcos0,6~
555
08
(0, = 0) are obtained
m~ cos eO
2m
[
position
systems
0
=-02 E: G -0.4
e,z 1
-0.6 -0
6 0 20
0 00
(10)
0 40 time
0 60
0 60
)
(set
I 00
Fig. 3. Angular displacement of the pendulum of Fig. I.
Scheme II There is another scheme for which we do not need to find the equilibrium position or nominal trajectory. In this scheme, perturbation about a current operating point at time T is represented as t=T+&
on the right hand side are assumed to be constant during the time step. The constant values for 6r, 6f and 6i: are found at 6t = 0, i.e. Gr(St = 0) = 0, &(6t = 0) = 6i(T = t) and 6f(& = 0) = iV(T = t). By this approach, the linearized equations of motion for the system shown in Fig. 1 are found as
[ ml cos e,
i =6i +
where, since linearization takes place about the operating point r* which is assumed to be constant, the differential of r * is zero. Using perturbation, Taylor’s expansion of function g(r, f, f) = 0 can then be written as g(r * + 6r, f, f) = g(r *, 63,6f) +~g(r,&,SP)(*br. a4
(11)
The higher-order terms resulting from this scheme are moved to the right side of the equations. Then, the variables 6r, 6i and 6i in the higher order terms
=
ml2
k
0
0
mgI ~0~8,
-ml
c5.f I{ 1
ml cos eO
2m r =r*+Sr
sine,@
se
- kx,
mgl sin eO
.. 1
(12)
For purposes of comparison and validation, the nonlinear equation in eqn (10) and the linearized equations of eqn (12) have been solved by the fourthorder Runge-Kutta method. The coefficients of the nonlinear equations, as well as of the linearized equations, were evaluated at every time step. In Figs 1 and 2, displacement x and 0 are obtained for both cases, and both cases show the same results. It has also been observed that the modal integration method described in the next section showed identical results for this example problem. MODAL TIME INTEGRATION METHOD It is assumed that the equations of motion have been obtained at the current time (T) and linearized about that operation point, as shown in eqn (7). The resulting linearized equations are represented as
0.6 04 d E
0.2
E
0
wfl(~ } + Pqa I+ Kl{x I= IF 17
H E-0.2 z -0.4 -0.6 ’ 0.00
0.20
’ 0 40 time
’
’ 0 60 (set)
’ 060
’
1
(13)
where matrices [Ml, [Cl, [K] and {F} are assumed to be constant over the coming (short) time interval. The above equations can now be transformed into firstorder form as
1 00
Fig. 2. Displacement of the moving cart of Fig. I.
[;
“C’](;}+[-;
;I{;}={;}.
(14)
Jung-Hwan Lee
556
By introducing obtain
A=[;
new
E]
Therefore, equations are
the
variables
and
{v} = [iixl’,
B=[-y
solutions
of
i].
the
we
the
initial
{Y(O)) = Pl{v)
(16)
PI {rl) = i Y(O)) - [Bl-’
The matrix {q} are obtained term [Q][[B] as follows:
Pl’P]P]l(rll
then the homogeneous generalized eigenvalue
(18)
form of eqn (15) becomes problem
(26)
+ P-’
; . 11
(27)
by premultiplying
by the
- [@]‘[B][B]-’ From the orthogonality
; 11
. (28)
of the eigenvectors,
the
(]9)
The solutions of eqn (19) give 2n eigenvalues (1) and 2n corresponding eigenvectors {rp}. With the eigenvalues and eigenvectors, the homogeneous solution becomes
the matrix [B] diagonalized
where
{II> = [~~,l-‘Pl’[4{Y(O))
Also rearranging
(21)
eL,’
by
t h e modal matrix [@I. Then eqn (27) becomes
-[\KJ’Pl’
[\B,l-‘[q’([Bl{l.(o)}
PI = [cpl(P2.. ~4ohl
above
= PlVl~Y(O)~
where [%I represents
G-LA1+ Wl){cp) = (0).
the
(17)
where { yC } and { y,} represent the homogeneous and particular solutions, respectively. If we assume the homogeneous solutions are of the form {JG> = icp} e”
conditions,
(15)
differential
{Y} = {VC) + {Y,f?
If we apply equation yields
- {;I).
(30)
eqn (29) yields
[B]-’ = [aq’[\8 \ I-‘[@].
(31)
0 ehl [\e*‘,]
=
.. 0 r
(22)
* e;.*r i
@I>= [rll 92. . ’ fl2J.
(23)
The column matrix fv} represents unknown modal amplitudes at time t = T which will be determined from the initial conditions for each time step. Since {F} is assumed constant over the short time interval, we have the particular solutions of eqn (15) as
Therefore
the general
solution
From eqns (30) and (31) the general solution eqn (24) becomes
for
1.~1 = Pl~\eA’,lPB\l-‘Pl’
(32) We can somewhat simplify the above equation by adding zero: -{y(O)} + {y(O)}. Since [El-‘[B] = I,
becomes (25)
- [W’[~l~YW~
+ {Y(O)).
(33)
557
Modal integration method applied to flexible multibody systems
Then, from eqn (31)
I
I -----
(~1 = ~~l~~~^‘,l~~B\l-‘~~l’PI{Y(O)I( (8
y
01
dt 0.0001 I
-
dt = 0.00
-
dt z 0.002
I
S z f a
- [~l[\B\l-‘[~l’[Bl{y(O)} + {Y(W~ (34)
0
5 .G -0. I
which reduces to {y} =[~][\e*‘-I,][\B,]-‘[~I
0.02
0.00
time
x([BllY(o)}-{~})+{Y(o)J.
(35)
If we explicitly express matrices [B] and {y(O)}, then
(rec.)
0.04
0 06
Fig. 5. Transverse displacement u2 of the coupler. the case. Denoting the modal vectors and eigenvalues corresponding to the reduced set of low frequencies by subscript 1, we have the closed form solutions as
{Yl = WWe:1- I,lN3,1-‘[@,I - {;} Since the initial conditions (35) can be simplified as
(36)
= {i:}.
x,, = x(8t = 0) = 0, eqn
X where [@,I=
{y} = [@][\e^’ - I,]@,]-‘[@I’ {
cp,l, W?,l
(37)
Equation (37) represents the exact solutions of the dynamic equations in eqn (13). By differentiating the above equation with respect to time, we can obtain the accelerations of x in the following equations:
= [Q] [A,] [\e*‘,][\B\]-‘[a]’
-yi” i
(38) 1
As mentioned in the second section, the modal coordinate reductions are possible for the time integration step. For design purposes, one may wish to see the wide range of the frequencies. However, only the lower range of frequencies are usually required for time history analysis with acceptable accuracy. One of the advantages in the final form of the solution in eqn (37) is that this equation exactly fits
J$+++ XC, a/ ’ wu, in its initial
and m is the number of the lower frequency modes used; m < n. We can recognize that [@,I, [\ey,], and [\&I are simply submatrices of [a], [\eAt,], and [G,], respectively. EXAMPLE
PROBLEMS
To demonstrate the validity of the simulation procedure described here, two representative examples of mechanisms-a slider-crank mechanism and a four-bar linkage-are selected. The slidercrank mechanism considered has been studied by several researchers [12] and is therefore chosen as a benchmark for comparison purposes. For the modeling of the mechanisms, the procedure shown in eqns (I)-(7) has been used. That is, each link has been separated from the mechanism and then modeled. Next, the dynamic equations for each link are linearized. At the same time, constraint equations are formulated. Finally, linearized system dynamic equations with minimum coordinates are obtained as shown in eqn (7). Slider -crank mechanism with massless slider
“I
Fig. 4. Slider-crank mechanism
.
1 + {Y(O)).
u.
[cp, v2..
-ri”
“.I configuration.
Figure 4 shows the slider-crank mechanism which consists of a rigid crank, a flexible connecting rod and a massless slider. The crank and the connecting rod are both assumed to have uniform circular cross-
Jung-Hwan
558
Lee
---dt:OOOl 01
--o-
dt T 0 002
k z 0 7
9 K 5
t_
0
i
-0
I
4% Fig. 8. Four-bar
-0 2 0 00
0 02
)
time (set Fig. 6. Displacement
0 04
u2 of the coupler angle is 90”.
0 06
when initial
crank
sectional area. The flexible connecting rod is modeled by two plane frame elements [13] which can deform both axially and in bending. Therefore, each node initially has three degrees of freedom as shown in Fig. 4. In Fig. 4, the body reference coordinate system, as well as the inertial coordinate system (X0 Y,,), are shown. The body reference coordinate systems are attached to the middle of each link in its initial configuration. These are attached to the bodies in such a way that x axes are always parallel to the line between the two pin joint axes during the motion. It is also assumed that the nodal displacements of the connecting rod are initially zero. Mechanical properties of the chosen mechanism are the same as those in the literature [12]. For the given mechanism with the initial positions shown in Fig. 4, the following cases are tested: dt =O.OOOl; dt = 0.001; and dt = 0.002. Figure 5 shows the transverse displacement of the node located at the middle of the connecting rod. The displacement is measured with respect to the reference coordinate system. The results with dt = 0.0001 show very good numerical agreement with the results in the literature. As we increase the step size, a noticeable time delay appears in the time history. This time delay is due to the initial configuration of the
02m
I __
full
tiechanism
in its initial position
a so-called “dead-center position”. mechanism, Figure 6 shows the results for the same node when the initial position of the crank is 90”. In this case, increase of the step size does not display time delay. Figure 7 is the result when frequency reduction at the time integration stage is performed. Eight modes of high frequencies are deleted from the first-order differential equations, as shown in eqn (39). This is equivalent to deleting four modes from the secondorder differential equations. As we can see in Fig. 7, this does not affect the displacement r2. Four-bar mechanism In this example, the connecting rod and output link (rocker) shown in Fig. 8 are considered to be flexible. Again, the crank is considered rigid. The data for the mechanism are: length of crank is 4in; length of coupler is 12 in; length of rocker is 6 in; and the distance between ground joints is 11 in. Figure 8 also shows the body reference coordinate systems attached to the middle of each body at its initial configuration. The flexible connecting rod and rocker are modeled with two plane frame elements the same as in the previous example. Therefore, each elastic link has six elastic degrees of freedom. All the links are assumed to be uniform circular rods. It is also assumed that the initial displacements of the nodes are zero. The properties of the four bar mechanism are: Young’s modulus of the rods is 30 x lOh lb in ‘; diameter of the rods is 0.25 in; mass density of the rods is 0.0007331 lb.~~~inP; and crank speed is 125 rad SC’ (constant).
I
modes
041i
.
I ---m-dt.00001 1 w(--dt:OOOl -02-
0 00
001
Fig. 7. Displacement
0 02
0 03 time (set
0 04
)
v2 of the connecting modes deleted.
0 05
0 06
-0
4
-
0 00 rod with eight Fig. 9. Transverse
dt : 0 002
001
0 02
time
displacement
1 0 03 (set
1
004
at the center
0 05
of coupler.
Modal integration method applied to flexible multibody
systems
559
a reasonably accurate time history with only 25 time steps for 360” rotation of crank has been obtained. It is also observed that we can obtain good approximate solutions with only a few low frequency modes. These results show that modal integration method in multibody dynamic simulations can be a useful tool.
REFERENCES
+ -04
0 00
16 modes
I 001
’ 0 02 time
Fig.
10. Displacement
deleted
’
(set
0 03 )
’
0 04
ur of the connecting modes deleted.
’
0 05
rod with
16
Transverse displacement of the node located at the middle of the connecting rod are shown in Fig. 9. These are tested with two finite elements on the coupler and output link. The following cases are tested: dt = 0.0001; dt = 0.001; and dt = 0.002. The results show good numerical agreement. Figure 10 shows the results when frequency reduction is performed at the time integration stage. It is observed that we can get a good approximate solution on displacement at the center of the connecting rod with 16 modes deleted in eqn (39). CONCLUSION
In this paper, the modal integration method has been presented for the dynamic analysis of multibody mechanical systems. The motive for use of the modal method was that we can avoid the problems caused by inherently stiff systems. Even though eigensystem calculations can be slower, there are several advantages as far as efficiency is concerned. The time step is only restricted by the range over which the linearization holds. In addition, an accurate time history analysis can usually be achieved using only a small range of low frequency information. The procedure proposed here was applied to the two representative mechanisms. The results obtained are comparable to the results in the literature. It has been observed that
1. E. J. Haug (ed.), Computer Aided Analysis and Optimization of Mechanical System Dynamics, NATO ASI Series F, Vol. 9. Springer, Berlin (1984). 2. A. Shabana, Dynamics of Multibody Systems. Wiley, New York (1989). 3. G.-P. Ostermeyer, On Baumgarte stabilization for differential algebraic equations. In: Real-Time fntegration Methods for Mechanical System Simulation (Edited by E. J. Haug and R. C. Deyo), NATO ASI Series F, Vol. 69, pp. 193-207. Springer, Berlin (1991). 4. W. Blajer. D. Bestle and W. Schiehlen, An orthogonal complement matrix formulation for constrained multibody systems. J. Mech. Des. Trans. ASME 116,423428 (1994). 5. C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, Englewood Cliffs, NJ (1971). 6. E. Isaacson and H. B. Keller, Analysis of Numerical Methods. Wiley, New York (1966). 7. K. J. Bathe and E. L. Wilson, Numerical Solution in Finite Element Analysis. Prentice-Hall, Englewood Cliffs, NJ (1971). 8. G. Peter and J. H. Wilkinson, Ax = IBx and generalized eigenproblem. SIAM J. numer. Anal. 7, 479492 (1970). 9. R. J. Cipra and J. J. Uicker Jr, On the dynamic simulation of large nonlinear mechanical system: part I -an overview of the simulation technique. Substructuring and frequency domain consideration. J. Mech. Des. Trans. ASME 103, 849-856 (1981). 10. P. N. Sheth and J. J. Uicker Jr, IMP (integrated mechanism program), a computer-aided design analysis system for mechanism and linkage. J. Engng Ind. Trans. ASME 94, 454464 (1972). Principles of Dynamics, 2nd edn. Il. D. T. Greenwood, Prentice-Hall, Englewood Cliffs, NJ (1988). and R. A. Wehage, Variable degree of 12. A. Shabana freedom component mode analysis of inertia variant flexible mechanical systems. J. Mech. Transmis. Automat. Des. 105, 371-378 (1983). 13. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd edn. Wiley, New York (1989).