- Email: [email protected]
Nonlinear dynamics analysis of flexible beams under large overall motions and the flexible manipulator simulation
Compurers & S~rurrur~~sVol. 56. No. I, pp I-14. 1995 Copyright 3 1995 Elsever Science Ltd Printed in Great Britain. All rights reserved 0045.7949195 $9.50 + 0.00
00457949(94)00542-7
NONLINEAR DYNAMICS ANALYSIS OF FLEXIBLE BEAMS UNDER LARGE OVERALL MOTIONS AND THE FLEXIBLE MANIPULATOR SIMULATION J. L. Meek and Hua Liu Department
of Civil Engineering,
The University (Received
of Queensland,
28 February
Queensland
4072, Australia
1994)
Abstract-The dynamic model of a flexible Timoshenko beam with geometrical nonlinearities subject to large overall motions has been derived by using the finite element method and multi-body system formulation. The equations of motion are derived using Lagrange’s equation based on the expression of the kinetic and potential energies of the flexible beam in terms of generalized coordinates. The nonlinear constraint equations are adjoined to the system equations of motion using Lagrange multipliers. An incremental-iterative method based on the Newmark direct integration together with Newton-Raphson iteration is employed here for the numerical studies. Numerical examples are presented to demonstrate the efficiency of the present method.
(1) reference coordinates which describe large rigidbody translation and large angular rotations of the body reference; and (2) elastic coordinates which characterize elastic deformations of the body, that is, relative translational and angular displacement of infinitesimal volumes at nodal points on the body. The location and angular orientation of every infinitesimal volume in each element can be approximated in terms of its elastic coordinates and the assumed shape functions by the finite element approach. Energy equations of each element are first written separately using the assumed shape functions. The elements of each body are then assembled using a standard finite element procedure. The degrees of freedom of each body may be reduced using a sub-structuring technique. Algebraic equations prescribing constraints between various bodies are formulated and coupled to the equations of motion by a Lagrange multiplier technique. An implicit numerical integration method is used to solve the resulting set of nodal equilibrium equations. The intention of this study is to produce analysis software suitable for a broad range of applications in contrast with other approaches mentioned above. Formulations and software for the in-plane case using the proposed approach have been implemented. Cubic functions for translational deformation and quadratic functions for angular rotation are used as shape functions which lead to results of both mass and stiffness matrices consistent with the pure bending theory. The Newmark implicit integration method is employed in combination with the Newton-Raphson iteration method to handle the nonlinear nature of the resulting dynamic equations. Some examples, including planar beam models of a manipulator, are presented to demonstrate the
1. INTRODUCTION dynamic analysis of flexible bodies has been the subject of many investigations. One of the fields of application is robotics, in which earlier work has concentrated on a strictly rigid-body model for a manipulator system. The demands for higher speeds and better system performance have caused the investigation of lightweight designs which require that the flexibility be included because of the increased deformations due to the bending and torsion effects of the flexible manipulator. A method for dynamic analysis of manipulators with flexible members is needed for the design and control of such systems. Currently, the most popular approach for this analysis is to develop a finite element model. In the finite element method the flexible members of a manipulator are modeled by beam elements. The resulting equations of motion are usually nonlinear and highly coupled in the inertia terms due to the presence of Coriolis and centrifugal effects as well as the inertia due to rotation of the beam. The inertia properties and eigenspectrum vary with time. This makes the problem different from general structural dynamics. Various approaches have been developed recently by several authors, such as the float frame method, inertial frame method, co-rotational formulation, etc. [I, 21. Unfortunately the majority of these studies have concentrated only on specific problems and are not easily employed for solution of more general and complex systems. Herein an alternative approach based on a multibody system formulation [3] and small deflection Timoshenko beam theory with the inclusion of the geometrical nonlinearity and the shear deformation is proposed. Two sets of coordinates have been defined:
The
1
J. L. Meek and Hua
Fig. 1. The configuration
efficiency of the proposed formulation. The last example is of a two link robot arm, in this case with only translational coupling between the links.
Liu
of deformation.
where
In this section the procedure for deriving the equations of motion is described. Firstly, the means of expressing the deformations of a body are summarized, and then, using the finite element formulation, the equations of motion are obtained.
Let the position of a material point P measured with respect to a reference axes, X = [X, X2 XjIT, be denoted by coordinates x = [x, x2 x3]’ in the undeformed state, as shown in Fig. 1, while for the material point P’ in the deformed geometry the coordinates 5 = [[, t2 &IT measured relative to the reference X are employed. The two coordinates are related using the displacement field u = [u, u2 us]‘, as follows, using the Lagrange description,
(1) d< = J dx,
where J is the Jacobian of the coordinate ation which can be written as
(2)
J=I+u,
Then substituting
(5)
ds: = dxT dx.
(6)
of the line element,
S,,
(7)
eqns (5) and (6) into eqn (7) S, = dxT[C ’ + c “‘]dx
(8)
in which 6’ = f(u: + ud), 6”’ = f(uldd). The superscripts 1, nl denote the linear and nonlinear components respectively. The total strain, t is written as t =C’+t”‘. By writing the total strain vector components order
the deferential (4)
ds; = dtTd<
Se = f(ds; - ds;).
ET = From eqn (1), we obtain
1
u3.2 u3.3
An indication of the extension can be written as
transform-
(3)
4.3
u2.2 u2.3
From Fig. 1, the lengths of the line element in the deformed and undeformed states, respectively, can be written by the expressions,
2.1. Large deformation analysis
5 = c(x);
4.1
L%I
2. FLEXIBLE BODY SYSTEM
u=(-x
4.I f4.2
% =
k,l~22~33~,2~234
operators
(9) in the
(10)
D, and D, give,
t = (0, + D,)u
(11)
Nonhear
dynamics
analysis of flexible beams
where
a
zj II:=
0
0
a dx,
0
0
1;
0 0 a
0
;;
-3
Ia2
Ia
-2 ax,
_2 dx,
0
0
Ia -2 ax2
Ia -2 ax,
dx,
(12)
and a ‘I.1 ax,
a 'I.2 ax,
a
2u
h.2
T$
','ax,
a
a u2.1 dx,
'I.3 a ax3
u2,2 ax,
2%
g
3
3
a u2,3 dx,
2UZ.l &
2u2.2 $
a
a u3.2 ax,
a '3.3 dx,
u3J ax,
Then
3
a2 2U3,, ax*
2u3.2 &
2u3.1 & 3
2.2. Finite element formulation E =Du
whereD=D,+D,. Using interpolation ment field
(14)
Consider a deformed
u=Nq where iV is an appropriate
assume
a displace-
(15)
set of shape function and coordinates. In the finite element method the shape functions represent a piecewise polynomial that describes the displacements over each element in the body. Now the total strain vector can be written as,
q is a set of time-dependent
L =
mq.
3
of the jkxible
bodies
the elastic deformation of the point P on body i, which can be written as u; = N’q;.
functions,
(13)
2U2.l gx3
(17)
In eqn (17) N’ is an appropriate shape function, and q,.is a set of time dependent elastic coordinates of body i. Using the finite element method, the body i is divided into a set of small elements. P is in jth element so that the point is denoted, Pg. As shown in Fig. 2, select (x, y, z) to be an inertial frame of reference, and (x’, y’, zi) to be a coordinate system of the ith body. Let R’ and 19’represent translational and rotational orientation of the (x’, y’, 9) frame with respect to the (x, y, t) system. This is called the set of reference coordinates.
Fig. 2. The coordinate
systems.
J. L. Meek and Hua Liu
4 q: = [P
PjT.
(18)
Substitute q’=
Let du and eg be the deformed and undeformed vectors that define the location of the point P with respect to the body (x’,y’, z’) coordinate system. Then U” is elastic deformation at point P, such that, d” =
[&”
eqn
(23) into eqn (24) Then the kinetic can be written as
BIT cj;TIT.
expression
where
e; + N”ql.
(19) LN
The global location
and let energy
of the point
P is given,
”
BTLN NTLT
r; = R’ + Lie8 + L’N”q;.
NTLTBT
(20)
dV”
NTN I
Where L’ is the transformation matrix from the body coordinate system to the inertial frame. The last term in eqn (20) represents the contribution from consideration of the flexibility of the body. or partitioned
2.3. Energy expressions Differentiation of eqn (20) once with respect time yields the velocity vector:
=
This equation as
L;d”
=
$ .
&I
=
$
(el
.
+
u”),
can also be written in partitioned
If the ith body has n elements,
4jT]‘, then
m,,(s) m&) ’4: I[ 1
to
m,(q)
where B”
as 4’= [4:
(22) form
(23)
m,,(q)
4;
(27)
In eqn (27) q; are the reference coordinates locating the x’, y’. z’ body fixed coordinate system with respect to the inertial frame, whereas the variable q; represents the elastic generalized coordinates of the body. Note that the off diagonal terms, m;, = m;, are the nonlinear matrices representing the inertia coupling between the reference and flexible coordinates. This, nonlinear nature of in addition to mrrr represents inertia terms in the kinetic energy expression, The strain energy expression of the ith body can be written as
then
(28) (24)
An isotropic
linear material
is assumed,
Neutral axis /
Fig. 3. Two-dimensional
beam element.
fu5
thus
Nonlinear dynamics analysis of flexible beams o’j=
,,&”
(29)
then eqn (28) can be written
U’ = 1 Substitute
as
c “TEQ ‘id v”. s “0
In eqn (37), F represents the quadratic velocity term that contains the Coriolis components and the gyroscopic force. The eqn (36) can be rewritten as
(30)
eqn (16) into eqn (30) so that, U’ = f9F@
(31)
=[;;]+[;]+[;;]A(38)
where
F? = c
[NTBTEbN]
dV
(32)
is a stiffness matrix which represents the sum of the conventional elastic stiffness K: and the geometric stiffness K!! resulting from including the quadratic terms of strains in the strain energy expression to account for large displacements.
3.
A TWO-DIMENSIONAL BEAM ELEMENT WITH SHEAR DEFORMATION AND GEOMETRIC NONLINEARITIES
The special case of the two-dimensional lation is derived in this section.
3.1. Element 2.4. The system equations of motion The composite vector of all system coordinates is denoted by: q =
[qlT
qZT
.
The nonlinear
constraint
of the bodies).
equations
@(q, t) = 0 (vector Using the Lagrange
=
equations
The dimensions of a beam element are shown in Fig. 3. Define the element coordinates corresponding to the body reference,
-$
(33)
(34)
The elastic displacement of any point on the neutral axis of element j is given by
(35)
(40)
of motion
matrix and 1 is it follows that
0
0 rl(l+pLr()l-_)) 0
configuration
can be written,
form).
where @J: = (a/84)@’ is the Jacobian the vector of Lagrange multipliers,
rl
generalized
formu-
qm~T
(N = number
FJ
5
0
&l+a(?-0) ?(l -3Pcr)
The shape function
5
0
0
{(l-W(P-+c-l+P(a-0)
0
(36) (37)
+tl
I?l, can be written,
0
! ‘I
.
<(I - 3W)
(41)
In eqn (41) 5 =x/l; rj = 1 - 5; p = l/(1 + 121); 1 = (EI)/(C12); and C is shear stiffness, which is equal to K . A . G, where K is Timoshenko shear coefficient and G is the shear modulus.
J. L. Meek and Hua Liu
6
A cubic function is used for tiJ and a quadratic one for 0: because the consistent results with pure bending (non shear deformation) can be made if i = 0 and p = 1 when G+co, c-+00.
Substituting
3.2. Elastic st@?kss and mass matrices The strain energy written in the form:
for the beam
element
I
u;i =
0
-y’-
01
(46)
0
eqn (40) into eqn (46) yields
can be
U; = AI?q:! = Nq;! where a new shape function
(47)
can be obtained
N = AN!‘.
(48)
Denoting the shape function for the neutral the element defined in eqn (41) as (42)
NV=
NI N [1 2
(49)
LKJ
Substituting eqn (40) into eqn (42) the strain energy expression can be obtained
axis of
leads to ,‘J’/ =
$ (q:T~;q’/
+
where Kt is linear stiffness
Ki,
_ c
q”T~;qt/)
matrix
(43)
Al 2rp2
0
0
0
6
31
0
31 (2-t 61)12
2Erp2 s
Al
-2
and Kg the geometric
stiffness
I
-6
0
31
matrix
0
0
2IP 0
-6
31
0
-31
(1 - 61)12
1
(50)
w 0
-31
0
6
-31
(1 - 61)12
0
- 31
(2 + 61)12
given by 0
0
6p21
0
o0 0
-6O-6p21 0 12p2
(5 + -6~~1 03/l2)12
0
-6~~1 60 +0 12p2
0
6p21
0
-6~~1
IO
Iv, - yiv, I$ [
1
Al
x0 21/l 2
0
0
2Iv 0
-i
O 60 + 12~~
r”
K;=p,
N=
given by
(-5+3/w
0
-6O-
0
12p2
6p21 (-5+3pY -6~~1 0
(45)
(5 + 3/w
where Substitute
this into eqn (26) to obtain
EA (ts, - C,) Pg =
6012
According to the beam theories, the displacement at any point in the element due to the elastic deformation u$= [u,~u,JT can be written as:
+ rT?Il fl:m, S’0
dr.
(51)
Nonlinear dynamics analysis of flexible beams Then the element
mass matrix
7
can be obtained 0
L+K_’
0
PI I P21( i
P212 I2
120
840+120
o
my,= til
420
24
120 420 24 p21= I’ --840 120
O
0
o
PI I
IP’
1
120
420
24
0
P212 I2 --840
0
0
/II - -----
l/l’
I
120
420
24
0
0
-;+!g
0
5
-f+?&
--
PI2
+
2
+5’
0 0 0
0
3/l?= -+j 10
0
12
/lI 31u= --_ 2 5 -- p12 3y=1= 2+10 -+z
12
12
0
0
0
!$
0
pl --2
f_!$
3.3. Integration
defines the undeformed sume
L=
cos 8,
[ sm 8,
-sin
0,
cos 8,
1
where u0 = [u,, position
90 = [xi
%=[’
1
-sin tIi
- cos 8,
cos 6,
-sin
I
u,],
uoJIT
of the element.
41T
8,
ti1,r
(53)
L
ir
;][;;]=%q,,.
(54) As-
(55)
SITLT
(56)
Now recalling eqn (26), the components of the mass matrix can be derived in the body reference:
f[““o*]+SJq,]
LS'
[[*fqT + qTs#T] L'T -y’ + [2qi + q:1hf;, q/g + q:ls'
1
+C 3
then
and
M’=
5 3/l?? 10
B = L’[Nq,+
In the planar case the transformation matrix between the global reference coordinates and the body fixed coordinates is
(52) 31u=
-$+-
From eqn (22),
matrix
-$+-+f 3j.S 0
The global stiffness and mass matrix can be derived by summing all the element matrices using standard finite element routines. of body mass and stzfness
7
10
0
-- 6~= 5 -- p’I= 3/l?= 2 +- 1o+z
f_Ag
0
-f+$-
6~~
--
5 2
840+120
0
2
-6P2
0
0
120
pl lp’ I -----120 420 24 /1212 12
s’T[40+ S/l
4,
1 -I
(57)
8
J.
L.
Meek and Hua Liu The second term of the forces in eqn (37) can also be obtained as
where x,. is the centre point of the beam measured relative to the body fixed coordinates; I, is the length of the beam; NE is the number of the element in beam i;
0 -4O-2/.1
s+
-51 -/.Ll 0 -20+2/l 51 - /.ll
0
1
(58)
40 + 2/l
51+ PI
0
0
0 20 - 2p
0 51- PI
- [i’L(S,
=
+ Mfl(q0 + q#
&sTL’Ti
0
0
0
0
and C is a transformation matrix between the element coordinates and body fixed coordinates. The body stiffness matrix is given as
0
20-2/l
-20+2/l -51+/L1 0
0
0 40 + 2/l
0 -51 -/lI
0
0
0
0
-40-2~ 51-t /ll
.
(63)
+ ST&] I
1
-51+$
0
(59)
Then, the total inertial forces due to the quadratic velocity terms will be F = F, + F,.
4. NUMERICAL
where
+ Sq,)8 - PL’sq/]
INTEGRATION
PROCEDURES
Newmark’s direct integration method has been used to integrate this nonlinear system equation: r(q) = M(q)q + C(q)4 + K(q)q - F(q, 4, t)
and
- Q(t) = 0 (W K, = E C=K;C. ,=I
(60) subject
to the nonlinear
constraints
3.4. Inertial forces due to the quadratic velocity terms Consider
@(q, t) = 0.
eqn (37). The first term can be written as F, = -A&j.
In the planar
(61)
case, the expression
for 1%4is derived,
0
n;/=
qTIO - L[Sb + S’,,] L’S’lT QIOL’S’O]T
(65)
Here it has been assumed that the damping matrix C is proportional to the mass and/or the stiffness matrix.
[O - L [Sb + S’q& 2q;DTM;,g,
+ 2q;kf;,q,
PTcj,
[OL’S’ O]Q 4fSi 0
. 1
(62)
Nonlinear
dynamics
analysis
of flexible beams
I *
0
EI=l,OOO
EA=lOO
I=10
m=l
Fig. 4. A flexible robot
arm model.
The integration procedures are as follows. (1) Predict acceleration, velocities and displacement at the time corresponding to step n + 1 by:
or, making
. 4n+r =qn 4/i+1=~Jf+(l %+I =qn+h4,+fh2(l
(67)
-2B)&inh28&+,,
(68)
where h is the time interval, and gration parameters. For Newmark’s j3 = l/4 in the average acceleration jl = l/6 in the linear acceleration (2) Obtain corrections to the Newton type iteration: S&n+,
y and fi are intemethod, y = l/2, case and y = l/2, case. displacements by
where S is the tangent iteration metric. In order to include the constraint equation the following derivations are made. A virtual displacement 6q is considered which verifies that:
-aq1
a@ T 6q=O.
[
(70)
(71)
Thus the Lagrange multipliers 1 may be introduced which transfer the constraint problem eqns (56) and (57) into the augmented system of n + m equations, where m is the number of constraint equations
r(q)+1
[
Ff a4
T
=o
1
@((I,t) = 0.
(6%
= r(q,+s)
.dq
use of eqn (57):
(66)
-~)~&l++M+,
1sq=o
'a@T
@(q + &7, t) = @(q, t) +
(72)
This system is then solved incrementally by a predicted solution 4 and assuming an approximate solution 1. From these values, (4, /Ir), an improved solution (q + Aq, /r + AA) is calculated. A first order expansion of eqn (64) yields the linear equation
5.0
Fig. 5. A given rotation
angle.
J. L. Meek and Hua Liu
10
Aq+i
=O
(73) The correction vector [Aq Ai]’ is the solution augmented tangent equation
Fig. 6. A displacement
driven
with the tangent
iteration
matrix
and the extended
residual
vector
of the
flexible beam (a)
L 1 r(q) +
g(& 8 =
In Newmark’s
1
E [ aq1
4
@(a
(76)
method
(77) ~=W+~C(q)+~Wq)-$and the iterated
Fig. 7. A displacement
driven
formulation
flexible beam (b).
with the convergence
criterion (79)
2.5
Fig. 8. A given driven
force.
Nonlinear
dynamics
analysis
of flexible beams
11
Fig. 9. A force driven flexible beam. velocities and (3) Modify the accelerations, displacements at time t,, ,: the following are used to obtain the variables at the end of a new time step: 1
(b) Force driven. The beam with the same properties as above is now driven by a prescribed torque T(t) in Fig. 8 applied at the axis of rotation. The simulation is shown in Fig. 9.
q.+r=pAy”+,
(81) ?.+I
=qn+Aqn+,.
(82)
5. NUMERICAL EXAMPLES
Here numerical examples which simulate the flexible robot arms are given. These are taken from Refs [l, 21, and results may be compared with the sources. 5.1. Flexible robot arm Consider the rotation of a flexible beam about a horizontal axis passing through one end of the beam (Fig. 4) to simulate the motion of a flexible robot arm. Five elements are used. The material properties are: EI = 1000, EA = 100, m = 1, and K =0.867. Two cases (a) and (b) are considered. (a) Displacement driven. The beam is driven by a given rotation angle $(t) shown in Fig. 5. The simulation results are given in Figs 6 and 7.
EL500
EA=lOOO
m=l
Fig. 10. A flexible rod.
J. L. Meek and Hua
Liu
k8.0
Fig.
11,The first two revolutions
of the flying rod.
b
0.5
t
Fig. 12. A given rotation t=1.4
angle.
t=1.5
t=f .6
-8
‘_ t=0.7
t=0.6
k2.1
t=0.4
\
k0.2
t=0.3 Fig. 13. Displacement
driven
multi-component
robot
arm.
t=0.1
Nonlinear
160
dynamics
analysis
of flexible beams
13
.
0.5
t
F(t)=T(t)/4
Fig. 14. The given driving
Fig. 15. A force driven
multibody
force.
system
in free flight.
J. L. Meek and Hua
14 5.2.
Fl_ving jlexible
rod
A flexible rod with free ends, initially inclined position shown in Fig. IO, is force and a torque applied simultaneously The time history of T(t) is shown in F(t) = T(t)/lO. The first two revolutions Fig, I I 5.3. A displacement arm
driven
placed in an subject to a at one end. Fig. 8, and are shown in
multi-component
robot
The robot arm in example I is now considered to be composed of two flexible components connected together by a hinge at mid-length. The two components were subjected to the prescribed rotation as shown in Fig. 12. The material properties are: EI = 100, EA = 10000, m = I, K = 0.867 for both links. The sequence of motion is shown in Fig. 13. Five elements are used in each links. 5.4. A jorce dricen multibody
s,vstem in free flight
A two body system consisting of two flexible links connected by a hinge is initially at an inclined position. The system is set into motion by applying a force and a torque at one end of the lower link, as shown in Fig. 14. The material properties are the same as in example 3. Subsequently the articulated beams undergo free flight that is shown in Fig. 15. 6. CONCLUDING
REMARKS
In the formulation used herein all Coriolis and centrifugal effects, as well as the other inertia effects due to rotation, are automatically accounted for. The simulation program can be used directly for the
Liu
dynamic analysis of a multi-link flexible robot arm in high speed motion. This approach can be carried over into the ful! three-dimensional case. This extension relies on a proper treatment of three-dimensional finite rotations in both the structural deformations of the beam and in the overall motions. All the simulation results reported here were calculated by a FORTRAN program named BODYFLEX written by the second author. This program has been tested on a few UNIX based systems (SUN-OS, NEWS-OS and ESIX). Graphics presentation was made through X-windows and a PostScript printer by IMSL Exponent Graphics system.
REFERENCES
1. J. C. Simo and L. Vu-Quot. beams under large overall 849-863 (1986). 2. K. M. Hsiao and J. Y. Jang, flexible mechanisms by Compul.
Mefh.
On the dynamics of flexible motions. J. appl. M&I. 53, Dynamic analysis of planar co-rotational formulation.
Appi. Mech. und Engng 87, IL14 (1991).
nonlin3. E. M. Bakr and A. A. Shabana, Geometrically ear analysis of multi-body systems. Compur. Struck. 23(6), 7399751 (1986). 4. J. Meek and H. Liu, Nonlinear dynamic analysis of flexible beams under large motions. In Compurutional Mechanics (Edited by Valhappan, Pulmano and TinLoi), Vol. 1, pp. 1755182. Balkema. Rotterdam (1993). modelling and simu5. J. Meek and H. Liu. Nonlinear lation of a flexible manipulator. Proc. lsr Pun-Pacific Conf. on Computational Engineering. Seoul. Korea (1993). Analytical Dvuzmics of Discrete 6. R. M. Rosenberg, Systems. Plenum Press, New York (1987). 1. J. Argyris and H. P. Mlejnek, D_wamics of structure. North-Holland, Amsterdam (1991).
00457949(94)00542-7
NONLINEAR DYNAMICS ANALYSIS OF FLEXIBLE BEAMS UNDER LARGE OVERALL MOTIONS AND THE FLEXIBLE MANIPULATOR SIMULATION J. L. Meek and Hua Liu Department
of Civil Engineering,
The University (Received
of Queensland,
28 February
Queensland
4072, Australia
1994)
Abstract-The dynamic model of a flexible Timoshenko beam with geometrical nonlinearities subject to large overall motions has been derived by using the finite element method and multi-body system formulation. The equations of motion are derived using Lagrange’s equation based on the expression of the kinetic and potential energies of the flexible beam in terms of generalized coordinates. The nonlinear constraint equations are adjoined to the system equations of motion using Lagrange multipliers. An incremental-iterative method based on the Newmark direct integration together with Newton-Raphson iteration is employed here for the numerical studies. Numerical examples are presented to demonstrate the efficiency of the present method.
(1) reference coordinates which describe large rigidbody translation and large angular rotations of the body reference; and (2) elastic coordinates which characterize elastic deformations of the body, that is, relative translational and angular displacement of infinitesimal volumes at nodal points on the body. The location and angular orientation of every infinitesimal volume in each element can be approximated in terms of its elastic coordinates and the assumed shape functions by the finite element approach. Energy equations of each element are first written separately using the assumed shape functions. The elements of each body are then assembled using a standard finite element procedure. The degrees of freedom of each body may be reduced using a sub-structuring technique. Algebraic equations prescribing constraints between various bodies are formulated and coupled to the equations of motion by a Lagrange multiplier technique. An implicit numerical integration method is used to solve the resulting set of nodal equilibrium equations. The intention of this study is to produce analysis software suitable for a broad range of applications in contrast with other approaches mentioned above. Formulations and software for the in-plane case using the proposed approach have been implemented. Cubic functions for translational deformation and quadratic functions for angular rotation are used as shape functions which lead to results of both mass and stiffness matrices consistent with the pure bending theory. The Newmark implicit integration method is employed in combination with the Newton-Raphson iteration method to handle the nonlinear nature of the resulting dynamic equations. Some examples, including planar beam models of a manipulator, are presented to demonstrate the
1. INTRODUCTION dynamic analysis of flexible bodies has been the subject of many investigations. One of the fields of application is robotics, in which earlier work has concentrated on a strictly rigid-body model for a manipulator system. The demands for higher speeds and better system performance have caused the investigation of lightweight designs which require that the flexibility be included because of the increased deformations due to the bending and torsion effects of the flexible manipulator. A method for dynamic analysis of manipulators with flexible members is needed for the design and control of such systems. Currently, the most popular approach for this analysis is to develop a finite element model. In the finite element method the flexible members of a manipulator are modeled by beam elements. The resulting equations of motion are usually nonlinear and highly coupled in the inertia terms due to the presence of Coriolis and centrifugal effects as well as the inertia due to rotation of the beam. The inertia properties and eigenspectrum vary with time. This makes the problem different from general structural dynamics. Various approaches have been developed recently by several authors, such as the float frame method, inertial frame method, co-rotational formulation, etc. [I, 21. Unfortunately the majority of these studies have concentrated only on specific problems and are not easily employed for solution of more general and complex systems. Herein an alternative approach based on a multibody system formulation [3] and small deflection Timoshenko beam theory with the inclusion of the geometrical nonlinearity and the shear deformation is proposed. Two sets of coordinates have been defined:
The
1
J. L. Meek and Hua
Fig. 1. The configuration
efficiency of the proposed formulation. The last example is of a two link robot arm, in this case with only translational coupling between the links.
Liu
of deformation.
where
In this section the procedure for deriving the equations of motion is described. Firstly, the means of expressing the deformations of a body are summarized, and then, using the finite element formulation, the equations of motion are obtained.
Let the position of a material point P measured with respect to a reference axes, X = [X, X2 XjIT, be denoted by coordinates x = [x, x2 x3]’ in the undeformed state, as shown in Fig. 1, while for the material point P’ in the deformed geometry the coordinates 5 = [[, t2 &IT measured relative to the reference X are employed. The two coordinates are related using the displacement field u = [u, u2 us]‘, as follows, using the Lagrange description,
(1) d< = J dx,
where J is the Jacobian of the coordinate ation which can be written as
(2)
J=I+u,
Then substituting
(5)
ds: = dxT dx.
(6)
of the line element,
S,,
(7)
eqns (5) and (6) into eqn (7) S, = dxT[C ’ + c “‘]dx
(8)
in which 6’ = f(u: + ud), 6”’ = f(uldd). The superscripts 1, nl denote the linear and nonlinear components respectively. The total strain, t is written as t =C’+t”‘. By writing the total strain vector components order
the deferential (4)
ds; = dtTd<
Se = f(ds; - ds;).
ET = From eqn (1), we obtain
1
u3.2 u3.3
An indication of the extension can be written as
transform-
(3)
4.3
u2.2 u2.3
From Fig. 1, the lengths of the line element in the deformed and undeformed states, respectively, can be written by the expressions,
2.1. Large deformation analysis
5 = c(x);
4.1
L%I
2. FLEXIBLE BODY SYSTEM
u=(-x
4.I f4.2
% =
k,l~22~33~,2~234
operators
(9) in the
(10)
D, and D, give,
t = (0, + D,)u
(11)
Nonhear
dynamics
analysis of flexible beams
where
a
zj II:=
0
0
a dx,
0
0
1;
0 0 a
0
;;
-3
Ia2
Ia
-2 ax,
_2 dx,
0
0
Ia -2 ax2
Ia -2 ax,
dx,
(12)
and a ‘I.1 ax,
a 'I.2 ax,
a
2u
h.2
T$
','ax,
a
a u2.1 dx,
'I.3 a ax3
u2,2 ax,
2%
g
3
3
a u2,3 dx,
2UZ.l &
2u2.2 $
a
a u3.2 ax,
a '3.3 dx,
u3J ax,
Then
3
a2 2U3,, ax*
2u3.2 &
2u3.1 & 3
2.2. Finite element formulation E =Du
whereD=D,+D,. Using interpolation ment field
(14)
Consider a deformed
u=Nq where iV is an appropriate
assume
a displace-
(15)
set of shape function and coordinates. In the finite element method the shape functions represent a piecewise polynomial that describes the displacements over each element in the body. Now the total strain vector can be written as,
q is a set of time-dependent
L =
mq.
3
of the jkxible
bodies
the elastic deformation of the point P on body i, which can be written as u; = N’q;.
functions,
(13)
2U2.l gx3
(17)
In eqn (17) N’ is an appropriate shape function, and q,.is a set of time dependent elastic coordinates of body i. Using the finite element method, the body i is divided into a set of small elements. P is in jth element so that the point is denoted, Pg. As shown in Fig. 2, select (x, y, z) to be an inertial frame of reference, and (x’, y’, zi) to be a coordinate system of the ith body. Let R’ and 19’represent translational and rotational orientation of the (x’, y’, 9) frame with respect to the (x, y, t) system. This is called the set of reference coordinates.
Fig. 2. The coordinate
systems.
J. L. Meek and Hua Liu
4 q: = [P
PjT.
(18)
Substitute q’=
Let du and eg be the deformed and undeformed vectors that define the location of the point P with respect to the body (x’,y’, z’) coordinate system. Then U” is elastic deformation at point P, such that, d” =
[&”
eqn
(23) into eqn (24) Then the kinetic can be written as
BIT cj;TIT.
expression
where
e; + N”ql.
(19) LN
The global location
and let energy
of the point
P is given,
”
BTLN NTLT
r; = R’ + Lie8 + L’N”q;.
NTLTBT
(20)
dV”
NTN I
Where L’ is the transformation matrix from the body coordinate system to the inertial frame. The last term in eqn (20) represents the contribution from consideration of the flexibility of the body. or partitioned
2.3. Energy expressions Differentiation of eqn (20) once with respect time yields the velocity vector:
=
This equation as
L;d”
=
$ .
&I
=
$
(el
.
+
u”),
can also be written in partitioned
If the ith body has n elements,
4jT]‘, then
m,,(s) m&) ’4: I[ 1
to
m,(q)
where B”
as 4’= [4:
(22) form
(23)
m,,(q)
4;
(27)
In eqn (27) q; are the reference coordinates locating the x’, y’. z’ body fixed coordinate system with respect to the inertial frame, whereas the variable q; represents the elastic generalized coordinates of the body. Note that the off diagonal terms, m;, = m;, are the nonlinear matrices representing the inertia coupling between the reference and flexible coordinates. This, nonlinear nature of in addition to mrrr represents inertia terms in the kinetic energy expression, The strain energy expression of the ith body can be written as
then
(28) (24)
An isotropic
linear material
is assumed,
Neutral axis /
Fig. 3. Two-dimensional
beam element.
fu5
thus
Nonlinear dynamics analysis of flexible beams o’j=
,,&”
(29)
then eqn (28) can be written
U’ = 1 Substitute
as
c “TEQ ‘id v”. s “0
In eqn (37), F represents the quadratic velocity term that contains the Coriolis components and the gyroscopic force. The eqn (36) can be rewritten as
(30)
eqn (16) into eqn (30) so that, U’ = f9F@
(31)
=[;;]+[;]+[;;]A(38)
where
F? = c
[NTBTEbN]
dV
(32)
is a stiffness matrix which represents the sum of the conventional elastic stiffness K: and the geometric stiffness K!! resulting from including the quadratic terms of strains in the strain energy expression to account for large displacements.
3.
A TWO-DIMENSIONAL BEAM ELEMENT WITH SHEAR DEFORMATION AND GEOMETRIC NONLINEARITIES
The special case of the two-dimensional lation is derived in this section.
3.1. Element 2.4. The system equations of motion The composite vector of all system coordinates is denoted by: q =
[qlT
qZT
.
The nonlinear
constraint
of the bodies).
equations
@(q, t) = 0 (vector Using the Lagrange
=
equations
The dimensions of a beam element are shown in Fig. 3. Define the element coordinates corresponding to the body reference,
-$
(33)
(34)
The elastic displacement of any point on the neutral axis of element j is given by
(35)
(40)
of motion
matrix and 1 is it follows that
0
0 rl(l+pLr()l-_)) 0
configuration
can be written,
form).
where @J: = (a/84)@’ is the Jacobian the vector of Lagrange multipliers,
rl
generalized
formu-
qm~T
(N = number
FJ
5
0
&l+a(?-0) ?(l -3Pcr)
The shape function
5
0
0
{(l-W(P-+c-l+P(a-0)
0
(36) (37)
+tl
I?l, can be written,
0
! ‘I
.
<(I - 3W)
(41)
In eqn (41) 5 =x/l; rj = 1 - 5; p = l/(1 + 121); 1 = (EI)/(C12); and C is shear stiffness, which is equal to K . A . G, where K is Timoshenko shear coefficient and G is the shear modulus.
J. L. Meek and Hua Liu
6
A cubic function is used for tiJ and a quadratic one for 0: because the consistent results with pure bending (non shear deformation) can be made if i = 0 and p = 1 when G+co, c-+00.
Substituting
3.2. Elastic st@?kss and mass matrices The strain energy written in the form:
for the beam
element
I
u;i =
0
-y’-
01
(46)
0
eqn (40) into eqn (46) yields
can be
U; = AI?q:! = Nq;! where a new shape function
(47)
can be obtained
N = AN!‘.
(48)
Denoting the shape function for the neutral the element defined in eqn (41) as (42)
NV=
NI N [1 2
(49)
LKJ
Substituting eqn (40) into eqn (42) the strain energy expression can be obtained
axis of
leads to ,‘J’/ =
$ (q:T~;q’/
+
where Kt is linear stiffness
Ki,
_ c
q”T~;qt/)
matrix
(43)
Al 2rp2
0
0
0
6
31
0
31 (2-t 61)12
2Erp2 s
Al
-2
and Kg the geometric
stiffness
I
-6
0
31
matrix
0
0
2IP 0
-6
31
0
-31
(1 - 61)12
1
(50)
w 0
-31
0
6
-31
(1 - 61)12
0
- 31
(2 + 61)12
given by 0
0
6p21
0
o0 0
-6O-6p21 0 12p2
(5 + -6~~1 03/l2)12
0
-6~~1 60 +0 12p2
0
6p21
0
-6~~1
IO
Iv, - yiv, I$ [
1
Al
x0 21/l 2
0
0
2Iv 0
-i
O 60 + 12~~
r”
K;=p,
N=
given by
(-5+3/w
0
-6O-
0
12p2
6p21 (-5+3pY -6~~1 0
(45)
(5 + 3/w
where Substitute
this into eqn (26) to obtain
EA (ts, - C,) Pg =
6012
According to the beam theories, the displacement at any point in the element due to the elastic deformation u$= [u,~u,JT can be written as:
+ rT?Il fl:m, S’0
dr.
(51)
Nonlinear dynamics analysis of flexible beams Then the element
mass matrix
7
can be obtained 0
L+K_’
0
PI I P21( i
P212 I2
120
840+120
o
my,= til
420
24
120 420 24 p21= I’ --840 120
O
0
o
PI I
IP’
1
120
420
24
0
P212 I2 --840
0
0
/II - -----
l/l’
I
120
420
24
0
0
-;+!g
0
5
-f+?&
--
PI2
+
2
+5’
0 0 0
0
3/l?= -+j 10
0
12
/lI 31u= --_ 2 5 -- p12 3y=1= 2+10 -+z
12
12
0
0
0
!$
0
pl --2
f_!$
3.3. Integration
defines the undeformed sume
L=
cos 8,
[ sm 8,
-sin
0,
cos 8,
1
where u0 = [u,, position
90 = [xi
%=[’
1
-sin tIi
- cos 8,
cos 6,
-sin
I
u,],
uoJIT
of the element.
41T
8,
ti1,r
(53)
L
ir
;][;;]=%q,,.
(54) As-
(55)
SITLT
(56)
Now recalling eqn (26), the components of the mass matrix can be derived in the body reference:
f[““o*]+SJq,]
LS'
[[*fqT + qTs#T] L'T -y’ + [2qi + q:1hf;, q/g + q:ls'
1
+C 3
then
and
M’=
5 3/l?? 10
B = L’[Nq,+
In the planar case the transformation matrix between the global reference coordinates and the body fixed coordinates is
(52) 31u=
-$+-
From eqn (22),
matrix
-$+-+f 3j.S 0
The global stiffness and mass matrix can be derived by summing all the element matrices using standard finite element routines. of body mass and stzfness
7
10
0
-- 6~= 5 -- p’I= 3/l?= 2 +- 1o+z
f_Ag
0
-f+$-
6~~
--
5 2
840+120
0
2
-6P2
0
0
120
pl lp’ I -----120 420 24 /1212 12
s’T[40+ S/l
4,
1 -I
(57)
8
J.
L.
Meek and Hua Liu The second term of the forces in eqn (37) can also be obtained as
where x,. is the centre point of the beam measured relative to the body fixed coordinates; I, is the length of the beam; NE is the number of the element in beam i;
0 -4O-2/.1
s+
-51 -/.Ll 0 -20+2/l 51 - /.ll
0
1
(58)
40 + 2/l
51+ PI
0
0
0 20 - 2p
0 51- PI
- [i’L(S,
=
+ Mfl(q0 + q#
&sTL’Ti
0
0
0
0
and C is a transformation matrix between the element coordinates and body fixed coordinates. The body stiffness matrix is given as
0
20-2/l
-20+2/l -51+/L1 0
0
0 40 + 2/l
0 -51 -/lI
0
0
0
0
-40-2~ 51-t /ll
.
(63)
+ ST&] I
1
-51+$
0
(59)
Then, the total inertial forces due to the quadratic velocity terms will be F = F, + F,.
4. NUMERICAL
where
+ Sq,)8 - PL’sq/]
INTEGRATION
PROCEDURES
Newmark’s direct integration method has been used to integrate this nonlinear system equation: r(q) = M(q)q + C(q)4 + K(q)q - F(q, 4, t)
and
- Q(t) = 0 (W K, = E C=K;C. ,=I
(60) subject
to the nonlinear
constraints
3.4. Inertial forces due to the quadratic velocity terms Consider
@(q, t) = 0.
eqn (37). The first term can be written as F, = -A&j.
In the planar
(61)
case, the expression
for 1%4is derived,
0
n;/=
qTIO - L[Sb + S’,,] L’S’lT QIOL’S’O]T
(65)
Here it has been assumed that the damping matrix C is proportional to the mass and/or the stiffness matrix.
[O - L [Sb + S’q& 2q;DTM;,g,
+ 2q;kf;,q,
PTcj,
[OL’S’ O]Q 4fSi 0
. 1
(62)
Nonlinear
dynamics
analysis
of flexible beams
I *
0
EI=l,OOO
EA=lOO
I=10
m=l
Fig. 4. A flexible robot
arm model.
The integration procedures are as follows. (1) Predict acceleration, velocities and displacement at the time corresponding to step n + 1 by:
or, making
. 4n+r =qn 4/i+1=~Jf+(l %+I =qn+h4,+fh2(l
(67)
-2B)&inh28&+,,
(68)
where h is the time interval, and gration parameters. For Newmark’s j3 = l/4 in the average acceleration jl = l/6 in the linear acceleration (2) Obtain corrections to the Newton type iteration: S&n+,
y and fi are intemethod, y = l/2, case and y = l/2, case. displacements by
where S is the tangent iteration metric. In order to include the constraint equation the following derivations are made. A virtual displacement 6q is considered which verifies that:
-aq1
a@ T 6q=O.
[
(70)
(71)
Thus the Lagrange multipliers 1 may be introduced which transfer the constraint problem eqns (56) and (57) into the augmented system of n + m equations, where m is the number of constraint equations
r(q)+1
[
Ff a4
T
=o
1
@((I,t) = 0.
(6%
= r(q,+s)
.dq
use of eqn (57):
(66)
-~)~&l++M+,
1sq=o
'a@T
@(q + &7, t) = @(q, t) +
(72)
This system is then solved incrementally by a predicted solution 4 and assuming an approximate solution 1. From these values, (4, /Ir), an improved solution (q + Aq, /r + AA) is calculated. A first order expansion of eqn (64) yields the linear equation
5.0
Fig. 5. A given rotation
angle.
J. L. Meek and Hua Liu
10
Aq+i
=O
(73) The correction vector [Aq Ai]’ is the solution augmented tangent equation
Fig. 6. A displacement
driven
with the tangent
iteration
matrix
and the extended
residual
vector
of the
flexible beam (a)
L 1 r(q) +
g(& 8 =
In Newmark’s
1
E [ aq1
4
@(a
(76)
method
(77) ~=W+~C(q)+~Wq)-$and the iterated
Fig. 7. A displacement
driven
formulation
flexible beam (b).
with the convergence
criterion (79)
2.5
Fig. 8. A given driven
force.
Nonlinear
dynamics
analysis
of flexible beams
11
Fig. 9. A force driven flexible beam. velocities and (3) Modify the accelerations, displacements at time t,, ,: the following are used to obtain the variables at the end of a new time step: 1
(b) Force driven. The beam with the same properties as above is now driven by a prescribed torque T(t) in Fig. 8 applied at the axis of rotation. The simulation is shown in Fig. 9.
q.+r=pAy”+,
(81) ?.+I
=qn+Aqn+,.
(82)
5. NUMERICAL EXAMPLES
Here numerical examples which simulate the flexible robot arms are given. These are taken from Refs [l, 21, and results may be compared with the sources. 5.1. Flexible robot arm Consider the rotation of a flexible beam about a horizontal axis passing through one end of the beam (Fig. 4) to simulate the motion of a flexible robot arm. Five elements are used. The material properties are: EI = 1000, EA = 100, m = 1, and K =0.867. Two cases (a) and (b) are considered. (a) Displacement driven. The beam is driven by a given rotation angle $(t) shown in Fig. 5. The simulation results are given in Figs 6 and 7.
EL500
EA=lOOO
m=l
Fig. 10. A flexible rod.
J. L. Meek and Hua
Liu
k8.0
Fig.
11,The first two revolutions
of the flying rod.
b
0.5
t
Fig. 12. A given rotation t=1.4
angle.
t=1.5
t=f .6
-8
‘_ t=0.7
t=0.6
k2.1
t=0.4
\
k0.2
t=0.3 Fig. 13. Displacement
driven
multi-component
robot
arm.
t=0.1
Nonlinear
160
dynamics
analysis
of flexible beams
13
.
0.5
t
F(t)=T(t)/4
Fig. 14. The given driving
Fig. 15. A force driven
multibody
force.
system
in free flight.
J. L. Meek and Hua
14 5.2.
Fl_ving jlexible
rod
A flexible rod with free ends, initially inclined position shown in Fig. IO, is force and a torque applied simultaneously The time history of T(t) is shown in F(t) = T(t)/lO. The first two revolutions Fig, I I 5.3. A displacement arm
driven
placed in an subject to a at one end. Fig. 8, and are shown in
multi-component
robot
The robot arm in example I is now considered to be composed of two flexible components connected together by a hinge at mid-length. The two components were subjected to the prescribed rotation as shown in Fig. 12. The material properties are: EI = 100, EA = 10000, m = I, K = 0.867 for both links. The sequence of motion is shown in Fig. 13. Five elements are used in each links. 5.4. A jorce dricen multibody
s,vstem in free flight
A two body system consisting of two flexible links connected by a hinge is initially at an inclined position. The system is set into motion by applying a force and a torque at one end of the lower link, as shown in Fig. 14. The material properties are the same as in example 3. Subsequently the articulated beams undergo free flight that is shown in Fig. 15. 6. CONCLUDING
REMARKS
In the formulation used herein all Coriolis and centrifugal effects, as well as the other inertia effects due to rotation, are automatically accounted for. The simulation program can be used directly for the
Liu
dynamic analysis of a multi-link flexible robot arm in high speed motion. This approach can be carried over into the ful! three-dimensional case. This extension relies on a proper treatment of three-dimensional finite rotations in both the structural deformations of the beam and in the overall motions. All the simulation results reported here were calculated by a FORTRAN program named BODYFLEX written by the second author. This program has been tested on a few UNIX based systems (SUN-OS, NEWS-OS and ESIX). Graphics presentation was made through X-windows and a PostScript printer by IMSL Exponent Graphics system.
REFERENCES
1. J. C. Simo and L. Vu-Quot. beams under large overall 849-863 (1986). 2. K. M. Hsiao and J. Y. Jang, flexible mechanisms by Compul.
Mefh.
On the dynamics of flexible motions. J. appl. M&I. 53, Dynamic analysis of planar co-rotational formulation.
Appi. Mech. und Engng 87, IL14 (1991).
nonlin3. E. M. Bakr and A. A. Shabana, Geometrically ear analysis of multi-body systems. Compur. Struck. 23(6), 7399751 (1986). 4. J. Meek and H. Liu, Nonlinear dynamic analysis of flexible beams under large motions. In Compurutional Mechanics (Edited by Valhappan, Pulmano and TinLoi), Vol. 1, pp. 1755182. Balkema. Rotterdam (1993). modelling and simu5. J. Meek and H. Liu. Nonlinear lation of a flexible manipulator. Proc. lsr Pun-Pacific Conf. on Computational Engineering. Seoul. Korea (1993). Analytical Dvuzmics of Discrete 6. R. M. Rosenberg, Systems. Plenum Press, New York (1987). 1. J. Argyris and H. P. Mlejnek, D_wamics of structure. North-Holland, Amsterdam (1991).