- Email: [email protected]
Dynamics of planar flexible mechanisms by finite element method with truss-type elements
Computers & Structures Vol. 39. No. l/2. pp. 13%140, Printed in Great Britain.
1991
0
MMS-7949/91 $3.00 + 0.00 1991 Pergmon Press plc
DYNAMICS OF PLANAR FLEXIBLE MECHANISMS BY FINITE ELEMENT METHOD WITH TRUSS-TYPE ELEMENTS M. HAM Institute of Machine Design Fundamentals, Warsaw Technical University, Narbutta 84, 02-524 Warszawa, Poland (Received 26 February 1990)
Abstract-The kinematics and dynamics of planar mechanisms with revolute joints and flexible links are studied by using the finite element method. The mechanism finite element is based on a planar truss element. The equations of motion of a mechanism are formulated by applying Lagrange’s equations to the kinetic and potential energy of the system. The analysis includes the higher order elastic acceleration terms. The presented method of solving equations of motion is based on a finite difference procedure and avoids the necessity of elimination of the singularity of the system stiffness matrix. A four-bar linkage is studied to illustrate the computational efficiency of the presented method. The results for flexible links are compared with those obtained for the rigid body.
1. INTRODUCTION
elements. In order to properly model a revolute joint of a mechanism? a distinct rotary freedom should be assigned to each of the elements connected by the node that replaces the revolute joint. This has been done by assuming not only longitudinal but also transverse displacement of the element nodes (i.e. it is possible for the element to move and rotate in any direction, but the stresses occur only due to the elongation of the element). Lagrange’s equations are used to obtain the equations of motion of a mechanism. The coefficient matrices are functions of the linkage geometry and vary as the input (crank) angle is varied. In order to solve the equations of motion the finite difference method with respect to time is used. In this method the problem of singularity of the system stiffness matrix does not exist. Contrary to most of the other publications, the variable input torque is considered here. As an example a kinematic and dynamic analysis of a four-bar linkage is presented. The element stiffness matrix is based on the planar truss element. It is asumed that maximum stresses occurring in links do not exceed critical values. Such an approach for the modelling of four-bar linkages was suggested in [6]. Other methodologists [4, 7,8] used a series of beam elements to model a four-bar linkage. for truss
In the past few years the number of publications which deal with elastic behaviour of mechanisms has increased considerably [ 1,2]. Among other methods, the finite element models have also been employed in more general cases of completely flexible mechanisms. In the finite element method permanent contact between bodies is not expressed by kinematic constraints, but by letting them have boundary points in common. In order to apply the finite element method to a linkage mechanism, the articulating system is generally modelled as a series of instantaneous structures. Thus the continuous motion of the system is replaced by a sequence of structures at discrete crank angle configurations. The obtained equations of motion are usually solved by adopting the well-known harmonic series technique [3]. This analysis requires the inversion of the system stiffness matrix, which for mechanisms is singular. To avoid this complication the rigid body degrees of freedom of the linkage must be removed from the model. In order to do this, usually the crank is considered as a cantilever beam [4] or a matrix decomposition is used [5]. The method presented in this paper is formulated for planar mechanisms with revolute joints. It is assumed that: (i) the nodes are revolute joints with no friction; (ii) externally applied forces are imposed to the joints; and (iii) the speed and accelerations of the mechanism are moderate. Due to the last assumption the effect of inertia torques on the dynamic response of the mechanism (i.e. bending of links) can be neglected and every link of the mechanism can be modelled as a bar in the prime state of tension or compression. Such uniaxial stresses are characteristic
2. FINITE ELEMENT FORMULATION
The finite element description is developed for planar mechanisms with revolute joints. Making an assumption that the material of column elements is linearly elastic (follows Hooke’s law) with uniform mass density p and uniform cross-sectional area A,
135
136
M
HAC
the element stiffness matrix [k], of the ith element is in the same form as for a planar truss element [9]
or
in the local coordinate is}:=
(1) where E is the Young’s modulus and 1, is the length of the ith finite element. A general planar truss element is shown in Fig. la in three frames of reference. The XY frame is the global (in&ally fixed) reference frame, ,V.Vis a relative frame and 51 is the local reference frame. The .XY frame is attached to the node (joint) of a truss element at the beginning of each time step discretization. Its axes are parallel to the global reference frame. The nodal displacements in the xy and XY frames are essentially the same during each discretized time interval. The displacement between the XY and .XYframes is calculated in every step of time discretization. The local reference frame is such that the g-axis is parallel to the truss element axis. It is assumed that during the movement the ith element of the mechanism has changed its position determined by endpoints 1 and 2 to the points I’ and 2’ (dotted). The components of nodal displacement vectors 11’ and 22’ can be expressed in the relative coordinate system SJ by the following vector:
system
{q,. t,.q:. l>],.
The vector {a}, is transformed
where [r], is a transformation r cos CI, -sin 2, u-1, = o
cos 2, o
lo
0
sin r,
(J
O
() cos x,
O sin 2,
-sin r,
cos r, I
7
(5)
where 2, is the angle between local and global coordinate systems for the ith element. Introducing the transformation matrices [T,], and [TJ,, one has
,T,
],
=
iu )! = (41%4>Ir=[7‘,l,ld;,
(6)
It),={t,.r,if=[r~],;d),
(7)
L ‘0;
a,
-sin a, o
si;xi
0
0
cos 2,
sin r,
cos z, o
I
0 O -sin g, COSCX, I
(8)
(9)
Treating the finite element presented in Fig. la as a part of a planar mechanism, it should be noted that due to the movement of the mechanism the angle ~1, is changing and in this way transformation matrices depend on the current position of the given
I
11
mechanism
link. 3. EQUATIONS
1
I
(4) matrix
[
0
into {dj, by
{6j,=[T],jd;.,.
[~*I, =
(a’“4
(3)
Wl
OF MOTION
X
*
X
The equations of motion of a mechanism can be derived from Lagrange’s equations. The energy dissipation and body forces are not taken into account. The strain energy U,, of the ith element is based on a one-dimensional stress analysis and is in the same form as for planar coordinates it is as follows:
trusses [9]. In global
u,.,=~jd)JtT,l~[kl,tT,l,{d),.
Fig. I. (a) Displacements of a finite element in different coordinate frames. (b) Generalized displacement: of a finite element’s point.
(10)
To obtain a kinetic energy T,, of the mechanism element both longitudinal and transverse displacements of the element are taken into account. Let G be any point on link i, as shown in Fig. lb, with the position determined by vector rO in the inertial frame XY. After deformation the point G takes the place denoted by G’ and may be located in the global coordinate system by vector r. From Fig. lb the position of the point G’ can be written as r = r,, + r .
(11)
137
Dynamics of planar flexible mechanisms r’ is a deformation vector in the local coordinate system &. For the velocity of the point G’, eqn (11) is differentiated with respect to time to obtain where
i = v,~+ w x r’ -i- r’,
[M,J is the inertia matrix
(12)
where v, is a velocity vector of point G in the inertial frame XI’, w = col[O, 0, w] is the absolute angular velocity of the & frame, and a dot over a variable denotes differentiation with respect to time. The vector product in eqn (12) can be expressed as 0 x r’ = w[Q]r’,
(13)
where w is the length of the vector o, and matrix [QJ is defined as
0 -1 PI = 1 () [ I The kinetic energy expressed as follows:
where the coefficient matrices of the ith element in local coordinates are as given below.
[IL,],
is the centrifugal matrix
[C,,], is the gyroscopic matrix (arises due to the Coriolis component of the longitudinal acceleration of the link in a rotational velocity field)
(14)
Tei of the ith element
is
{Ffi is known in purely rotational centifugal force vector
motion
as the
and the vector (Dj, is given by =- :
’ (v:v,, + i’*i’ + w*r’r[Q]Tf2]r’ I0
+
2wi’r@]r’+2wr’~R]rv,,+2i’T~,,)dm.
(15)
Intr~ucing a matrix notation essential in the finite element formulation, the displacement vector r’ can be expressed in terms of the displacements of the element’s endpoints using the following formula (the subscript i indicates that the given term corresponds to the ith element of a mechanism): V)i =
k%I,~~Ii.
where W,= dai/dt is the angular velocity of the ith finite element. The kinetic energy T, of the ith element can also be expressed in terms of the planar truss element inertia matrix [ml, by substituting the following formulae into eqn (18):
(16)
Taking into account that all components of vectors (6}i change linearly with the length of the finite element, the shape function matrix [N,], is given as
1 -(iiri) @
O (14 0 l -(Ciri) O ii/i
1 (17
where
’
where 0 < [ C I,. Note that in our case the shape function represents both longitudinal and transverse displacements of the element. Employing eqns (16) and (4) in eqn (Is), and omitting a term not dependent on generalized coordinates, the kinetic energy 2’+,i in the global coordinate system is given by
(18)
(26)
Substituting expressions on energies of all elements into Lagrange’s equations, the equations of motion are obtained in the form
where [Ml, [C,], [K], and [A;‘,]are assembled forms of the element matrices [A&&, [C&Ii, f/c], and [lyceli, respectively, (P> is the vector of external nodal loads, {FJ is an assembled form of vectors {F),, and {Z}, {zt}, {x} represent acceleration velocity and
138
M.
displacement vectors (in nodal points). Vector {x} is an assembled form of vectors {d}, with regard to boundary conditions. In deriving eqn (27) from eqn (18) the constancy of the matrices [M,],, [C,,], and vector {D}, has been implicitly assumed. This is consistent with the common assumption in numerical methods that the timedependent system parameters are held constant over each discretized time interval. Making use of Lagrange’s equations the properties of the coefficient matrices have been taken into account, i.e. the matrices [Ml, [K]. and [Kc] are symmetric, and [C,] is antisymmetric. Contrary to the classical finite element analysis, the coefficient matrices in eqn (27) are functions of linkage geometry and vary as the input (crank) angle is varied. A somewhat similar situation occurs in the finite element analysis of heat conduction problems where the heat conductivity matrix depends on time [lo]. The solution {x} of the matrix eqn (27) is obtained by using a finite difference method. Each matrix is approximated by its finite difference with respect to time t. Thus the matrix eqn (27) is discretized for successive periods of time interval At. Introducing a notation {u 1 = {.t}, eqn (27) becomes a system of two equations, which for time t + (Afj2) (midpoint rule) are as follows:
HAM
step of time ACltis complete when the difference between the former and the latest estimation of {x},,~, is lower than the given tolerance. It should be noted that the force vector {P) contains externally applied forces imposed on the nodes only in the first step of time discretization. In the next calculations the force vector contains also generalized stresses corresponding to the elongation of flexible column elements calculated in the previous step of time discretization (these forces become external in the current step of calculation). For the ith element they can be expressed in the form (in local coordinates)
where the superscripts (n) and (n - 1) indicate current step of time discretization. 4. NUMERICAL
EXAMPLE
As an example, the kinematics and dynamics of a four-bar linkage are presented. The mechanism is divided into three finite elements (Fig. 2). For planar motion the nodal coordinates are elements of an eight-dimensional vector space. Taking into account the boundary conditions, the vector of the unknown functions (nodal displacement vector {.x}) consists of four elements: displacements M’and v in the X and Y directions, respectively, of two moving nodes 2 and 3:
’ I r = [M.z) L’?.11’3, I.J]. 1x1
M, = M,,, + M,, sin 2,.
(.ri,,
=i(iPj,+i\r + {P),, +f(jF), /A,+ {F),)
(31)
The input torque M, is applied to the crank of the mechanism and is assumed to be
+ ~([C:‘l,,~I+I(.ll,)(i-~~ird,-(lli)
- m, 1,+A,+K,l,)(l.~j,*,,+
the
(29)
with appropriate initial conditions for {x) and (u). The subscripts in eqns (28) and (29) indicate the time at which the corresponding term must be evaluated. by an iteration The solution {x}, + A, is evaluated method. The first iteration of the solution {.~}j”,‘~, must be estimated. The estimation can be very rough and may be easily determined by predicting a direction of the movement of a mechanism. The scheme of the iteration procedure for solving eqn (27) is as follows: for initial conditions imposed on ix), and {u}, the matrices [Ml,, [C,],, [K],, [K,],, and {F), are determined as assembled forms of appropriate matrices in eqns (19). (21), (I), (20). and (22). Then the solution {x}!? *, is estimated and again from the referenced equations the first estimations of [M]j’: *, , [C,]lOJA,, [K]jOJ*, , [Kc]:“: Ar, and {F}I”i A, are obtained. Next, by using the system of eqns (28) and (29), The iteration of one lullOjar and ix Z’i dl are obtained.
(32)
The force vector (Pi consists of two components: the first is the externally applied forces imposed on node 2 (due to the existence of the input torque M,) and the second is the generalized forces P,, [eqn (30)] corresponding to the elongation of elements. The coefficient matrices of eqn (27) are calculated in every step of time discretization. The structure stiffness [K],
Y
t
1
Fig. 2. Four-bar
linkage.
139
Dynamics of planar flexible mechanisms inertia [Ml, and gyroscopic [C,] matrices are given in Appendix A. Solving the system of equations (28) and
(29) the nodal displacement vector {x} and velocity vector {a} are obtained for successive periods of time. To estimate the accuracy of the obtained results, a kinematic and dynamic analysis of the four-bar linkage with rigid links is made. In this case the system has one degree of freedom and the crank’s angle displacement A, is taken as an independent variable (in the rigid body model variables are taken as capital Greek letters of the appropriate variables for the flexible link mechanism). The angles A,, Aj, Fig. 4. The velocity i, of nodal point 2 with time for flexible (-) and rigid (---) links. the angular velocities R, , R2, R,, and angular accelerations c,, t2, t, of links 1, 2, and 3, respectively, describe the rigid-body motion of the mechanism [l 11. From these figures it is seen that the motion of the The governing equations are derived by using mechanism is slower for the flexible body model. Lagrange’s equations of motion and can be written as Bt, = C,
5. CONCLUSIONS
(33)
where B and C are functions of the linkage geometry, angular velocities and accelerations, and are given in Appendix B. The differential equation (33) is solved numerically by using a standard procedure based on the Runge-Kutta method. The input data for kinematic and dynamic analysis of a four-bar linkage fulfills the Grashof category. The mechanism is made of steel. Its geometric and physical properties are as follows. Length of links (Fig. 2): crank 1, = 1 m, coupler fr = 2 m, rocker I3 = 2.5 m, ground l., = 3 m; cross-sectional area of links A = lo-’ m2, density p = 7800 kg/m3, elastic modulus E = 2.1 x 10” N/m’, input torque [eqn (32)] M,,, = 5 Nm, M,, = 2 Nm. The initial conditions in link notation are: A, = 1.1 rad, the rigid R, = 0 rad/sec. Simulation time is from 0 to 1.85 sec. During this time the crank of the mechanism is able to do a full rotation. The displacement and velocity in the X direction of the nodal point 2 of the mechanism are presented in Figs 3 and 4, respectively. These figures contain graphs for flexible links as well as for rigid
The analysis presented in this paper has two main aims. Firstly, a procedure based on the finite difference technique has been employed to solve the equations of motion of flexible mechanisms. The advantages of this method are as follows. (a) It avoids the complication connected with the singularity of the system stiffness matrix. (b) The variant input torque of any form can be very easy incorporated in the solution procedure. Secondly, the truss-type finite element has been applied for modelling moderate speed mechanisms. Since the truss-type finite element does not inhabit rotary freedoms at its nodes, a special shape function has been built to properly model a revolute joint of a mechanism. The knowledge of the rigid body motion is not required. In view of the presented results the method is considered reliable and the differences between the flexible body model and the rigid model are acceptable. REFERENCES
links.
1. G. G. Lowen and C. Chassapis, The elastic behavior of linkages: an update. Mech. Mach. Theory 21, 3342 (1986). 2. B. S. Thompson and C. K. Sung, A survey of finite element technics for mechanism design. Mech. Much. Theory 21, 351-359 (1986). 3. B. M. Bahgat and K. D. Willmert, Finite element vibrational analysis of planar mechanisms. Mech. Much. Theory 11, 47-71 (1976). 4. A. Midha, A. G. Erdman and D. A. Frohrib, Finite element approach to mathematical modeling of highspeed elastic linkages. Mech. Much. Theory 13, 603 (1978).
Fig. 3. The displacement wr of nodal point 2 with time for flexible (-) and rigid (---) links.
P. K. Nath and A. Ghosh,
Kineto-elastodynamic analysis of mechanisms by finite element method. Mech. Much. Theory 15, 179-197 (1980). 6. J. F. Besseling, J. B. Jonker and A. L. Schwab, Kinematics and dynamics of mechanisms. Delfr Prog. Rep. 10, 16&172 (1985). 5.
140
M. HAG
7. J. 0. Song and E. J. Haug, Dynamic analysis of planar flexible mechanism. Comput. Meth. appl. Mech. Engng 24, 359-381 (1980). 8. X. Gao, Z. King and Q. Zhang, A hybrid beam element for mathematical modelling of high-speed flexible linkages. Mech. Mach. Theory 24, 29 (1989). 9. C. S. Desai, Elementary Finire Element Method. Prentice-Hall, Englewood Cliffs, NJ (1979). 10. S. Orivuori, Efficient method for solution of nonlinear heat conduction problems. In/. J. Numer. Meth. Engng 14, 1461-1476 (1979). 1 I. H. M. Mabie and F. W. Ocvirk. Mechanisms and Dynamics of Machiner,v. John Wiley, New York (1966).
[ac,+c,) [M]=pA/6
O c2
0
0
I
Symm. 2(c, fc,)
0
[CR1 = WA I
‘3
The coefficients A B
Introducing
the following
=
notation
a,=(I,)~“cosu,
I
+
2l2l
2(Czf(.d
0
0
- (‘2
c2
0
0 2(1,,+c,)
0
B
B and C in eon (33) are as follows:
sin’(A, -- A,) sfn’(A,- Ad f 6ljl, + 21; I22-.-.----’ ‘sln2(A,-A2) sin (A, - Al) sin(A, - ~4,)
+ 61;12 cos(A, - A,) --;--- --~ -I~stn(A, - .4,)
h, = (I,)-’ 2 sin a, c,=lf(af+hy),
213
0 Antisymm.
2(c,+c,)
APPENDIX APPENDIX
2(cz+c,)
(‘7
sin(A, - .4,) sin(A, - A,) C = 6M,,‘p + 2Sl:l, 7-m- --~- - 2l;[,R 7---m-sln(A, - ,4,) s!n(A, - .4:)
i = 1,2,3,
the stiffness [K], inertia [M] and gyroscopic global coordinates are as follows:
[C,] matrices
- 31,I;R - 31,l;R cos(A; -- AZ).
in where
R = [I&
cos(A, - A,) + I@, + lfn, cos(A, - A,)]/[12 sin(A, - A:)]
(K] = EA
S = [l,Q; + l$2: cos(A, - A?) + I, Qy cos(A, - A2)]/[13sin(A, - A,)].
1991
0
MMS-7949/91 $3.00 + 0.00 1991 Pergmon Press plc
DYNAMICS OF PLANAR FLEXIBLE MECHANISMS BY FINITE ELEMENT METHOD WITH TRUSS-TYPE ELEMENTS M. HAM Institute of Machine Design Fundamentals, Warsaw Technical University, Narbutta 84, 02-524 Warszawa, Poland (Received 26 February 1990)
Abstract-The kinematics and dynamics of planar mechanisms with revolute joints and flexible links are studied by using the finite element method. The mechanism finite element is based on a planar truss element. The equations of motion of a mechanism are formulated by applying Lagrange’s equations to the kinetic and potential energy of the system. The analysis includes the higher order elastic acceleration terms. The presented method of solving equations of motion is based on a finite difference procedure and avoids the necessity of elimination of the singularity of the system stiffness matrix. A four-bar linkage is studied to illustrate the computational efficiency of the presented method. The results for flexible links are compared with those obtained for the rigid body.
1. INTRODUCTION
elements. In order to properly model a revolute joint of a mechanism? a distinct rotary freedom should be assigned to each of the elements connected by the node that replaces the revolute joint. This has been done by assuming not only longitudinal but also transverse displacement of the element nodes (i.e. it is possible for the element to move and rotate in any direction, but the stresses occur only due to the elongation of the element). Lagrange’s equations are used to obtain the equations of motion of a mechanism. The coefficient matrices are functions of the linkage geometry and vary as the input (crank) angle is varied. In order to solve the equations of motion the finite difference method with respect to time is used. In this method the problem of singularity of the system stiffness matrix does not exist. Contrary to most of the other publications, the variable input torque is considered here. As an example a kinematic and dynamic analysis of a four-bar linkage is presented. The element stiffness matrix is based on the planar truss element. It is asumed that maximum stresses occurring in links do not exceed critical values. Such an approach for the modelling of four-bar linkages was suggested in [6]. Other methodologists [4, 7,8] used a series of beam elements to model a four-bar linkage. for truss
In the past few years the number of publications which deal with elastic behaviour of mechanisms has increased considerably [ 1,2]. Among other methods, the finite element models have also been employed in more general cases of completely flexible mechanisms. In the finite element method permanent contact between bodies is not expressed by kinematic constraints, but by letting them have boundary points in common. In order to apply the finite element method to a linkage mechanism, the articulating system is generally modelled as a series of instantaneous structures. Thus the continuous motion of the system is replaced by a sequence of structures at discrete crank angle configurations. The obtained equations of motion are usually solved by adopting the well-known harmonic series technique [3]. This analysis requires the inversion of the system stiffness matrix, which for mechanisms is singular. To avoid this complication the rigid body degrees of freedom of the linkage must be removed from the model. In order to do this, usually the crank is considered as a cantilever beam [4] or a matrix decomposition is used [5]. The method presented in this paper is formulated for planar mechanisms with revolute joints. It is assumed that: (i) the nodes are revolute joints with no friction; (ii) externally applied forces are imposed to the joints; and (iii) the speed and accelerations of the mechanism are moderate. Due to the last assumption the effect of inertia torques on the dynamic response of the mechanism (i.e. bending of links) can be neglected and every link of the mechanism can be modelled as a bar in the prime state of tension or compression. Such uniaxial stresses are characteristic
2. FINITE ELEMENT FORMULATION
The finite element description is developed for planar mechanisms with revolute joints. Making an assumption that the material of column elements is linearly elastic (follows Hooke’s law) with uniform mass density p and uniform cross-sectional area A,
135
136
M
HAC
the element stiffness matrix [k], of the ith element is in the same form as for a planar truss element [9]
or
in the local coordinate is}:=
(1) where E is the Young’s modulus and 1, is the length of the ith finite element. A general planar truss element is shown in Fig. la in three frames of reference. The XY frame is the global (in&ally fixed) reference frame, ,V.Vis a relative frame and 51 is the local reference frame. The .XY frame is attached to the node (joint) of a truss element at the beginning of each time step discretization. Its axes are parallel to the global reference frame. The nodal displacements in the xy and XY frames are essentially the same during each discretized time interval. The displacement between the XY and .XYframes is calculated in every step of time discretization. The local reference frame is such that the g-axis is parallel to the truss element axis. It is assumed that during the movement the ith element of the mechanism has changed its position determined by endpoints 1 and 2 to the points I’ and 2’ (dotted). The components of nodal displacement vectors 11’ and 22’ can be expressed in the relative coordinate system SJ by the following vector:
system
{q,. t,.q:. l>],.
The vector {a}, is transformed
where [r], is a transformation r cos CI, -sin 2, u-1, = o
cos 2, o
lo
0
sin r,
(J
O
() cos x,
O sin 2,
-sin r,
cos r, I
7
(5)
where 2, is the angle between local and global coordinate systems for the ith element. Introducing the transformation matrices [T,], and [TJ,, one has
,T,
],
=
iu )! = (41%4>Ir=[7‘,l,ld;,
(6)
It),={t,.r,if=[r~],;d),
(7)
L ‘0;
a,
-sin a, o
si;xi
0
0
cos 2,
sin r,
cos z, o
I
0 O -sin g, COSCX, I
(8)
(9)
Treating the finite element presented in Fig. la as a part of a planar mechanism, it should be noted that due to the movement of the mechanism the angle ~1, is changing and in this way transformation matrices depend on the current position of the given
I
11
mechanism
link. 3. EQUATIONS
1
I
(4) matrix
[
0
into {dj, by
{6j,=[T],jd;.,.
[~*I, =
(a’“4
(3)
Wl
OF MOTION
X
*
X
The equations of motion of a mechanism can be derived from Lagrange’s equations. The energy dissipation and body forces are not taken into account. The strain energy U,, of the ith element is based on a one-dimensional stress analysis and is in the same form as for planar coordinates it is as follows:
trusses [9]. In global
u,.,=~jd)JtT,l~[kl,tT,l,{d),.
Fig. I. (a) Displacements of a finite element in different coordinate frames. (b) Generalized displacement: of a finite element’s point.
(10)
To obtain a kinetic energy T,, of the mechanism element both longitudinal and transverse displacements of the element are taken into account. Let G be any point on link i, as shown in Fig. lb, with the position determined by vector rO in the inertial frame XY. After deformation the point G takes the place denoted by G’ and may be located in the global coordinate system by vector r. From Fig. lb the position of the point G’ can be written as r = r,, + r .
(11)
137
Dynamics of planar flexible mechanisms r’ is a deformation vector in the local coordinate system &. For the velocity of the point G’, eqn (11) is differentiated with respect to time to obtain where
i = v,~+ w x r’ -i- r’,
[M,J is the inertia matrix
(12)
where v, is a velocity vector of point G in the inertial frame XI’, w = col[O, 0, w] is the absolute angular velocity of the & frame, and a dot over a variable denotes differentiation with respect to time. The vector product in eqn (12) can be expressed as 0 x r’ = w[Q]r’,
(13)
where w is the length of the vector o, and matrix [QJ is defined as
0 -1 PI = 1 () [ I The kinetic energy expressed as follows:
where the coefficient matrices of the ith element in local coordinates are as given below.
[IL,],
is the centrifugal matrix
[C,,], is the gyroscopic matrix (arises due to the Coriolis component of the longitudinal acceleration of the link in a rotational velocity field)
(14)
Tei of the ith element
is
{Ffi is known in purely rotational centifugal force vector
motion
as the
and the vector (Dj, is given by =- :
’ (v:v,, + i’*i’ + w*r’r[Q]Tf2]r’ I0
+
2wi’r@]r’+2wr’~R]rv,,+2i’T~,,)dm.
(15)
Intr~ucing a matrix notation essential in the finite element formulation, the displacement vector r’ can be expressed in terms of the displacements of the element’s endpoints using the following formula (the subscript i indicates that the given term corresponds to the ith element of a mechanism): V)i =
k%I,~~Ii.
where W,= dai/dt is the angular velocity of the ith finite element. The kinetic energy T, of the ith element can also be expressed in terms of the planar truss element inertia matrix [ml, by substituting the following formulae into eqn (18):
(16)
Taking into account that all components of vectors (6}i change linearly with the length of the finite element, the shape function matrix [N,], is given as
1 -(iiri) @
O (14 0 l -(Ciri) O ii/i
1 (17
where
’
where 0 < [ C I,. Note that in our case the shape function represents both longitudinal and transverse displacements of the element. Employing eqns (16) and (4) in eqn (Is), and omitting a term not dependent on generalized coordinates, the kinetic energy 2’+,i in the global coordinate system is given by
(18)
(26)
Substituting expressions on energies of all elements into Lagrange’s equations, the equations of motion are obtained in the form
where [Ml, [C,], [K], and [A;‘,]are assembled forms of the element matrices [A&&, [C&Ii, f/c], and [lyceli, respectively, (P> is the vector of external nodal loads, {FJ is an assembled form of vectors {F),, and {Z}, {zt}, {x} represent acceleration velocity and
138
M.
displacement vectors (in nodal points). Vector {x} is an assembled form of vectors {d}, with regard to boundary conditions. In deriving eqn (27) from eqn (18) the constancy of the matrices [M,],, [C,,], and vector {D}, has been implicitly assumed. This is consistent with the common assumption in numerical methods that the timedependent system parameters are held constant over each discretized time interval. Making use of Lagrange’s equations the properties of the coefficient matrices have been taken into account, i.e. the matrices [Ml, [K]. and [Kc] are symmetric, and [C,] is antisymmetric. Contrary to the classical finite element analysis, the coefficient matrices in eqn (27) are functions of linkage geometry and vary as the input (crank) angle is varied. A somewhat similar situation occurs in the finite element analysis of heat conduction problems where the heat conductivity matrix depends on time [lo]. The solution {x} of the matrix eqn (27) is obtained by using a finite difference method. Each matrix is approximated by its finite difference with respect to time t. Thus the matrix eqn (27) is discretized for successive periods of time interval At. Introducing a notation {u 1 = {.t}, eqn (27) becomes a system of two equations, which for time t + (Afj2) (midpoint rule) are as follows:
HAM
step of time ACltis complete when the difference between the former and the latest estimation of {x},,~, is lower than the given tolerance. It should be noted that the force vector {P) contains externally applied forces imposed on the nodes only in the first step of time discretization. In the next calculations the force vector contains also generalized stresses corresponding to the elongation of flexible column elements calculated in the previous step of time discretization (these forces become external in the current step of calculation). For the ith element they can be expressed in the form (in local coordinates)
where the superscripts (n) and (n - 1) indicate current step of time discretization. 4. NUMERICAL
EXAMPLE
As an example, the kinematics and dynamics of a four-bar linkage are presented. The mechanism is divided into three finite elements (Fig. 2). For planar motion the nodal coordinates are elements of an eight-dimensional vector space. Taking into account the boundary conditions, the vector of the unknown functions (nodal displacement vector {.x}) consists of four elements: displacements M’and v in the X and Y directions, respectively, of two moving nodes 2 and 3:
’ I r = [M.z) L’?.11’3, I.J]. 1x1
M, = M,,, + M,, sin 2,.
(.ri,,
=i(iPj,+i\r + {P),, +f(jF), /A,+ {F),)
(31)
The input torque M, is applied to the crank of the mechanism and is assumed to be
+ ~([C:‘l,,~I+I(.ll,)(i-~~ird,-(lli)
- m, 1,+A,+K,l,)(l.~j,*,,+
the
(29)
with appropriate initial conditions for {x) and (u). The subscripts in eqns (28) and (29) indicate the time at which the corresponding term must be evaluated. by an iteration The solution {x}, + A, is evaluated method. The first iteration of the solution {.~}j”,‘~, must be estimated. The estimation can be very rough and may be easily determined by predicting a direction of the movement of a mechanism. The scheme of the iteration procedure for solving eqn (27) is as follows: for initial conditions imposed on ix), and {u}, the matrices [Ml,, [C,],, [K],, [K,],, and {F), are determined as assembled forms of appropriate matrices in eqns (19). (21), (I), (20). and (22). Then the solution {x}!? *, is estimated and again from the referenced equations the first estimations of [M]j’: *, , [C,]lOJA,, [K]jOJ*, , [Kc]:“: Ar, and {F}I”i A, are obtained. Next, by using the system of eqns (28) and (29), The iteration of one lullOjar and ix Z’i dl are obtained.
(32)
The force vector (Pi consists of two components: the first is the externally applied forces imposed on node 2 (due to the existence of the input torque M,) and the second is the generalized forces P,, [eqn (30)] corresponding to the elongation of elements. The coefficient matrices of eqn (27) are calculated in every step of time discretization. The structure stiffness [K],
Y
t
1
Fig. 2. Four-bar
linkage.
139
Dynamics of planar flexible mechanisms inertia [Ml, and gyroscopic [C,] matrices are given in Appendix A. Solving the system of equations (28) and
(29) the nodal displacement vector {x} and velocity vector {a} are obtained for successive periods of time. To estimate the accuracy of the obtained results, a kinematic and dynamic analysis of the four-bar linkage with rigid links is made. In this case the system has one degree of freedom and the crank’s angle displacement A, is taken as an independent variable (in the rigid body model variables are taken as capital Greek letters of the appropriate variables for the flexible link mechanism). The angles A,, Aj, Fig. 4. The velocity i, of nodal point 2 with time for flexible (-) and rigid (---) links. the angular velocities R, , R2, R,, and angular accelerations c,, t2, t, of links 1, 2, and 3, respectively, describe the rigid-body motion of the mechanism [l 11. From these figures it is seen that the motion of the The governing equations are derived by using mechanism is slower for the flexible body model. Lagrange’s equations of motion and can be written as Bt, = C,
5. CONCLUSIONS
(33)
where B and C are functions of the linkage geometry, angular velocities and accelerations, and are given in Appendix B. The differential equation (33) is solved numerically by using a standard procedure based on the Runge-Kutta method. The input data for kinematic and dynamic analysis of a four-bar linkage fulfills the Grashof category. The mechanism is made of steel. Its geometric and physical properties are as follows. Length of links (Fig. 2): crank 1, = 1 m, coupler fr = 2 m, rocker I3 = 2.5 m, ground l., = 3 m; cross-sectional area of links A = lo-’ m2, density p = 7800 kg/m3, elastic modulus E = 2.1 x 10” N/m’, input torque [eqn (32)] M,,, = 5 Nm, M,, = 2 Nm. The initial conditions in link notation are: A, = 1.1 rad, the rigid R, = 0 rad/sec. Simulation time is from 0 to 1.85 sec. During this time the crank of the mechanism is able to do a full rotation. The displacement and velocity in the X direction of the nodal point 2 of the mechanism are presented in Figs 3 and 4, respectively. These figures contain graphs for flexible links as well as for rigid
The analysis presented in this paper has two main aims. Firstly, a procedure based on the finite difference technique has been employed to solve the equations of motion of flexible mechanisms. The advantages of this method are as follows. (a) It avoids the complication connected with the singularity of the system stiffness matrix. (b) The variant input torque of any form can be very easy incorporated in the solution procedure. Secondly, the truss-type finite element has been applied for modelling moderate speed mechanisms. Since the truss-type finite element does not inhabit rotary freedoms at its nodes, a special shape function has been built to properly model a revolute joint of a mechanism. The knowledge of the rigid body motion is not required. In view of the presented results the method is considered reliable and the differences between the flexible body model and the rigid model are acceptable. REFERENCES
links.
1. G. G. Lowen and C. Chassapis, The elastic behavior of linkages: an update. Mech. Mach. Theory 21, 3342 (1986). 2. B. S. Thompson and C. K. Sung, A survey of finite element technics for mechanism design. Mech. Much. Theory 21, 351-359 (1986). 3. B. M. Bahgat and K. D. Willmert, Finite element vibrational analysis of planar mechanisms. Mech. Much. Theory 11, 47-71 (1976). 4. A. Midha, A. G. Erdman and D. A. Frohrib, Finite element approach to mathematical modeling of highspeed elastic linkages. Mech. Much. Theory 13, 603 (1978).
Fig. 3. The displacement wr of nodal point 2 with time for flexible (-) and rigid (---) links.
P. K. Nath and A. Ghosh,
Kineto-elastodynamic analysis of mechanisms by finite element method. Mech. Much. Theory 15, 179-197 (1980). 6. J. F. Besseling, J. B. Jonker and A. L. Schwab, Kinematics and dynamics of mechanisms. Delfr Prog. Rep. 10, 16&172 (1985). 5.
140
M. HAG
7. J. 0. Song and E. J. Haug, Dynamic analysis of planar flexible mechanism. Comput. Meth. appl. Mech. Engng 24, 359-381 (1980). 8. X. Gao, Z. King and Q. Zhang, A hybrid beam element for mathematical modelling of high-speed flexible linkages. Mech. Mach. Theory 24, 29 (1989). 9. C. S. Desai, Elementary Finire Element Method. Prentice-Hall, Englewood Cliffs, NJ (1979). 10. S. Orivuori, Efficient method for solution of nonlinear heat conduction problems. In/. J. Numer. Meth. Engng 14, 1461-1476 (1979). 1 I. H. M. Mabie and F. W. Ocvirk. Mechanisms and Dynamics of Machiner,v. John Wiley, New York (1966).
[ac,+c,) [M]=pA/6
O c2
0
0
I
Symm. 2(c, fc,)
0
[CR1 = WA I
‘3
The coefficients A B
Introducing
the following
=
notation
a,=(I,)~“cosu,
I
+
2l2l
2(Czf(.d
0
0
- (‘2
c2
0
0 2(1,,+c,)
0
B
B and C in eon (33) are as follows:
sin’(A, -- A,) sfn’(A,- Ad f 6ljl, + 21; I22-.-.----’ ‘sln2(A,-A2) sin (A, - Al) sin(A, - ~4,)
+ 61;12 cos(A, - A,) --;--- --~ -I~stn(A, - .4,)
h, = (I,)-’ 2 sin a, c,=lf(af+hy),
213
0 Antisymm.
2(c,+c,)
APPENDIX APPENDIX
2(cz+c,)
(‘7
sin(A, - .4,) sin(A, - A,) C = 6M,,‘p + 2Sl:l, 7-m- --~- - 2l;[,R 7---m-sln(A, - ,4,) s!n(A, - .4:)
i = 1,2,3,
the stiffness [K], inertia [M] and gyroscopic global coordinates are as follows:
[C,] matrices
- 31,I;R - 31,l;R cos(A; -- AZ).
in where
R = [I&
cos(A, - A,) + I@, + lfn, cos(A, - A,)]/[12 sin(A, - A:)]
(K] = EA
S = [l,Q; + l$2: cos(A, - A?) + I, Qy cos(A, - A2)]/[13sin(A, - A,)].