Journal of Sound and Vibration (1995) 181(2), 261278

DYNAMICS OF FLEXIBLE BEAMS UNDERGOING

OVERALL MOTIONS
H. H. Y
Department of Mechanical Design and Production Engineering, Hanyang University,
Haengdang-Dong 17, Sungdong-Gu, Seoul, Korea, 133-791
R. R. R
Mechanical Dynamics, Inc., 2301 Commonwealth Boulevard, Ann Arbor, Michigan 48105,
U.S.A.

R. A. S
Department of Mechanical Engineering and Applied Mechanics, The University of Michigan
at Ann Arbor, Michigan 48109, U.S.A.

(Received 19 April 1993, and in final form 4 January 1994)

A modelling method for straight beams undergoing large overall motions as well as small
elastic deformations is presented in this paper. Different from the classical linear Cartesian
modelling method which employs three Cartesian deformation variables, the present
modelling method uses a non-Cartesian variable along with two Cartesian variables to
describe the elastic deformation. A quadratic form of the strain energy expressed with the
hybrid set of deformation variables is used to obtain the generalized active forces and a
geometric constraint equation relating the non-Cartesian variable and Cartesian variables
is used to obtain the generalized inertia forces in the equations of motion. The present
modelling method not only provides accurate simulation results but also clarifies the limit
of validity of the classical linear Cartesian modelling method.

1. INTRODUCTION
A dynamic modelling method is presented for straight beams undergoing large overall
motions as well as small strain elastic deformations. Flexible structures are often idealized
as beams and their dynamic characteristics are represented by those of beams. Large
overall motions consist of translational and rotational motions. They are frequently called
rigid body motions. There are several engineering examples which have beam-like shapes
and undergo large overall motions. Satellite antennas and helicopter blades are such
examples. Engineering examples of flexible structures will appear more frequently, since
the general structural design trend is toward lighter systems with increasingly faster
dynamic response and minimal power requirements.
Structural flexibility can cause problems for the standard operational motions of a
system. For instance, in 1958 the Explorer I spacecraft experienced a radical instability,
which was blamed on energy dissipation induced by the vibration of long antennas of the
satellite [1]. To avoid this kind of problem, the dynamic characteristics of structure should
be identified accurately so that one can properly design or control the system.
The first modelling approach for beams, which hereafter will be referred to as the C.L.C.
(classical linear Cartesian) approach, was introduced [2, 3] in the 1970s when the speed of
261
0022460X/95/120261 + 18 $08.00/0 7 1995 Academic Press Limited
262 . . .
computer and numerical methods was progressing rapidly. This approach is based on the
classical linear elastic modelling, where geometric as well as material linearity is assumed.
It has merits such as simplicity of formulation, ease of implementation in finite element
methods, and available co-ordinate reduction techniques [4], which are often critically
important for dynamic analyses of structures. Because of these merits, the C.L.C. approach
has been widely used by many engineers. However, this approach displays a critical flaw
when the structures undergo large overall rotational motions. In this case, the approach
produces erroneous results in numerical simulations.
To resolve the problem of the C.L.C. approach, several non-linear methods [58] have
been introduced. Since the problem is believed to originate from the assumption of
geometric linearity, those methods are based on geometric non-linear relations between
strains and displacements. With the non-linear methods, the accuracy problem was shown
to be resolved. However, other problems were raised due to the non-linearity of the
modelling methods. Actually, these modelling methods lack the merits of the C.L.C.
approach, including the capability of co-ordinate reduction. Without this capability, a
large amount of computational effort is needed for dynamic analyses. Furthermore, a
question still remains for the C.L.C. approach: Exactly under what conditions does the
approach fail? If one knows the answer, one does not have to rely all the time on the
inefficient non-linear approach.
The purposes of this paper are to overcome the disadvantage of the non-linear approach
and to identify explicitly the limit of validity of the C.L.C. approach. The present work
is based on the work in reference [9]. Compared to the work in reference [9], the present
work simplifies the procedure of formulating the derivation of the equations of motion,
clarifies its relation to non-linear modelling, provides a simple explanation of the modelling
effectiveness, and presents more numerical examples, including a three-dimensional one.
In the following sections, the procedure for deriving the equations of motion is
described. To focus on the effect of large overall motions on the elastic behavior of flexible
beams, large overall motions are prescribed. Numerical examples show the importance of
the effect. In the penultimate section, a simple discrete example is introduced. The discrete
system is designed to be analogous to a rotating beam so that critical physical phenomena
involved in the beam can be easily understood through the system.

2. SYSTEM AND VARIABLE DESCRIPTION

The following assumptions are made. The beam has homogeneous and isotropic
material properties. The elastic and centroidal axes in the cross-section of the beam
coincide so that no eccentricity effect needs to be considered. The beam has a slender shape
so that shear and rotary inertia effects are neglected. Large overall motions are prescribed
as functions of time. No external forces are given in the system of concern. These
assumptions are made to simplify the formulation and to focus on the critical effects
induced by large overall motions which make differences between the C.L.C. and the
present modelling methods. However, they can be removed without having many
difficulties if necessary. Such cases have been well described in reference [10].
In Figure 1 are shown the configurations of a straight beam before and after
deformation. The deformation can be represented by three scalar variables in 3-D space.
Conventionally, Cartesian variables, as shown in the figure, are widely used. They are used
in the C.L.C. modelling method. In the present modelling method, a non-Cartesian
variable along with two Cartesian variables is used instead of three Cartesian variables.
s, denoting the arc length stretch, is used instead of u1 , which denotes the Cartesian
distance measure of a generic point in the axial direction of the undeformed configuration.
 263

Figure 1. The configuration of a beam before and after deformation.

u2 and u3 , which denote the lateral displacements, are also used in the present modelling
method.
There is a geometric relation between the arc length stretch s and the Cartesian variables.
This relation is later used to derive generalized inertia forces in the equations of motion.
In reference [9], the relation is given as

g $ 0 1 0 1%
1u2 1u3
x + u1 2 2 1/2

x+s= 1+ + ds. (1)

0
1s 1s
In words, this equation expresses the fact that the arc length distance from the reference
point O of the beam B to some generic particle P of the beam can be expressed as either
(1) the sum of the distance x (from the reference point O to the undeformed position of
a generic point P0 ) and the neutral axis stretch s, or (2) the integral, from zero to x + u1
of a differential arc length of the beam. This formulation has analytical and numerical
difficulties associated with it: namely, the determination of u1 . The explicit dependence of
the integral on u1 complicates the procedure of deriving the equations of motion and makes
it difficult to obtain the numerical value of u1 . In the present work, a different strategy
is adopted to attack the difficulties. The drawbacks of the form presented in equation (1)
can be overcome by changing the domain of integration to include only the undeformed
length instead of the deformed length and altering the integrand accordingly. Thus, as an
equivalent to equation (1), one can write

g $0 1 0 1 0 1%
1u1 1u2 1u3
x 2 2 2 1/2

x+s= 1+ + + ds (2)
0
1s 1s 1s
and, by using a binomial expansion of the integrand of equation (2), this can be shown
to give

g $0 1 0 1 %
1u2 1u3
x 2 2

s = u1 + 12 + ds + (higher order terms). (3)

0
1s 1s
One of the eventual goals of the present modelling method is the development of equations
of motion which are linear in the deformation variables s, u2 and u3 . Therefore, it seems
sensible to search for an approximate form of s which will yield only those linear and
quadratic terms appearing explicitly in equation (3). Such an approximate form is

g $0 1 0 1 %
1u2 1u3
x 2 2

s = u1 + 12 + ds. (4)
0
1s 1s
264 . . .
Differentiation of this expression for s with respect to x yields

0 1 $0 1 0 1 %
1s 1u1 1u2 1u3
2 2

= + 12 + . (5)
1x 1x 1x 1x

This equation shows that the differentiation of the approximate s is equivalent to the well
known von Karman strain measure for beams. Actually, without any approximation,
1s/1x represents the exact strain for the stretching of a beam. The use of equation (4) not
only simplifies the derivation of the equations of motion but also provides a simple way
of determining u1 numerically when the variables s, u2 and u3 are known.

3. APPROXIMATION AND STRAIN ENERGY

In the previous section, it was mentioned that the deformation variables, s, u2 and u3
are used in the present modelling method. This implies that these continuous variables
should be approximated by using spatial functions and corresponding coordinates to
derive the ordinary differential equations of motion. In the present work, the Rayleigh
Ritz method is used to approximate the variables as follows:
m m m
s(x, t) = s f1j (x)qj (t), u2 (x, t) = s f2j (x)qj (t), u3 (x, t) = s f3j (x)qj (t). (68)
j=1 j=1 j=1

Here f1j , f2j and f3j are the spatial functions. These functions can be obtained from the
vibration analysis of the beam undergoing no rigid body motion. f1j is the modal function
of the longitudinal vibration, and f2j and f3j are modal functions of the bending vibrations
in the two lateral directions. The qj s are generalized co-ordinates and m is the total number
of modal co-ordinates. The same modal functions are often used for u1 , u2 , and u3 in the
C.L.C. modelling method. For the convenience of the formalism, s, u2 and u3 explicitly have
the same number of co-ordinates m. However, they are not actually coupled. For instance,
f1j is not zero only if j E m1 , f2j is not zero only if m1 Q j E m1 + m2 and f3j is not zero only
if m1 + m2 Q j E m1 + m2 + m3 . In other words, m1 , m2 and m3 denote the actual numbers of
generalized co-ordinates for s, u2 and u3 , respectively. m is the total sum of m1 , m2 and m3 .
In the C.L.C. modelling method, u1 (instead of s in equation (6)) is approximated and
the strain energy expression is

g 0 1 g 0 1 g 0 1
1u1 1 2u2 1 2u3
L 2 L 2 L 2

U = 12 EA dx + 12 EIzz dx + 12 EIyy dx, (9)

0
1x 0
1x 2 0
1x 2

where E denotes Youngs modulus, A is the cross-sectional area of the beam, Izz and Iyy
are the second area moments of the cross-section, and L is the undeformed length of the
beam. For the present modelling method, however, the following strain energy expression
is used:

g 0 1 g 0 1 g 0 1
1s 1 2u2 1 2u3
L 2 L 2 L 2

U = 12 EA dx + 12 EIzz dx + 12 EIyy dx. (10)

0
1x 0
1x 2 0
1x 2

Equations (9) and (10) differ only in their first terms, which represent the stretching energy
of the beam. Since s represents the exact stretch of the beam, the first term of equation
(10) represents the exact stretching energy of the beam. On the other hand, since 1u1 /1x
of equation (9) is the linear approximation of 1s/1x, the first term of equation (9) does
not represent the exact stretching energy. With the three Cartesian variables, the exact
 265
stretching energy can be expressed in a non-quadratic form which results in non-linear
active forces in the equations of motion. Actually, this is how the non-linear modelling
methods in references [6] and [7] express the strain energy forms. For example, if the right
side of equation (5), which is the approximation of 1s/1x, is substituted into equation (10),
one can obtain the strain energy expression based on the von Karman strain measure,
which results in non-linear generalized active force terms. On the other hand, with the
quadratic form of strain energy given in equation (10), the generalized active force terms
become linear and the exact stretching energy can be considered. However, the use of s
complicates the formulation of generalized inertia forces in the equations of motion and
results in non-linear generalized inertia force terms. The non-linear inertia force terms are
linearized to obtain the final equations of motion for the present modelling method.

4. EQUATIONS OF MOTION
In the present section, Kanes method is employed to derive the equations of motion.
This method provides a direct and systematic way of deriving ordinary differential
equations of motion for continua such as beams, plates and shells. With the assumptions
given in section 2, the equations of motion can be obtained for beams as follows:

g 0 1
1v P P 1U
L

r a dx + = 0, i = 1, . . . , m. (11)
0
1qi 1qi

r is the mass per unit length of the beam, the qi s are the time derivatives of the generalized
co-ordinates, and v P and a P are the inertial velocity and acceleration of the center point
P of a differential element of the beam (see Figure 2). The acceleration can be obtained
simply by differentiating the expression for the inertial velocity,
v P = v O + v A (r + u ) + Av P, (12)
where O is a reference point identifying a point fixed in a floating frame A in the system,
v O is the inertial velocity of point O, v A is the inertial angular velocity of the floating
reference frame, r is a position vector from O to the location of P in the undeformed body,
u is the displacement vector of point P, and Av P is the relative velocity of P in frame A
obtained by taking the time derivative of u in frame A. In component notation, the various
terms can be written as
r = xa 1 , u = u1 a 1 + u2 a 2 + u3 a 3 (13, 14)
v O = v1 a 1 + v2 a 2 + v3 a 3 , v A = v1 a 1 + v2 a 2 + v3 a 3 . (15, 16)

Figure 2. The kinematics of a beam undergoing large overall motion.

266 . . .
By using equations (4), (6)(8) and (13)(16), the partial derivative of the velocity of
P with respect to the generalized speed qi can be obtained as

$ 0g 1 0g 1%
1v P m m
x x

= f1i s f2i,h f2j,h dh qj s f3i,h f3j,h dh qj a 1 + f2i a 2 + f3i a 3 , (17)

1qi j=1 0 j=1 0

where h denotes partial differentiation with respect to the dummy spatial variable h.
Non-linear equations of motion can be derived from equation (11) by using equation (17),
along with the expression for the system strain energy and the acceleration of P. The final
procedure in the present modelling method is to linearize the equations with respect to the
generalized co-ordinates for the elastic deformations s, u2 and u3 . The results are described
in equations (18)(20), which will be presented in what follows. The whole procedure is
straightforward except for the use of the integration by parts technique employed to obtain
some motion-induced stiffness terms in the equations of motion. Those terms appear in
the second and third lines of equations (19) and (20). If the terms are deleted, the resulting
equations become equivalent to those of the C.L.C. modelling method. Equations
(18)(20) are as follows:

$0g 1 0g 1 0g 1%
m L L L

s rf1i f1j dx qj (v22 + v32) rf1i f1j dx qj + EAf1i, x f1j, x dx qj

j=1 0 0 0

$ 0g 1 0g 1%
m L L

s 2v3 rf1i f2j dx qj (v1 v2 v 3 ) rf1i f2j dx qj

j=1 0 0

$ 0g 1 0g 1%
m L L

+s 2v2 rf1i f3j dx qj + (v1 v3 + v 2 ) rf1i f3j dx qj

j=1 0 0

g g
L L

= (v22 + v32) rxf1i dx (v1 + v2 v3 v3 v2 ) rf1i dx; (18)

0 0

$0g 1 0g 1 0g 1%
m L L L

s rf2i f2j dx qj (v12 + v32) rf2i f2j dx qj + EIzz f2i, xx f2j, xx dx qj

j=1 0 0 0

$ 0g 1%
m L

s (v1 + v2 v3 v3 v2 ) r(L x)f2i, x f2j, x dx qj

j=1 0

$ 0g 1%
m
r 2
L

+s (v22 + v32) (L x 2)f2i, x f2j, x dx qj

j=1 0
2

$ 0g 1 0g 1%
m L L

s 2v1 rf2i f3j dx qj (v2 v3 v 1 ) rf2i f3j dx qj

j=1 0 0

$ 0g 1 0g 1%
m L L

+s 2v3 rf2i f1j dx qj + (v1 v2 + v 3 ) rf2i f1j dx qj

j=1 0 0

g g
L L

= (v1 v2 + v 3 ) rxf2i dx (v2 + v3 v1 v1 v3 ) rf2i dx; (19)

0 0
 267

$0g 1 0g 1 0g 1%
m L L L

s rf3i f3j dx qj (v12 + v22) rf3i f3j dx qj + EIyy f3i, xx f3j, xx dx qj

j=1 0 0 0

$ 0g 1%
m L

s (v1 + v2 v3 v3 v2 ) r(L x)f3i, x f3j, x dx qj

j=1 0

$ 0g 1%
m
r 2
L

+s (v22 + v32) (L x 2)f3i, x f3j, x dx qj

j=1 0
2

$ 0g 1 0g 1%
m L L

s 2v2 rf3i f1j dx qj (v1 v3 v 2 ) rf3i f1j dx qj

j=1 0 0

$ 0g 1 0g 1%
m L L

+s 2v1 rf3i f2j dx qj + (v2 v3 + v 1 ) rf3i f2j dx qj

j=1 0 0

g g
L L

= (v1 v3 v 2 ) rxf3i dx (v3 + v1 v2 v2 v1 ) rf3i dx. (20)

0 0

One of the merits of the present modelling method is that it gives a clear indication of
the factors involved in motion-induced stiffness variations. For example, as one may see
from the second and third terms in square brackets of equations (19) and (20),
(v1 + v2 v3 v3 v2 ) and (v22 + v32) are the overall motions involved in the stiffness variations.
Such overall motion terms do not appear in any non-linear modelling methods (see
references [5, 6]). Instead, such effects are compensated by non-linear elastic force terms
in the non-linear modelling methods. However, in the C.L.C. modelling method, neither
the overall motion terms nor non-linear elastic force terms appear in the equations of
motion. As such overall motion terms become significant for the system, those terms
should be included to obtain accurate system responses. However, if the overall motion
terms are not significant to the system, the C.L.C. modelling method can be used without
any trouble. In other words, by inspecting the explicit factors discussed above, one may
easily determine whether or not the C.L.C. modelling method can be safely used.

5. NUMERICAL EXAMPLES
A cantilever beam attached to a rigid base is shown in Figure 3. The base performs a
planar rotational motion around the vertical axis which is given as

6 7
Vs [6t 5 15t 4 + 10t 3] if 0 E t E 1
V= , (21)
Vs if t e 1
where
t ,t/Ts . (22)

Figure 3. Abrupt planar spin-up motion of a cantilever beam. Mass per unit length r = 550 kg/m; Youngs
modulus E = 689 E 9 N/m2, length L = 305 m, cross-sectional area A = 929 E 2 m2, area moment of inertia
I = 719 E 4 m2, time constant Ts = 10 s, steady state angular speed Vs = 031 rad/s.
268 . . .

Figure 4. Simulation results for the abrupt spin-up example. , Present; , C.L.C.

Here Vs and Ts are the steady state angular speed and the time to reach the angular speed,
respectively. This example was introduced and solved by the non-linear modelling method
in reference [5]. All the numerical data are converted to SI units in this paper. The
simulation results are shown in Figure 4. The axial and transverse tip deflections (denoted
as u1 and u2 ) were obtained by using the present modelling method (solid lines) and the
C.L.C. modelling method (dotted line). It is clearly shown that the two methods produce
quite different numerical results. However, it is found that the results obtained by the
present modelling method are almost identical to those obtained by the non-linear
modelling method in reference [5]. The foreshortening effect (reduction of u1 due to the
change of lateral deflections u2 or u3 ) can be observed from the result obtained by the
present modelling method. As the lateral deflection becomes large, large negative axial
deflection results (as shown in Figure 4). As mentioned in section 2, u1 can be determined
numerically by using equation (4) once s, u2 and u3 have been obtained by numerical
integration. It should be noted that the maximum lateral deflection in the plot is quite
large. Thus, the present modelling method is as effective as the non-linear modelling
method for solving large deflection problems of beams.
Another similar numerical example is shown in Figure 5. For this example, the
prescribed base motion is smooth and is given as

6 7
Vs [t (1/2p) sin (2pt)] if 0 E t E 1
V= , (23)
Vs if t e
where t is defined same as in equation (22). This motion smoothly increases the angular
speed until it reaches the steady state and the steady state angular speed is sustained. It
is so smooth and slow that the lateral oscillation after reaching the steady state remains
quite small. In Figure 6 are shown the simulation results (lateral deflections) by the present

Figure 5. Smooth planar spin-up motion of a cantilever beam. Mass per unit length P = 12 kg/m, Youngs
modulus E = 70 E 10 N/m2, length L = 10 m, cross-sectional area A = 40 E 4 m2, area moment of inertia
I = 20 E 7 m4, time constant Ts = 120 s, steady state angular speed Vs = 60 rad/s.
 269

Figure 6. Simulation results for the smooth spin-up example. , Present; , C.L.C.

modelling method (solid line) and the C.L.C. modelling method (dotted line). The figure
shows that even for the relatively small deflection case the C.L.C. modelling method fails
to provide reasonable results. It is often believed that the C.L.C. modelling method can
be safely used unless the displacement becomes relatively large. This example shows that
such an idea is not always correct. The effective range of the C.L.C. modelling method
can be identified only by inspecting the factors involved in the motion induced stiffness
terms (second and third terms in square brackets of equations (19) and (20)). The terms
should not be ignored if they are significant enough in comparison with the elastic stiffness
terms in the equations of motion, whereas the C.L.C. modelling method can be as effective
as the present modelling method if they are relatively small.
For the purpose of exhibiting the three-dimensional capability of the formulation
derived in this paper, a problem of a rotating cantilever beam attached to a rigid base with
a taper angle a has been solved. The prescribed function in equation (23) is used for the
numerical analysis. With the numerical data given in Figure 7, simulation results were
obtained and are plotted in Figure 8. The results show the difference between the results
of the present modelling method (solid line) and those of the C.L.C. modelling method
(dotted line). The steady state lateral deflection u2 does not converge to zero since the
centrifugal inertia force acting perpendicular to the rotating axis deforms the beam
outwards. The magnitudes and frequencies of the steady state oscillation of u3 obtained
by the two modelling methods are also shown to be quite different from each other.
The steady state deformed shape of the beam (in the plane of a 1 and a 2 ) is shown in
Figure 9.

Figure 7. Three-dimensional motion example of a cantilever beam. Mass per unit length P = 12 kg/m,
Youngs modulus E = 7 E 10 N/m2, length L = 10 m, cross-sectional area A = 4 E 4 m2, area moments of
inertia Iyy = 20 E 7 m4 and Izz = 40 E 7 m4, time constant Ts = 15 s, steady state angular speed Vs = 3 rad/s.
270 . . .

Figure 8. Simulation results for the three dimensional motion example. (a) Lateral deflection u2 ; (b) lateral
deflection u3 . , Present; , C.L.C.

6. ANALOGY STUDY
With the results obtained in the previous section, the following questions may be raised.
Why does the C.L.C. modelling method lack proper stiffness variation terms? How can
the problem of the C.L.C. modelling method be resolved? How is the stiffness variation
effect captured by the present modelling method? Does the choice of the deformation
variable play a role in obtaining a proper modelling method? What is the key ingredient
of the non-linear modelling method? These questions need to be addressed to provide
better insight into dynamic modelling of flexible structures undergoing large overall
motions.

Figure 9. The steady state deformed shape of the 3-D example.

 271

Figure 10. A discrete system analogous to a rotating cantilever beam. Link A is rigid; inner wall of B is
frictionless.

Due to the complexities involved in structural theories, it may not be easy to obtain
proper answers to the questions raised in the previous paragraph. Therefore, a simple
discrete system which is analogous to a flexible structure undergoing large overall motion
was devised. Due to the analogy, the discrete system exhibits characteristic behavior similar
to that of flexible structures. Since the discrete system possesses only two degrees of
freedom, the exact equations of motion can be easily derived and solved numerically.
Furthermore, the mathematical modelling analogy (between the discrete system and
continuous flexible structures) enables one to obtain the answers to the previously raised
questions. The mathematical models developed for the discrete system have their
analogous counterparts in the continuous modelling methods. Therefore, the conclusions
drawn from the study of discrete system can shed light on dynamic modelling strategy for
flexible structures.
The discrete mechanical system depicted in Figure 10 is representative of the class of
problems of interest, in that it involves large overall motions as well as small elastic
deformations. The system consists of the following elements: a rigid link A of the length
d, attached to ground through a revolute joint at point C; a smooth massless tube B,
connected to A via a revolute joint at point Q; a torsional spring of modulus kr acting
between A and B; a particle P, of mass m, which slides freely along the smooth inner wall
of tube B; and a translational spring, of modulus kT , which connects P to point Q. When
the translational spring is undeformed, the distance between P and Q is denoted by the
symbol L. The rotational motion of A, defined in terms of the magnitude V of its inertial
angular velocity vector, is prescribed as a function of time. Thus, the system possesses two
degrees of freedom and the equations of motion can be derived in terms of any two of
the four variables u, s, x and y, introduced in the figure. The symbol u refers to the radian
measure between the line of action of the translational spring force and the line extending
through points C and Q. The stretch in the translational spring is represented by the
symbol s, and the variables x and y refer to the axial and transverse displacements,
respectively, of P relative to the line CQ. The rotational and translational springs are
undeformed when u = 0 and s = 0, respectively. This model is similar to a rotating
cantilever which exhibits both extensional and bending flexibility. The extensional and
bending stiffnesses are represented by the translational and the rotational springs,
respectively.
The main objective of the present section is to illustrate how different choices of
deformation co-ordinates coupled with linearization result in different mathematical
models which lead to radically different predictions of system behavior. This can provide
272 . . .
T 1
Five mathematical models for the discrete system
NC set: fully non-linear in x and y:
mx mV y 2mVy mV 2(d + L + x) + 1U1 /1x + 1U2 /1x = 0,
my + mV x + 2mVx mV 2y + 1U1 /1y + 1U2 /1y = 0,
U1 = 12 kT {[(L + x)2 + y 2]1/2 L}2, U2 = 12 kr [tan1 {y/(L + x)}]2.
NP set: fully non-linear in s and u:
ms + mV d sin u m(L + s)(u + V)2 mdV 2 cos u + 1U1 /1s = 0,
(L + s)[m(L + s)u + mV (L + s + d cos u) + 2ms (u + V) + mdV 2 sin u] + 1U2 /1u = 0,
U1 = 12 kT s 2, U2 = 12 kr u 2.
LC set: linear in x and y:
mx mV y 2mVy mV 2(d + L + x) + kT x = 0,
my + mV x + 2mVx + (mV 2 + kr /L 2)y = 0,
U1 = 12 kT x 2, U2 = 12 kr (y/L)2.
LP set: linear in s and u:
ms + mV du 2mLVu mV 2(d + L + s) + kT s = 0,
mL u + mV L(L + d) + mV (2L + d)s + 2mLVs + (mdLV 2 + kr )u = 0,
2

U1 = 12 kT s 2, U2 = 12 kr u 2.
NS set: non-linear in x and y, up to second degree stretch related to kT :
mx mV y 2mVy mV 2(d + L + x) + kT x + 12 kT (y 2/L) = 0,
my + mV x + 2mVx + (mV 2 + kr /L 2)y + kT (xy/L) + kT (y 3/2L 2) = 0,
U1 = 12 kT (x + y 2/2L)2, U2 = 12 kr (y/L)2.

important insight into any modelling methods intended for the treatment of
motion-induced structural stiffness variations. The two sets of co-ordinates to be studied
in this example are the Cartesian co-ordinate set x and y and the non-Cartesian (polar)
co-ordinate set s and u.
In Table 1 are listed five sets of dynamic equations relevant to the system represented
in Figure 10. The table also contains five sets of strain energy expressions used to obtain
the five dynamical equation sets. U1 and U2 represent the strain energies stored in the
translational and the torsional springs, respectively.
The first two sets of equations, labelled NC (non-linear Cartesian) and NP (non-linear
polar), represent the fully non-linear equations of motion expressed in terms of (x, y) and
(s, u), respectively. One can easily demonstrate their mathematical equivalence. Thus,
when numerical simulations are performed, the two models always produce equivalent
results. Actually, any fully non-linear exact model can produce equivalent results.
Now, if straightforward mathematical linearization of set NC in x and y and set NP
in s and u is performed (all terms of degree two and higher in x, x , y, y , s, s , u and u
are dropped), one obtains the third and fourth equation sets, respectively. These two sets
of equations, labelled LC (linear Cartesian) and LP (linear polar) are similar in form, with
one crucial difference. The second equation in set LP contains a stiffness term which
increases in magnitude with increasing V, thus exhibiting the well known centrifugal
stiffening effect. In contrast to this, the second equation of set LC involves a stiffness term
 273
which decreases in magnitude with increasing V, thereby effectively softening the rotational
spring. This is in direct conflict with experimentally observed behaviors. As will be shown
later in this section, the solution of set LP matches closely that of its non-linear counterpart
NP for a wide range of V values; such is not the case, however, for LC and NC. The
important fact is that the two linearization procedures based on different co-ordinate sets
results in two qualitatively different mathematical models. It is possible to obtain such a
result since the exact models (NC or NP) are linearized only in co-ordinates, but not in
V, which is not generally small.
Since the development of equation set LC involves no linearization in inertia force
contributions, the problem of the LC set should be caused by dropping non-linear terms
of the active force in the equations of motion. In Figure 11(a) are shown the centrifugal
force F c and the translational spring force F k acting on the particle P. Also shown in the
axial and transverse directional components of the two forces, F cx , F cy F kx and F ky . Note
that the transverse direction component of the centrifugal force has a destabilizing effect
on the particle P (destabilizing simply means that the force component pushes the
particle away from the undeformed position: i.e., x = y = s = u = 0). This leads to the

Figure 11. The resolution of centrifugal and translational spring forces in Cartesian and polar components.
(a) Resolution of forces in Cartesian (axial and transverse) components (e x and e y are unit vectors fixed in A);
(b) resolution of forces in polar (radial and tangential) components (e r and e u are unit vectors fixed in B).
274 . . .
softening effect on the rotational spring in the LC models. Thus, the destabilizing
motion-induced stiffness term mV 2y appears in the LC set. However, the actual problem
of the LC model lies in the linearization of the translational spring force. Actually, the
same destabilizing effect is included in the NC set. However, in the NC set, the destabilizing
effect is compensated by non-linear terms in the translational spring force. The
translational spring force is dependent on the extension in the spring. Hence, in the
mathematical expression for the spring force, system co-ordinates appear at least to the
first degree (linearly) or higher. When this spring force is then resolved into two
components, parallel to e x and e y (see the sketch) respectively, the mathematical expression
for the e y component of the force will be of the second degree or higher in the system
co-ordinates due to the trigonometric sine of the angle between F k and F kx . F ky is, therefore,
not included in the LC model. The truncation of F ky would not be significant if F ky were
small compared to other linear terms in the equations of motion. Unfortunately, F ky is
significant, since it involves the translational spring coefficient kT , which is very large
compared to any other parameters of the system. For instance, kT is much larger than kr
since the discrete system is devised to imitate a flexible structure in which the stretching
stiffness is much larger than the bending stiffness.
Unlike the equation set LC, the equation set LP involves no linearization of the active
force contributions. It involves linearization only in the inertia force contributions.
Therefore the LP set does not have the problem linearizing the translational spring force
since it expresses the force in an exact linear form. The linearization of the inertia force,
however, does not cause any serious problems in the LP set. Since the inertia force has
an expression of zero degree and higher, the LP set contains all the zero and first degree
terms of inertia force with the linearization. Furthermore, the magnitude of the parameter
(m) involved in the inertia force is not very large compared to other parameters of system.
Therefore, the truncation of the inertia force terms which are second degree and higher
does not affect the fidelity of the LP set (to the NP set). In Figure 11(b) are shown the
centrifugal force F c and the translational spring force F k acting on the particle P. Also
shown are the radial and tangential direction components, F cs , F cu , F ks and F ku . Unlike the
transverse direction component of the centrifugal force in Figure 11(a), the tangential
direction component in Figure 11(b) has a stabilizing effect. The stabilizing term mdLV 2y
appears in the LP set. Therefore, the combination of linearized forces (by F cu and F ku ) has
a stabilizing effect. This leads to the fundamental difference (stiffness variations due to
large overall motion) between the LC and LP models.
Since the failure of the LC model results from the truncation of non-linear terms in the
active force contribution, one can overcome the difficulty by keeping all non-linear active
terms (as in the NC set). However, is it necessary to keep all the non-linear terms to develop
an accurate model? To find the answer to this question, the NS set was generated as
follows. The stretch length of the translational spring is approximated up to second degree
in terms of x and y and retained in the strain energy expression associated with the
translational spring. The effectiveness of this set, in comparison with NC set, was then
evaluated by using numerical simulations in order to determine whether such a model does
indeed capture the motion-induced stiffness variation effect.
It is worthwhile to note the analogies between the discrete models and the structural
modelling method discussed in the previous sections. The fully non-linear NC model is
analogous to the non-linear modelling method (introduced in reference [5]), which is based
on fully non-linear exact strain measures for flexible structures. The LC model is analogous
to the C.L.C. modelling method, since they both employ Cartesian deformation measures
and linear quadratic strain energy expressions. The LP model is analogous to the modelling
method developed in this paper, since they are both based on non-Cartesian deformation
 275
T 2
Data of the discrete system used for numerical simulations
Notation Numerical data Description
m 1 kg Mass of particle P
L 1m Unstretched length of translational spring
d 1m Distance from C to Q
kT 600 000 N/m Translational spring stiffness
kr 300 Nm/rad Torsional spring stiffness
Vs 10 rad/s Steady state angular speed of spin-up case 1
Vs 20 rad/s Steady state angular speed of spin-up case 2
Ts 1s Steady-state time constant for spin-up cases
x(0) 0m Initial value of x
y(0) 0m Initial value of y
x (0) 0 m/s Initial value of x
y (0) 0 m/s Initial value of y
t 02 s Simulation time interval

measures with which the exact strain energy (for stretching) can be expressed in quadratic
forms. Lastly, the NS model is analogous to the non-linear modelling method based on
von Karman strain measures (see, for instance, reference [11]) since they are both based
on the second degree strain measures expressed in Cartesian variables.
Simulations were performed by using the parameters listed in Table 2. For the large
overall motion (represented by V) of the system, a spin-up function (given in equation (23))

Figure 12. A comparison of the numerical results obtained by using the NC () and LC () models.
(a) Angular speed Vs = 10 rad/s; (b) angular speed Vs = 20 rad/s.
276 . . .
was used. Since the spin-up function varies with time, it is quite effective in representing
general cases of overall motions.
In Figures 12 and 13 is provided a qualitative idea of the fidelity of the linear equation
sets LC and LP relative to their non-linear counterparts NC and NP. The lateral
displacements y are plotted in the figures. In Figure 12 it is illustrated how poorly the
solution of the LC set approximates that of NC. For the slower spin rate, Vs = 10 rad/s,
the errors are substantial. With the faster spin rate, Vs = 20 rad/s, the linear equations
incorrectly predict an unbounded response. The former response (with Vs = 10 rad/s) is
actually more disturbing than the latter (with Vs = 20 rad/s), because the former response
might appear reasonable in the absence of the correct response while the latter response
is obviously in conflict with observed physical behavior. Figure 13, on the other hand,
provides a glimpse of the impressive accuracy of the linear equation LP set in tracking
the true solution even at high spin rates. For convenience of comparison, output results
of all models are converted to the transverse displacement y. For instance, once s and u
had been obtained numerically in the LP set, y was calculated by using the trigonometric
relation.
Additional simulations using the NS set were performed. The results are depicted in
Figure 14. In each plot, the solid line represents the solution obtained by using the NC
equation set and the dashed curve represents the solution obtained by using the NS
equation set. As shown in the plot, the two curves are close to each other. Note that,
however, the numerical results given by the NS equation set are no better than those given
by the LP equation set, despite the additional effort of computation for high degree terms.
If non-linear terms based on third or fourth degree stretch (or strain) are kept in the

Figure 13. A comparison of the numerical results obtained by using the NP () and LP () models. (a)
Angular speed Vs = 10 rad/s; (b) angular speed Vs = 20 rad/s.
 277

Figure 14. A comparison of the numerical results obtained by using the NC () and NS () models.
(a) Angular speed Vs = 10 rad/s; (b) angular speed Vs = 20 rad/s.

equations of motion, the accuracy of the mathematical model will be improved. However,
such models are not investigated here, since the counterparts (structural modelling) are not
as meaningful as the von Karman modelling (counterpart for the NS set) is.
Through the introduction and analysis of a simple discrete system, a number of
important modelling issues such as motion-induced stiffness variations, their contribution
to linear and non-linear equations of motion, their dependence on co-ordinate choice,
and the influence of linearization procedures in accounting for such effects have been
investigated. If one chooses a special set of deformation variables such that the strain
energy of the structure is accurately characterized in a quadratic form in these variables,
then it is possible to develop differential equations of motion which are totally linear in
the deformation variables and yet still exhibit proper stiffness variations. This is clearly
shown by the performance of the LP set. This is the approach taken in the modelling
method developed in this paper, wherein the strain energy of the 3-D beam is expressed
in a quadratic form by using the stretch variable of the beam. If one chooses the Cartesian
set of deformation variables, the stiffness variation due to large overall motion cannot be
properly predicted by a linear model. This is shown by the performance of the LC set.
The counterpart (the C.L.C. modelling method) has shown spurious performances. If
proper motion-induced stiffness variation is to be accommodated, at least up to second
degree terms in the stretch strain should be retained. This is shown by the performance
of the NS set. This is consistent with the results of the non-linear modelling method [11]
based on von Karman strain measures.
278 . . .
7. CONCLUSIONS
The present paper mainly serves to propose an accurate linear modelling method for
flexible beams undergoing large overall motion and concomitant small linear elastic
deformation. The crucial part of the present modelling method is the introduction of a
stretch deformation variable by which the strain energy function can be expressed in a
quadratic form (effectively resulting in linear generalized active forces in the equations of
motion). In this modelling method, motion-induced stiffness variations are captured
accurately and expressed in analytical forms in which large overall motion component
functions are explicitly shown. This enables one to predict easily the relation between the
amount of stiffness variations and large overall motions. By observing the stiffness
variation terms shown in the equations of motion (these terms are absent in the equations
of motion given by the C.L.C. modelling method), one may determine the range of validity
of the C.L.C. modelling method. In other words, if the stiffness variation terms are small
compared to the rest of the terms in the equations of motion, the C.L.C. modelling method
may provide reasonable simulation results.

ACKNOWLEDGMENT
This work has been partially supported by the Korea Ministry of Education through
Mechanical Engineering Research Fund (ME93-D-12), for which the authors are grateful.

REFERENCES
1. P. W. L and H. K. B 1971 Astronautics and Aeronautics 9(5), 6471. Attitude control
of non-rigid spacecraft.
2. J. Y. L. H 1977 Journal of Spacecraft and Rockets 14, 102110. Direct path method for flexible
multibody spacecraft dynamics.
3. C. S. B, A. D. D, A. C. P and H. P. F 1978 NASA Technical Paper 1219
1 & 2. A digital computer program for the dynamic interaction simulation of controls and
structure (DISCOS).
4. W. H, J. C and G. H 1971 Computers and Structures 1, 535563. Dynamic
analysis of large structures by modal synthesis techniques.
5. E. R. C and S. W. L 1986 Computers and Structures 23, 819829. Nonlinear finite
element modeling of the dynamics of unrestrained flexible structures.
6. J. C. S and L. V-Q 1986 Journal of Applied Mechanics 53, 849863. On the dynamics
of flexible beams under large overall motionsthe plane case: part I and part II.
7. T. B and B. H 1973 International Journal of Numerical Methods in Engineering
26, 255271. Nonlinear transient finite element analysis with convected coordinates.
8. W. J. H, R. R. R and R. A. S 1992 33rd Structures, Structural Dynamics, and
Material Conference, Dallas, Texas, AIAA 92-1162. A new flexible body dynamic formulation
for beam structures undergoing large overall motion.
9. T. R. K, R. R. R and A. K. B 1987 Journal of Guidance, Control, and Dynamics
10(2), 139151. Dynamics of a cantilever beam attached to a moving base.
10. R. R. R 1986 Ph.D. Dissertation, Stanford University. Flexibility modeling method in
multibody dynamics.
11. H. H. Y 1989 Ph.D. Dissertation, The University of Michigan at Ann Arbor. Dynamic
modeling of flexible bodies in multibody systems.