- Email: [email protected]
A NUMERICAL METHOD FOR STABILITY ANALYSIS OF PINNED FLEXIBLE MECHANISMS
Journal of Sound and Vibration (1996) 192(5), 941–957
A NUMERICAL METHOD FOR STABILITY ANALYSIS OF PINNED FLEXIBLE MECHANISMS D. G. B and S. W. L Department of Mechanical Engineering, 202 Ross Hall, Auburn University, AL 36849-5341, U.S.A. (Received 7 September 1994, and in final form 11 September 1995) A technique is presented to investigate the stability of mechanisms with pin-jointed flexible members. The method relies on a special floating frame from which elastic link co-ordinates are defined. Energies are easily developed for use in a Lagrange equation formulation, leading to a set of non-linear and mixed ordinary differential-algebraic equations of motion with constraints. Stability and bifurcation analysis is handled using a numerical procedure (generalized co-ordinate partitioning) that avoids the tedious and difficult task of analytically reducing the system of equations to a number equalling the system degrees of freedom. The proposed method was then applied to (1) a slider–crank mechanism with a flexible connecting rod and crank of constant rotational speed, and (2) a four-bar linkage with a flexible coupler with a constant speed crank. In both cases, a single pinned–pinned beam bending mode is employed to develop resonance curves and stability boundaries in the crank length–crank speed parameter plane. Flip and fold bifurcations are common occurrences in both mechanisms. The accuracy of the proposed method was also verified by comparison with previous experimental results [1]. 7 1996 Academic Press Limited
1. INTRODUCTION
Researchers and engineers have recognized that the dynamic behavior of modern-day, lightweight and high speed mechanisms are influenced by the elastic nature of the links. For example the slider–crank mechanisms, found in automotive engines and gas compressors, is a commonly investigated flexible mechanism. The focus of past investigations in flexible mechanisms has been for response (using either finite element or mode methods) and elastic stability, the latter being the concern of this paper. Most studies to date concerning elastic stability have dropped the non-linearities from the analysis, producing a system of ordinary differential equations with periodic coefficients such as the Mathieu equation; stability boundaries separated stable zones and unstable zones of exponentially increasing solutions. Stability zones were plotted in a parameter plane; for example, crank length versus crank speed. However, recent experimental [1, 2] and analytical studies [2–6] have indicated that the non-linearities are important contributors to accurate stability and bifurcation analysis. For example, period doubling —a distinctly non-linear phenomenon—occurs in both experiment and non-linear analysis. The primary purpose of this paper is to present a method for elastic stability in closed loop mechanisms with flexible beam links of constant cross-sectional area and undergoing in-plane motion, and with non-linearities in the differential equations of motion. One link of the mechanism rotates at a constant angular velocity. The technique is a numerical approach based on an assumed mode method. A floating frame between pin joints is 941 0022–460X/96/200941 + 17 $18.00/0
7 1996 Academic Press Limited
942
. . . .
chosen from which to easily describe energies for use in Lagrange’s equation. Constraints are developed from closed loop vector equations around the mechanism. The equations of motion are interfaced with the software AUTO [7] to determined bifurcations and resonance curves. The method has its advantages over existing approaches used to study stability in flexible mechanisms. It does not have limited accuracy for large parameters because, unlike perturbation methods, there is no dependence on small parameters. Prior publications that studied stability have all invested a substantial effort to reduce the equations that describe the constrained mechanical system to a minimal set of equations without constraints. For example, Viscomi and Ayres [8] and Badlani and Kleinhenz [9] presented a minimal set of unconstrained ODE’s to describe flexible rod slider–cranks, but in so doing introduce a number of physical assumptions to make those equations tractable; the technique presented here is a real time saver for the analyst because a minimal set of unconstrained ODE’s is not required. The method also produces more accurate stability boundaries because, as mentioned, fewer physical assumptions are required, and without removing non-linearities. When compared to the equations presented in reference [8], the equations presented here include the fuller effect of foreshortening on response, by not constraining the rod angle to be the same as would be generated by rigid body kinematics. This should, therefore, lead to more accurate response and stability boundaries than the equations of Viscomi and Ayres would produce, particularly for larger deflection. The aforementioned comparison is detailed in reference [2]. The technique is extendable to mechanisms with more than one flexible link, which has yet to be considered in the mechanism stability literature. Another important feature is the ease with which higher order deformations and their energies can be included in the analysis (e.g., foreshortening, defined here as the axial displacement of the beam tip as a result of the transverse displacement of the inextensible beam neutral axis, which is very important for stability analysis). Stability boundaries of a slider–crank compare favorably to experimental results. Lastly, the non-linear equation stability of a four-bar mechanism with a flexible coupler is presented for the first time in the literature. 2. EQUATIONS OF MOTION
A floating frame is defined to be located between the pinned ends of any flexible link. These floating frames are shown in Figures 1 and 2 as rigid massless links, pinned at one end and pinned-and-slotted at the other end, so that the slotted pin can translate along an imaginary rigid link axis. By attaching frames in a similar manner to successive links, it is possible to study systems with more than one flexible link. Note that such frames do not effect the motion of the flexible rod because they are massless and do not alter the system degrees of freedom, yet do provide a convenient frame from which elastic deformations can be described. Asada [10] proposed a pinned–pinned frame without the slot for the response of flexible robots in simple bending. In a flexible mechanism investigation, Beale and Lee [11] included link foreshortening by incorporating the slot, because foreshortening was a necessary component for parametric resonance studies. The presence of the slot also allows for the inclusion of other beam-in-plane deformations; for example, axial deformation of the link centroidal axis. Because the slot permits the length of the flexible member projected onto the floating frame to vary with load, simple trigonometric relationships derived from closed loop, rigid link mechanism geometry and kinematics are not valid. It should be noted that other floating frames of reference exist for different end conditions; see, for example, reference [12].
943
Figure 1. Fleixble beams and floating frames.
We first define the co-ordinates to be used in this paper in a very similar notation used by Shabana [13]. Bold face letters represent column vectors or matrices. Here, the column vector qi is the generalized co-ordinate of the ith flexible member which includes two sets of co-ordinates as shown in Figure 2: reference frame co-ordinates qir which includes the location Ri and orientation u i of a selected member reference frame xi − yi with respect iT T to global frame X − Y; elastic modal co-ordinates qif = (qiT fa qft ) , which describe the flexible deformation with respect to the body reference frame and includes the axial contribution qifa and the transverse component qift . The vector of co-ordinates qi of the ith flexible iT T iT iT iT iT T member can then be written as qi = (qiT qfa qft ) . r qf ) = (R u i The elastic deformation on the ith body is w = (u i, v i )T, where u i is the longitudinal displacement and v i is the transverse displacement. The longitudinal deformation u i is
Figure 2. The co-ordinates of a flexible link.
. . . .
944 i s
a result of foreshortening u (x, t) of the ith flexible member, and also axial deformation uai : u i = uai − usi ,
(1)
i s
and u (x, t) is usi = 12
g
x
(1v i/1x)2 dx.
(2)
0
uai and v i on the ith body by assumed modes method are i uai (x, t) = ciT a (x)qfa (t)
(3)
i v i(x, t) = ciT t (x)qft (t),
(4)
and where cia and cit are the mode shape functions for axial and transverse displacement and t is time. We then denote the composite vector of the system generalized co-ordinates by q = [q1T , q2T , . . . , qNT]T; here N is the total number of flexible members in the system. In Figure 1 is shown the ith floating frame and flexible member connected to the (i + 1)th frame and elastic member. The kinematic constraints between adjacent bodies can be written by the non-linear algebraic constraint equations f(q, t) = 0.
(5)
In Figure 2, the distance between O i and O i + 1 is L i + u i(L, t) = L i + uai (L, t) − usi (L, t), and L i is the undeformed length of the ith member. The non-linear dynamic equations governing the motion of flexible mechanisms are developed by a direct variational approach. A global position vector ri for an arbitrary point on flexible member i, as shown in Figure 2, is ri = Ri + Ai(x + wi ),
(6)
where Ai is the transformation matrix from components in the local co-ordinate system x i − y i to the global fixed co-ordinate system. x is the local position vector of the arbitrary point on the elastic member before deformation. The kinetic energy T i is T i = 12
g
r i(r˙ i )Tr˙ i dV i,
(7)
Vi
where r i and V i are the mass density and volume of the ith flexible body, respectively. The strain energy U i stored in the flexible beam is determined from U i = 12
g
i E(exx )2 dV i,
(8)
Vi
i where exx is the normal strain for large deflection beam bending (in this paper, we neglect the effect of shear strain) and is given as
i exx =
0 1
1uai 1 1v i + 1x 2 1x
2
−y
1 2v i , 1x 2
(9)
where y is the distance from the neutral axis of the beam. By substituting equation (9) into equation (8), one obtains the energy U i contributed by the transverse bending stiffness, linear axial stiffness, and higher order stiffness terms [14], as shown in Table 1. Other types
945
T 1 Expressions for typical potential energies Potential energy, U
Mathematical description
g 0 1
1ua dx 1x
L
1 2
Axial strain energy
EA
0
g 0 1
1 2v dx 1x 2 2
L
1 2
Bending strain energy
EI
0
1 2
Higher order term 1
g
01
1ua 1v dx 1x 1x 2
L
EA
0
EA
g 01 1 4
0
Periodic gas force
1v dx 1x 4
L
1 2
Higher order term 2
2
−Ps
of potential energy created by forces such as a gas force (or applied periodic force) on a piston in a yet to be discussed slider–crank example are also shown in Table 1. Lagrange’s equation of motion for a constrained mechanism is given by d (1T i/1q˙i )T − (1T i/1qi )T + (1U i/1qi )T = Fie − FiT q l, dt
i = 1, . . . , N,
(10)
where Fiq is the constraint Jacobian matrix, l is the vector of Lagrange multipliers, and Fie is the generalized non-conservative forces. Then, substituting T i in equation (7) and U i in Table 1 into equation (10) and combining for all N flexible members, one obtains the following equation of motion for the system: Mq¨ + Kq = FNL − FTq l + Fe + FV ,
(11)
where M is the system mass matrix, K is the system stiffness matrix, FNL is the vector of the non-linear elastic forces, Fe is the vector of all Fie forces and FV is the vector of the centripetal and Coriolis force components. FNL and FV are by-products of the left side of equation (10). The combination of equation (11) and the twice differentiated constraint equations (5) can be written in the form
$
M Fq
FTq 0
%$ % $
%
q¨ −Kq + FNL + Fe + Fv = , l g
(12)
where g = −(Fq q˙)q q˙ − 2Fqt q˙ − Ftt .
(13)
Equation (12) can be further streamlined to the first order form for use in AUTO:
&
I 0 0
0 M Fq
'&' &
0 FTq 0
'
q˙ Q Q = −Kq + FNL + Fe + Fv , l s g
(14)
946
. . . .
where l s = l and I is the identity matrix. Equation (14) is an unusual description of the equations of motion, at least for those who work in non-linear dynamics and are accustomed to differential equations in terms of a minimum number of co-ordinates, without constraints and without Lagrange multipliers. However, this form has distinct advantages. First, the energies and constraint equations can be easily written and hence the equations of motion easily formed. Second, the equations of motion are tractable; this is often not the case when written in minimum co-ordinates [11]. Lastly, no physical insight is required to produce tractable equations through techniques such as linearization, perturbation, arguments about the order of term, etc. These techniqus typically place limits on the applicable range and accuracy of the analysis. The stability analysis formulation presented here and based on equation (14) can therefore avoid some pitfalls that can occur when dealing with a minimum number of differential equations. 3. STABILITY ANALYSIS
The flexible mechanisms considered here are driven by actuators with a motion that is a prescribed function of time, leading to systems of equations with explicitly time-dependent and periodic coefficients. The stability treatment of such a system is based on the application of Floquet theory, after linearization. General n-degree-of-freedom non-linear, time periodic systems are governed by 2n non-linear first order differential equations, written as z˙(t) = f(z, t),
(15)
where f(t) = f(t + T) is an 2n × 1 column vector of non-linear functions, t is time, and T is the period of the time-dependent coefficients. Equations (15) are in a form that is most easily combined with AUTO, through a user-defined subroutine. Notice that equations (15) have no constraints, and hence describe a mechanical system entirely in terms of an independent or minimal co-ordinate set z. To use AUTO to determine the stability of equations (15), a user-defined subroutine calculates and returns to AUTO the vector z˙ (equalling the calculated right side of equation (15)); this calculation must be performed using values of z and t that AUTO has returned to the subroutine. Hence user-defined subroutines interact with AUTO in the same ways as user-defined subroutines interact with a Runge–Kutta numerical integration routine. Since the equations developed by the general technique and presented in this paper are in the form of constrained equations (14), the co-ordinates q and Q are not independent. Hence, constraint equations (5) must be used to find the dependent co-ordinates from in independent co-ordinates z values, that had been returned from AUTO to the user-defined subroutine. z is a subset of q and Q, chosen here as the flexible link mode co-ordinates positions and velocities; the remaining co-ordinates of q and Q are the dependent co-ordinates. Next, the matrix and right side of equations (14) can now be formed. Equations (14) are now solved for the unknown vector [q˙ Q l s ] using, for example, Gaussian elimination. Lastly, the independent co-ordinates derivatives z˙ are chosen from among the q˙ and Q , and are returned to AUTO. This aforementioned procedure—used numerically to reduce the system within the user-defined subroutine and isolate the independent co-ordinates and their derivatives—is called numerically implicit generalized co-ordinate partitioning [15]. By using and satisfying the constraint equations (through a Newton–Raphson solver) to determine the co-ordinates q and Q in the user-defined subroutine, this also assures that kinematic constraints are not violated [15]. This procedure is also discussed later in the paper.
947
Generally, the stability of an orbit is investigated by perturbing the motion from a periodic solution. Defining the periodic steady state solution as zss (t) and the perturbed motion as zp (t), z(t) = zss (t) + zp (t) are substituted into differential equations (15). After ignoring the higher order terms, 2n first order linear differential equations with periodic coefficients of the following form are obtained: z˙p = A(t)zp ,
(16)
where A(t) is a 2n × 2n matrix including periodic coefficients, and A(t) = A(t + T). C(t) is called the state transition matrix of equation (16), and satisfies C = A(t)C,
(17)
where C(0) = I, an identity matrix, and hence C(T) = C,
(18)
where C is called the Floquet transition matrix. The stability of equation (16) is completely determined by the eigenvalues jk (also called the characteristic multipliers) of matrix C. The characteristic exponents of the system are defined as ak 2 ibk , where ak =
1 ln =jk =, T
bk =
0 1
j 1 tan−1 kI , jkR T
(19, 20)
where jkI and jkR represent the imaginary and real parts of the eigenvalue. The stability criteria of the system are determined by the real part of the characteristic exponent ak . If ak Q 0 or =jk = Q 1, for k = 1, . . . , 2n, then the system is stable. For linear systems with periodic coefficients, stability zones are bordered in a parameter plane (e.g., crank length and crank speed) by boundaries with solutions of period T and 2T. Non-linear system stability is also dependent upon the parameter values, the solution bifurcating when a parameter combination changes value and crosses a boundary between a stable and unstable zone. Three different types of non-linear equation bifurcation can occur. If one of the characteristic multipliers jk leaves the unit circle in the complex plane through (+1,0), then a fold (also known as saddle node or jump) bifurcation occurs. If one of the jk leaves from (−1,0), then a flip (or period-doubling) bifurcation occurs. Otherwise, a pair complex conjugate multipliers of magnitude equaling unity produce a Hopf bifurcation when crossing the unit circle (see details in reference [16]). In this paper, the software AUTO [7] is used as a tool to analyze the stability and bifurcations in flexible mechanisms. Several major capabilities of the AUTO program are as follows: (1) it continues branches of stable and unstable solutions in parameter space; (2) the program detects and locates all types of bifurcations; (3) the program can switch branches at bifurcation points. Equation (14), including its Lagrange multipliers, is now coupled with AUTO to study stability. The following steps describe a co-ordinate partitioning algorithm that was incorporated in user-supplied subroutine FUNC, which is called numerous times by AUTO as it incrementally constructs a resonance curve. The algorithm is limited at this point in time to mechanical systems with degrees of freedom solely from the NT T elastic links. Elastic link variables qf and Qf , which are defined as [q1T and f , . . . , qf ] contained in the q vector, enter subroutine FUNC through the subroutine statement argument list. q˙f and Q f are returned to AUTO through the same subroutine FUNC argument list.
. . . .
948
Step 1. Given values qf and Qf , reference frame generalized co-ordinates qr and Qr (also called the dependent variables) are obtained from constraint equations (5) and their derivatives using a Newton–Raphson solver or explicitly derivable mathematical relationship. Step 2. Using q and Q calculated from step 1 in the left side matrix and the right side of equation (14), a Gaussian elimination routine is called to solve for [q˙T, Q˙T, l Ts ]T in equation (14). Step 3. q˙f and Q f are then returned through FUNC argument list. In this scheme, only those co-ordinates associated with the degree of freedom (the independent co-ordinates) in the system are handled by the AUTO program; the Lagrange multipliers and the dependent co-ordinates are confined to subroutine FUNC. 4. NUMERICAL EXAMPLES
The method is tested in the next section on a slider–crank mechanism with a flexible connecting rod and a four-bar linkage with a flexible coupler. Each mechanism is limited to planar motion, and the contributions to flexible member response of axial deformation of the centroidal axis length, shear deformation and rotatory inertia are not considered. One pinned–pinned beam bending mode is used to approximate rod deformation: previous experimental and analytical studies [1, 2] attest to the accuracy of this approximation in the slider–crank at crank speeds below the first natural frequency of the rod in bending, and in the presence of such non-linear phenomena as jumps and period-doubling bifurcations. 4.1. 1: – The slider–crank mechanism in Figure 3 has a flexible connecting rod of undeformed length, L, and mass per unit length, m¯ . It consists of a rigid crank with length (a), running at constant speed (v), and with a rigid sliding block of mass (Ms ). One rigid pinned and pinned-and-slotted floating frame is required to permit the flexible rod to freely bend and foreshorten, and is shown in Figure 3. Because both ends of the rod are pinned, the boundary conditions for bending are v(0, t) = 0,
v(L, t) = 0,
vxx (0, t) = 0,
vxx (L, t) = 0,
Figure 3. The flexible rod slider–crank mechanism with a floating frame.
(21)
949
where v is the rod transerse deflection measured from the rigid floating frame, and vxx is v differentiated twice with respect to x. Using the assumed mode method (with a pinned–pinned beam first mode shape) and separation of variables, v can be approximated as v(x, t) = q1(t) sin (px/L).
(22)
Three generalized co-ordinates are utilized to describe the motion of the system; they are q1 (transverse deflection of the bending rod), f (rod angle), and s (position of the sliding block). The position vector to any material point on the rod is, including foreshortening,
0 g
x
r¯ = aaˆ + x −
0
1
vx2 dx eˆ1 + veˆ2 , 2
(23)
where aˆ, eˆ1 and eˆ2 are unit vectors along vector OA, the x-axis and the y-axis of the floating frame, as shown in Figure 3. The kinetic energy T is the sum of the contributions from the flexible rod in equation (5), the crank and the piston, respectively: T=
m¯ 2
g
L
0
01
2
dr¯ dr¯ ds dx + Tcrank + 12 Ms , dt dt dt
(24)
where, dr¯ /dt is to be described as
0
1 0
1
dr¯ df df = av sin (f − vt) − v eˆ1 + av cos (f − vt) + x + vt eˆ2 . dt dt dt
(25)
Note that the foreshortening velocity has been removed because it has historically been considered insignificant; see, for example reference [6]. The potential energy U of the system includes the potential energy in bending and a conservative gas pressure force P, and is expressible as U=
EI 2
g
L 2 vxx dx − Ps.
(26)
0
The system is subject to the following two constraints, F1 and F2 , derivable by forming dot products of a vector loop equation around the mechanism in Figure 3 with unit vectors ıˆ and j :
0 g
1
(27)
vx2 dx +s − a − L = 0. 2
(28)
L
F1 0 a sin (vt) + sin (f) L −
0
0 g
L
F2 0 a cos (vt) + cos (f) L −
0
vx2 dx = 0, 2
1
Equations (22) are substituted into equation (24)–(28), and then into the Lagrange equation (10), to obtain the following equations of motion: 1 2
Lm¯
0 1
d2q1 L 2m¯ d2f p 2q1 sin f p 2q1 cos f df p 4EI l1 + l2 − 12 Lm¯ q1 + q 2 + 2 + dt p dt 2L 2L dt 2L 3 1 2
2m¯av 2 sin L (f − vt) = 0, p
+
(29)
. . . .
950
0 1
L 2m¯ d2q1 df dq + q1 Lm¯q1 1 + dt p dt 2 dt
0
+ L sin f −
0
1 3
L 3m¯ +
1
0
1
Lm¯ 2 d2f p 2q1 cos f q1 l1 2 − L cos f − 2 dt 4L
1
p 2q1 sin f m¯av 2L 2 sin (f − vt) 2m¯av 2Lq1 cos (f − vt) l2 + + = 0, 4L 2 p Ms
d2s − P − l2 = 0, dt 2
(30) (31)
with constraint equations
0
0
F2 0 a cos (vt) + cos (f) L −
1
p 2q12 = 0, 4L
(32)
p 2q12 + s − a − L = 0. 4L
(33)
F1 0 a sin (vt) + sin (f) L −
1
Several non-dimensional parameters are defined to non-dimensionalize equations (29)–(33): g = q1 /L,
R1 = a/L,
V 2 = m¯L 4v 2/p 4EI, L1 = l1 /m¯Lv 2,
u = vt,
M = Ms /m¯L,
F = P/m¯Lv 2,
S = s/L,
L2 = l2 /m¯Lv 2,
(34)
and the non-dimensional equations, now including non-dimensional damping m1 , are 2 1 4 g¨ + m1 g˙ + f + L1 p 2g sin f + L2 p 2g cos f − f 2g + 2 g + R1 sin (f − u) = 0, p V p g¨ + pf gg˙ +
0
1 0
1 0
(35)
1
p p 2 p3 p3 + g f + −p cos f + g cos f L1 + p cos f − g 2 sin f L2 3 2 4 4
p +2gR1 cos (f − u) + R1 sin (f − u) = 0, 2
(36)
MS − F − L2 = 0,
(37)
F1 = R1 sin u + (1 − p 2/4g 2 ) sin f = 0,
(38)
0
F2 = S + 1 −
1
p2 2 g cos f + R1 cos u − (1 + R1 ) = 0. 4
(39)
Equations (38) and (39) are, next, differentiated twice with respect to time. These equations, combined with equations (35)–(37), are now written in a non-dimensional form of equation (14) i.e., eight first order non-dimensionalized equations of the form z˙ = f(z, u),
(40)
= f , S = S , G = g˙ , L 1s = L1 , L 2s = L2 , and f is a where z = (f, S, g, f , S , G , L 1s , L 2s )T, f vector of non-linear functions. Bifurcation analyses are performed based on the algorithm described in section 3. Given flexible link co-ordinates g and G , constraint equation (38) and (39) and their derivations
951
are used to solve for rigid co-ordinates f, f , S and S . Their values are substituted into equation (40) and Gaussian elimination is employed to obtain the values of g˙ and G , which are returned to AUTO through the subroutine FUNC argument list. Two cases— M = 2·0 and M = 9·2—are presented to show the effect of M on stability boundaries in the R1–V plane. A non-dimensional damping ratio m1 = 0·085 was selected in each case. Although M = 9·2 is very large in a practical sense, this case is studied because it models experimental data.
4.1.1. Non-dimensionalized mass M = 2·0 Three resonance curves (with non-dimensional crank lengths of 0·04, 0·08 and 0·2) of non-dimensional amplitude versus non-dimensional crank speed are shown in Figure 4. For the 0·04 crank and increasing speed, the response jumps at V = 0·8034 to another stable period-1 solution, and then eventually period doubles through a flip bifurcation when V = 1·0714. For speeds greater than 1·0714, the curve displays only the unstable period-1 solution amplitude; the actual response is an unspecified sequence of period-doubled solutions, i.e., a period-doubling cascade of unspecified (but large) amplitude and of period 2n (T. In this study non-dimensional amplitudes beyond 0·15 will violate the assumptions of small deformation applicable to the Euler–Bernoulli beam theory used to derive the equations of motion. However, resonance curves will sometimes be presented that are of larger amplitude, simply to demonstrate the mathematical nature of the equations’ response. For the 0·08 crank, the response at V = 0·7322 jumps to an unstable period-double cascade instead of a stable period-1 solution. The 0·2 crank resonance curve presents lower speed bifurcations. A jump (saddle node) bifurcation occurs between V = 0·4203 and V = 0·3641. Like the aforementioned jumps, this one also has a softening spring characteristic. A zone of period doubling occurs between V of 0·4941 and 0·5785. The
Figure 4. Resonance curves with instability zones for M = 2. ———, Stable period-1; r– – – r, jump region; q – – – q, period-doubling cascade region.
952
. . . .
Figure 5. Non-linear equation instability zones in the non-dimensional crank length (R1 ) versus non-dimensional crank speed (V) plane, M = 2.
response of other solutions was not investigated here, the period-1 solution being shown as unstable within this region. Another means of presenting instability zones is shown in Figure 5, in which zones of instability between bifurcation boundaries are delineated in the (R1–V) parameter plane; i.e., a summary of resonance curve bifurcation points. This is similar to defining linear parametric resonance zones of instability, but based upon the non-linear equations. The instability boundaries are designated with the prefix P for period-doubling and J for jump. Bifurcations are numbered in order of occurrence (regardless of the letter prefix designation) as the backbone curve is traced along the period-1 resonance curve; imagine trekking the path of the resonance curve from low speed to high speed, along stable and unstable solution paths, and without jumping. For example, the sequence of bifurcations encountered while trekking the backbone curve at R1 = 0·15 is J1–J2–P3–P4–J5–J6–P7. 4.1.2. Non-dimensionalized mass M = 9·2 This case was the subject of an experimental and analytical study of response presented in references[1] and [2] that showed good correlation between them for a wide range of speeds. Here, those experimental results are compared to the analytically derived resonance curve and its bifurcations. In practice, this non-dimensional mass ratio would be excessive and makes for very complicated non-linear dynamics, but was a consequence of the desire for an inexpensive and easily constructed mechanism. The parameters of that apparatus were measured to be as follows: length of aluminum rod L = 11·5 in; rod cross-section 0·75 in × 0·0603 in; damping ratio = 0·02 (m1 = 0·085 in equations (35)); piston mass Ms = 0·478 lb, rod mass per unit length m¯ = 0·0045 lb/inch (M = 9·2); simply supported rod first natural frequency v = 38 Hz. In Figure 6, experimental resonance curves are shown from reference [1], with the analytical resonance boundaries superposed for comparison. For the 0·125 in crank (R1 = 0·0109) the first and only
953
Figure 6. Resonance curves of mid-span rod deflections from the non-linear equation and experiment, M = 9·2; crank lengths 0·125 in, 0·25 in, 0·5 in and 1 in. ———, Experimental; –––, analytical; r–––r, jump region; q–––q, period-doubling cascade region.
analytical instability is a jump at 1908 rpm; experimentally, the rod yielded just beyond the same speed. For the 0·25 in crank, the first analytical instability is a jump at 1075 rpm, and the experimental jump occurred just below that. Increasing speed caused both analytical and experimental response to decrease, until a period-doubling bifurcation occurred between 1400 rpm and 1467 rpm. In this region, the period-2 solutions exist with a large response for both analytical and experimental results. For the 0·5 in crank, the first analytical instability is a jump at 994 rpm, and the experimental jump occurred just below that speed. Increasing speed caused both the analytical and experimental response to decrease, until a period-doubling bifurcation at 1138 rpm. At 1080 rpm (amplitude not shown) a recorded experimental trace did become transient or chaotic, at which point the experiment was stopped. For the 1 in crank, the first encountered instability is at 690 rpm, where the response, both experimentally and analytically did jump to the top of the backbone curve. Experimental data collection ended at 750 rpm and before 780 rpm, where a period-doubling cascade is expected. 4.2. 2: - The response from the non-linear equations of a four-bar with flexible links has previously been studied by Sadler and Sandor [17] and Sung and Thompson [18]. Numerous parametric resonance studies based upon linear equation approximations have been performed, for example in references [19, 20], but no stability and bifurcation analysis has been attempted when non-linearities are not truncated. The geometry of a four-bar with a flexible coupler and floating frame is shown in Figure 7, with crank and follower assumed rigid. By application of the same aforementioned procedures as were used for the slider–crank, it was found that five equations (three second order dynamic equations in f2 , f3 and v, and two constraint equations formed from a vector loop around the mechanism dotted with ıˆ and j ) also described the system. For brevity, the derivation
954
. . . .
Figure 7. The flexible coupler four-bar linkage with a floating frame.
of equations of motion is not shown here. The interested reader can find the details in reference [21]. AUTO was again used to construct resonance curves and locate bifurcations. In Figure 8 are shown two resonance curves that demonstrate the effect of non-dimensional crank length. In this example, the four-bar has non-dimensional crank length R1 ( = a1 /L) equalling 0·04 and 0·08, R3 ( = a3 /L) = 0·70, R4 ( = a4 /L) = 0·90, where a1 , L, a3 and a4 are the lengths of the crank, undeformed coupler, follower and
Figure 8. Resonance curves of mid-span rod deflections with instability zones, R1 = 0·04 and 0·08, R3 = 0·7, R4 = 0·9, M = 2.
955
Figure 9. Non-linear equation instability zones in the non-dimensional crank length (R1 ) versus non-dimensional crank speed (V) plane.
baseline, respectively. The ratio of the mass of the follower to the coupler (M) is 2·0, and the non-dimensional damping m1 = 0·01. The material and cross-section properties of the flexible coupler are the same as the connecting rod in the slider–crank. The non-dimensional crank speed V is the ratio of the crank speed to a flexible coupler’s first natural frequency of 38 hz. The non-dimensional response g is the coupler mid-span deflection divided by the coupler length. The follower is treated as a rigid and slender rod, the mass moment of inertia of which about its center of mass is 1/12 of the follower mass times a32. The system response g and its stability are totally described through the equation of motion by non-dimensional parameters R1 , R3 , R4 , M, m1 and V. As would be expected, the larger crank size induces a larger bending response in the coupler at a given non-dimensional crank speed (defined as the ratio of crank speed to first natural frequency of the pinned–pinned coupler in bending). The jump instability boundaries with a softening spring, and a period-doubling boundary are shown in Figure 8, and also are compiled in the R1–V parameter plane in Figure 9. In Figure 10 are shown three resonance curves that demonstrate the effect of follower length. R3 has values of 0·5, 0·7 and 0·964286, with R1 = 0·06, R4 = 0·90714 and M = 0·8143. The figure again shows a softening spring character of the jump, and a response that amplifies with followers of larger length. 5. CONCLUSIONS
An easily implemented method to derive equations of motion for flexible linkage mechanisms (such as constant crank speed planar slider–cranks and four-bars) is combined with a numerical scheme to study interesting non-linear effects, such as bifurcations and stability boundaries in the crank length–crank speed parameter plane. The method hinges on both a simple floating frame scheme and a generalized co-ordinate partitioning
956
. . . .
Figure 10. Resonance curves for R3 = 0·5, 0·7 and 0·964286; R1 = 0·06, R4 = 0·90714, M = 0·8143. ———, Stable period-1; r. . . . . r, unstable period-1, jump; q— – — –q, unstable period-1, period-doubling cascade.
algorithm that was specifically designed to interact with software (AUTO) that is used for continuation, bifurcation and construction of resonance curves. The floating frame is pinned on one end, and pinned-and-slotted on the other to permit flexible link length changes which are so important for stability studies. Geometric and dynamic non-linearities are easily included by having sufficient energies as input to the Lagrange equations. These energies are easily derived because of the ease with which they are described with respect to the floating frame. The technique does not require that the system be reduced to a minimum number of differential equations; that is the constraint equations and Lagrange multipliers need not be removed from the analysis. This fact considerably simplifies the determination of instabilities and improves and accuracy of stability boundaries by eliminating the need to apply physical reasoning, intuition and/or assumptions to remove complicating terms that stand in the way of reducing the number and size of equations to a tractable set. Stability boundaries and resonance curves are presented for a slider–crank with a flexible rod and four-bar mechanism with a flexible coupler. In each case a single mode is used to describe the response and stability. Bifurcations included jumps with a softening spring character and period-doubling cascades. In the case of the slider–crank, the results compared favorably with the experiments. Future work will extend this method to situations with more than one flexible link and more than one elastic mode per flexible link. REFERENCES 1. D. H and D. G. B 1995 Nonlinear Dynamics 7, 365–384. Experimental observations of a flexible slider crank mechanism at very high speeds. 2. D. G. B and S. W. L 1993 ASME Design Conference DE-Vol. 56, 123–130. Linear and nonlinear instabilities in a flexible mechanism. 3. D. G. B and S. W. L 1996 Mechanism and Machine Theory 31, 215–227. Nonlinear equation instability boundaries in flexible mechanisms.
957
4. S. W. L and D. G. B 1995 Zeitschrift fur Angewandte Mathematik und Mechanik 75, 883–888. Periodic solutions and stability of a nonlinear flexible slider crank mechanism by HBM/FFT. 5. S. H and S. W. S 1993 Nonlinear Dynamics 4, 573–603. The dynamic stability and nonlinear resonance of a flexible connecting rod: continuous parameter model. 6. S. H and S. W. S 1994 Journal of Sound and Vibration 170, 25–49. The dynamic stability and nonlinear resonance of a flexible connecting rod: single-mode model. 7. E. D and H. B. K 1991 International Journal of Bifurcation and Chaos 1, 493–520. Numerical analysis and control of bifurcation problems (I) bifurcation in finite dimensions. 8. B. V. V and R. S. A 1971 Transactions of the American Society of Mechanical Engineers, Journal of Engineering Industry 93, 251–262. Nonlinear dynamic response of elastic slider–crank mechanism. 9. M. B and W. K 1979 Transactions of the American Society of Mechanical Engineers, Journal of Mechanical Design 101, 149–153. Dynamic stability of elastic mechanisms. 10. H. A, Z. D. M and H. T 1990 Transactions of the American Society of Mechanical Engineers, Journal of Dynamic Systems, Measurement, and Control 112, 177–185. Inverse dynamics of flexible robot arms: modeling and computation for trajectory control. 11. D. G. B and S. W. L 1991 ASME Design Conference DE-Vol. 37, 161–166. Investigation of parametric resonance instability in a flexible rod slider crank mechanism. 12. O. P. A and A. A. S 1985 Computers and Structure 21, 1303–1312. Dynamic analysis of multibody systems using component modes. 13. A. A. S 1989 Dynamics of Multibody Systems. New York: John Wiley. 14. J. S. P 1968 Theory of Matrix Structural Analysis. New York: McGraw-Hill. 15. E. J. H 1989 Computer-aided Kinematics and Dynamics of Mechanical Systems. Boston, Massachussetts: Allyn and Bacon. 16. J. G and P. H 1983 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. 17. J. P. S and G. N. S 1974 Transactions of the American Society of Mechanical Engineers, Journal of Engineering for Industry 96, 411–419. Nonlinear vibration analysis of elastic four-bar linkages. 18. C. K. S, B. S. T, T. M. X and C. H. W 1986 Mechanism and Machine Theory 21, 121–133. An experimental study on the nonlinear elastodynamic response of linkage mechanisms. 19. M. R. S and L. M 1971 Journal of Mechanical Engineering Science 13, 237–242. Stability of a four-bar linkage with flexible coupler. 20. W. L. C, B. T and R. G. F 1984 Mechanism and Machine Theory 19, 307–317. Critical running speeds and stability of high-speed flexible mechanisms. 21. S. W. L 1994 Ph.D. Thesis, Auburn University. Elastic stability and active control of a flexible rod mechanism.
A NUMERICAL METHOD FOR STABILITY ANALYSIS OF PINNED FLEXIBLE MECHANISMS D. G. B and S. W. L Department of Mechanical Engineering, 202 Ross Hall, Auburn University, AL 36849-5341, U.S.A. (Received 7 September 1994, and in final form 11 September 1995) A technique is presented to investigate the stability of mechanisms with pin-jointed flexible members. The method relies on a special floating frame from which elastic link co-ordinates are defined. Energies are easily developed for use in a Lagrange equation formulation, leading to a set of non-linear and mixed ordinary differential-algebraic equations of motion with constraints. Stability and bifurcation analysis is handled using a numerical procedure (generalized co-ordinate partitioning) that avoids the tedious and difficult task of analytically reducing the system of equations to a number equalling the system degrees of freedom. The proposed method was then applied to (1) a slider–crank mechanism with a flexible connecting rod and crank of constant rotational speed, and (2) a four-bar linkage with a flexible coupler with a constant speed crank. In both cases, a single pinned–pinned beam bending mode is employed to develop resonance curves and stability boundaries in the crank length–crank speed parameter plane. Flip and fold bifurcations are common occurrences in both mechanisms. The accuracy of the proposed method was also verified by comparison with previous experimental results [1]. 7 1996 Academic Press Limited
1. INTRODUCTION
Researchers and engineers have recognized that the dynamic behavior of modern-day, lightweight and high speed mechanisms are influenced by the elastic nature of the links. For example the slider–crank mechanisms, found in automotive engines and gas compressors, is a commonly investigated flexible mechanism. The focus of past investigations in flexible mechanisms has been for response (using either finite element or mode methods) and elastic stability, the latter being the concern of this paper. Most studies to date concerning elastic stability have dropped the non-linearities from the analysis, producing a system of ordinary differential equations with periodic coefficients such as the Mathieu equation; stability boundaries separated stable zones and unstable zones of exponentially increasing solutions. Stability zones were plotted in a parameter plane; for example, crank length versus crank speed. However, recent experimental [1, 2] and analytical studies [2–6] have indicated that the non-linearities are important contributors to accurate stability and bifurcation analysis. For example, period doubling —a distinctly non-linear phenomenon—occurs in both experiment and non-linear analysis. The primary purpose of this paper is to present a method for elastic stability in closed loop mechanisms with flexible beam links of constant cross-sectional area and undergoing in-plane motion, and with non-linearities in the differential equations of motion. One link of the mechanism rotates at a constant angular velocity. The technique is a numerical approach based on an assumed mode method. A floating frame between pin joints is 941 0022–460X/96/200941 + 17 $18.00/0
7 1996 Academic Press Limited
942
. . . .
chosen from which to easily describe energies for use in Lagrange’s equation. Constraints are developed from closed loop vector equations around the mechanism. The equations of motion are interfaced with the software AUTO [7] to determined bifurcations and resonance curves. The method has its advantages over existing approaches used to study stability in flexible mechanisms. It does not have limited accuracy for large parameters because, unlike perturbation methods, there is no dependence on small parameters. Prior publications that studied stability have all invested a substantial effort to reduce the equations that describe the constrained mechanical system to a minimal set of equations without constraints. For example, Viscomi and Ayres [8] and Badlani and Kleinhenz [9] presented a minimal set of unconstrained ODE’s to describe flexible rod slider–cranks, but in so doing introduce a number of physical assumptions to make those equations tractable; the technique presented here is a real time saver for the analyst because a minimal set of unconstrained ODE’s is not required. The method also produces more accurate stability boundaries because, as mentioned, fewer physical assumptions are required, and without removing non-linearities. When compared to the equations presented in reference [8], the equations presented here include the fuller effect of foreshortening on response, by not constraining the rod angle to be the same as would be generated by rigid body kinematics. This should, therefore, lead to more accurate response and stability boundaries than the equations of Viscomi and Ayres would produce, particularly for larger deflection. The aforementioned comparison is detailed in reference [2]. The technique is extendable to mechanisms with more than one flexible link, which has yet to be considered in the mechanism stability literature. Another important feature is the ease with which higher order deformations and their energies can be included in the analysis (e.g., foreshortening, defined here as the axial displacement of the beam tip as a result of the transverse displacement of the inextensible beam neutral axis, which is very important for stability analysis). Stability boundaries of a slider–crank compare favorably to experimental results. Lastly, the non-linear equation stability of a four-bar mechanism with a flexible coupler is presented for the first time in the literature. 2. EQUATIONS OF MOTION
A floating frame is defined to be located between the pinned ends of any flexible link. These floating frames are shown in Figures 1 and 2 as rigid massless links, pinned at one end and pinned-and-slotted at the other end, so that the slotted pin can translate along an imaginary rigid link axis. By attaching frames in a similar manner to successive links, it is possible to study systems with more than one flexible link. Note that such frames do not effect the motion of the flexible rod because they are massless and do not alter the system degrees of freedom, yet do provide a convenient frame from which elastic deformations can be described. Asada [10] proposed a pinned–pinned frame without the slot for the response of flexible robots in simple bending. In a flexible mechanism investigation, Beale and Lee [11] included link foreshortening by incorporating the slot, because foreshortening was a necessary component for parametric resonance studies. The presence of the slot also allows for the inclusion of other beam-in-plane deformations; for example, axial deformation of the link centroidal axis. Because the slot permits the length of the flexible member projected onto the floating frame to vary with load, simple trigonometric relationships derived from closed loop, rigid link mechanism geometry and kinematics are not valid. It should be noted that other floating frames of reference exist for different end conditions; see, for example, reference [12].
943
Figure 1. Fleixble beams and floating frames.
We first define the co-ordinates to be used in this paper in a very similar notation used by Shabana [13]. Bold face letters represent column vectors or matrices. Here, the column vector qi is the generalized co-ordinate of the ith flexible member which includes two sets of co-ordinates as shown in Figure 2: reference frame co-ordinates qir which includes the location Ri and orientation u i of a selected member reference frame xi − yi with respect iT T to global frame X − Y; elastic modal co-ordinates qif = (qiT fa qft ) , which describe the flexible deformation with respect to the body reference frame and includes the axial contribution qifa and the transverse component qift . The vector of co-ordinates qi of the ith flexible iT T iT iT iT iT T member can then be written as qi = (qiT qfa qft ) . r qf ) = (R u i The elastic deformation on the ith body is w = (u i, v i )T, where u i is the longitudinal displacement and v i is the transverse displacement. The longitudinal deformation u i is
Figure 2. The co-ordinates of a flexible link.
. . . .
944 i s
a result of foreshortening u (x, t) of the ith flexible member, and also axial deformation uai : u i = uai − usi ,
(1)
i s
and u (x, t) is usi = 12
g
x
(1v i/1x)2 dx.
(2)
0
uai and v i on the ith body by assumed modes method are i uai (x, t) = ciT a (x)qfa (t)
(3)
i v i(x, t) = ciT t (x)qft (t),
(4)
and where cia and cit are the mode shape functions for axial and transverse displacement and t is time. We then denote the composite vector of the system generalized co-ordinates by q = [q1T , q2T , . . . , qNT]T; here N is the total number of flexible members in the system. In Figure 1 is shown the ith floating frame and flexible member connected to the (i + 1)th frame and elastic member. The kinematic constraints between adjacent bodies can be written by the non-linear algebraic constraint equations f(q, t) = 0.
(5)
In Figure 2, the distance between O i and O i + 1 is L i + u i(L, t) = L i + uai (L, t) − usi (L, t), and L i is the undeformed length of the ith member. The non-linear dynamic equations governing the motion of flexible mechanisms are developed by a direct variational approach. A global position vector ri for an arbitrary point on flexible member i, as shown in Figure 2, is ri = Ri + Ai(x + wi ),
(6)
where Ai is the transformation matrix from components in the local co-ordinate system x i − y i to the global fixed co-ordinate system. x is the local position vector of the arbitrary point on the elastic member before deformation. The kinetic energy T i is T i = 12
g
r i(r˙ i )Tr˙ i dV i,
(7)
Vi
where r i and V i are the mass density and volume of the ith flexible body, respectively. The strain energy U i stored in the flexible beam is determined from U i = 12
g
i E(exx )2 dV i,
(8)
Vi
i where exx is the normal strain for large deflection beam bending (in this paper, we neglect the effect of shear strain) and is given as
i exx =
0 1
1uai 1 1v i + 1x 2 1x
2
−y
1 2v i , 1x 2
(9)
where y is the distance from the neutral axis of the beam. By substituting equation (9) into equation (8), one obtains the energy U i contributed by the transverse bending stiffness, linear axial stiffness, and higher order stiffness terms [14], as shown in Table 1. Other types
945
T 1 Expressions for typical potential energies Potential energy, U
Mathematical description
g 0 1
1ua dx 1x
L
1 2
Axial strain energy
EA
0
g 0 1
1 2v dx 1x 2 2
L
1 2
Bending strain energy
EI
0
1 2
Higher order term 1
g
01
1ua 1v dx 1x 1x 2
L
EA
0
EA
g 01 1 4
0
Periodic gas force
1v dx 1x 4
L
1 2
Higher order term 2
2
−Ps
of potential energy created by forces such as a gas force (or applied periodic force) on a piston in a yet to be discussed slider–crank example are also shown in Table 1. Lagrange’s equation of motion for a constrained mechanism is given by d (1T i/1q˙i )T − (1T i/1qi )T + (1U i/1qi )T = Fie − FiT q l, dt
i = 1, . . . , N,
(10)
where Fiq is the constraint Jacobian matrix, l is the vector of Lagrange multipliers, and Fie is the generalized non-conservative forces. Then, substituting T i in equation (7) and U i in Table 1 into equation (10) and combining for all N flexible members, one obtains the following equation of motion for the system: Mq¨ + Kq = FNL − FTq l + Fe + FV ,
(11)
where M is the system mass matrix, K is the system stiffness matrix, FNL is the vector of the non-linear elastic forces, Fe is the vector of all Fie forces and FV is the vector of the centripetal and Coriolis force components. FNL and FV are by-products of the left side of equation (10). The combination of equation (11) and the twice differentiated constraint equations (5) can be written in the form
$
M Fq
FTq 0
%$ % $
%
q¨ −Kq + FNL + Fe + Fv = , l g
(12)
where g = −(Fq q˙)q q˙ − 2Fqt q˙ − Ftt .
(13)
Equation (12) can be further streamlined to the first order form for use in AUTO:
&
I 0 0
0 M Fq
'&' &
0 FTq 0
'
q˙ Q Q = −Kq + FNL + Fe + Fv , l s g
(14)
946
. . . .
where l s = l and I is the identity matrix. Equation (14) is an unusual description of the equations of motion, at least for those who work in non-linear dynamics and are accustomed to differential equations in terms of a minimum number of co-ordinates, without constraints and without Lagrange multipliers. However, this form has distinct advantages. First, the energies and constraint equations can be easily written and hence the equations of motion easily formed. Second, the equations of motion are tractable; this is often not the case when written in minimum co-ordinates [11]. Lastly, no physical insight is required to produce tractable equations through techniques such as linearization, perturbation, arguments about the order of term, etc. These techniqus typically place limits on the applicable range and accuracy of the analysis. The stability analysis formulation presented here and based on equation (14) can therefore avoid some pitfalls that can occur when dealing with a minimum number of differential equations. 3. STABILITY ANALYSIS
The flexible mechanisms considered here are driven by actuators with a motion that is a prescribed function of time, leading to systems of equations with explicitly time-dependent and periodic coefficients. The stability treatment of such a system is based on the application of Floquet theory, after linearization. General n-degree-of-freedom non-linear, time periodic systems are governed by 2n non-linear first order differential equations, written as z˙(t) = f(z, t),
(15)
where f(t) = f(t + T) is an 2n × 1 column vector of non-linear functions, t is time, and T is the period of the time-dependent coefficients. Equations (15) are in a form that is most easily combined with AUTO, through a user-defined subroutine. Notice that equations (15) have no constraints, and hence describe a mechanical system entirely in terms of an independent or minimal co-ordinate set z. To use AUTO to determine the stability of equations (15), a user-defined subroutine calculates and returns to AUTO the vector z˙ (equalling the calculated right side of equation (15)); this calculation must be performed using values of z and t that AUTO has returned to the subroutine. Hence user-defined subroutines interact with AUTO in the same ways as user-defined subroutines interact with a Runge–Kutta numerical integration routine. Since the equations developed by the general technique and presented in this paper are in the form of constrained equations (14), the co-ordinates q and Q are not independent. Hence, constraint equations (5) must be used to find the dependent co-ordinates from in independent co-ordinates z values, that had been returned from AUTO to the user-defined subroutine. z is a subset of q and Q, chosen here as the flexible link mode co-ordinates positions and velocities; the remaining co-ordinates of q and Q are the dependent co-ordinates. Next, the matrix and right side of equations (14) can now be formed. Equations (14) are now solved for the unknown vector [q˙ Q l s ] using, for example, Gaussian elimination. Lastly, the independent co-ordinates derivatives z˙ are chosen from among the q˙ and Q , and are returned to AUTO. This aforementioned procedure—used numerically to reduce the system within the user-defined subroutine and isolate the independent co-ordinates and their derivatives—is called numerically implicit generalized co-ordinate partitioning [15]. By using and satisfying the constraint equations (through a Newton–Raphson solver) to determine the co-ordinates q and Q in the user-defined subroutine, this also assures that kinematic constraints are not violated [15]. This procedure is also discussed later in the paper.
947
Generally, the stability of an orbit is investigated by perturbing the motion from a periodic solution. Defining the periodic steady state solution as zss (t) and the perturbed motion as zp (t), z(t) = zss (t) + zp (t) are substituted into differential equations (15). After ignoring the higher order terms, 2n first order linear differential equations with periodic coefficients of the following form are obtained: z˙p = A(t)zp ,
(16)
where A(t) is a 2n × 2n matrix including periodic coefficients, and A(t) = A(t + T). C(t) is called the state transition matrix of equation (16), and satisfies C = A(t)C,
(17)
where C(0) = I, an identity matrix, and hence C(T) = C,
(18)
where C is called the Floquet transition matrix. The stability of equation (16) is completely determined by the eigenvalues jk (also called the characteristic multipliers) of matrix C. The characteristic exponents of the system are defined as ak 2 ibk , where ak =
1 ln =jk =, T
bk =
0 1
j 1 tan−1 kI , jkR T
(19, 20)
where jkI and jkR represent the imaginary and real parts of the eigenvalue. The stability criteria of the system are determined by the real part of the characteristic exponent ak . If ak Q 0 or =jk = Q 1, for k = 1, . . . , 2n, then the system is stable. For linear systems with periodic coefficients, stability zones are bordered in a parameter plane (e.g., crank length and crank speed) by boundaries with solutions of period T and 2T. Non-linear system stability is also dependent upon the parameter values, the solution bifurcating when a parameter combination changes value and crosses a boundary between a stable and unstable zone. Three different types of non-linear equation bifurcation can occur. If one of the characteristic multipliers jk leaves the unit circle in the complex plane through (+1,0), then a fold (also known as saddle node or jump) bifurcation occurs. If one of the jk leaves from (−1,0), then a flip (or period-doubling) bifurcation occurs. Otherwise, a pair complex conjugate multipliers of magnitude equaling unity produce a Hopf bifurcation when crossing the unit circle (see details in reference [16]). In this paper, the software AUTO [7] is used as a tool to analyze the stability and bifurcations in flexible mechanisms. Several major capabilities of the AUTO program are as follows: (1) it continues branches of stable and unstable solutions in parameter space; (2) the program detects and locates all types of bifurcations; (3) the program can switch branches at bifurcation points. Equation (14), including its Lagrange multipliers, is now coupled with AUTO to study stability. The following steps describe a co-ordinate partitioning algorithm that was incorporated in user-supplied subroutine FUNC, which is called numerous times by AUTO as it incrementally constructs a resonance curve. The algorithm is limited at this point in time to mechanical systems with degrees of freedom solely from the NT T elastic links. Elastic link variables qf and Qf , which are defined as [q1T and f , . . . , qf ] contained in the q vector, enter subroutine FUNC through the subroutine statement argument list. q˙f and Q f are returned to AUTO through the same subroutine FUNC argument list.
. . . .
948
Step 1. Given values qf and Qf , reference frame generalized co-ordinates qr and Qr (also called the dependent variables) are obtained from constraint equations (5) and their derivatives using a Newton–Raphson solver or explicitly derivable mathematical relationship. Step 2. Using q and Q calculated from step 1 in the left side matrix and the right side of equation (14), a Gaussian elimination routine is called to solve for [q˙T, Q˙T, l Ts ]T in equation (14). Step 3. q˙f and Q f are then returned through FUNC argument list. In this scheme, only those co-ordinates associated with the degree of freedom (the independent co-ordinates) in the system are handled by the AUTO program; the Lagrange multipliers and the dependent co-ordinates are confined to subroutine FUNC. 4. NUMERICAL EXAMPLES
The method is tested in the next section on a slider–crank mechanism with a flexible connecting rod and a four-bar linkage with a flexible coupler. Each mechanism is limited to planar motion, and the contributions to flexible member response of axial deformation of the centroidal axis length, shear deformation and rotatory inertia are not considered. One pinned–pinned beam bending mode is used to approximate rod deformation: previous experimental and analytical studies [1, 2] attest to the accuracy of this approximation in the slider–crank at crank speeds below the first natural frequency of the rod in bending, and in the presence of such non-linear phenomena as jumps and period-doubling bifurcations. 4.1. 1: – The slider–crank mechanism in Figure 3 has a flexible connecting rod of undeformed length, L, and mass per unit length, m¯ . It consists of a rigid crank with length (a), running at constant speed (v), and with a rigid sliding block of mass (Ms ). One rigid pinned and pinned-and-slotted floating frame is required to permit the flexible rod to freely bend and foreshorten, and is shown in Figure 3. Because both ends of the rod are pinned, the boundary conditions for bending are v(0, t) = 0,
v(L, t) = 0,
vxx (0, t) = 0,
vxx (L, t) = 0,
Figure 3. The flexible rod slider–crank mechanism with a floating frame.
(21)
949
where v is the rod transerse deflection measured from the rigid floating frame, and vxx is v differentiated twice with respect to x. Using the assumed mode method (with a pinned–pinned beam first mode shape) and separation of variables, v can be approximated as v(x, t) = q1(t) sin (px/L).
(22)
Three generalized co-ordinates are utilized to describe the motion of the system; they are q1 (transverse deflection of the bending rod), f (rod angle), and s (position of the sliding block). The position vector to any material point on the rod is, including foreshortening,
0 g
x
r¯ = aaˆ + x −
0
1
vx2 dx eˆ1 + veˆ2 , 2
(23)
where aˆ, eˆ1 and eˆ2 are unit vectors along vector OA, the x-axis and the y-axis of the floating frame, as shown in Figure 3. The kinetic energy T is the sum of the contributions from the flexible rod in equation (5), the crank and the piston, respectively: T=
m¯ 2
g
L
0
01
2
dr¯ dr¯ ds dx + Tcrank + 12 Ms , dt dt dt
(24)
where, dr¯ /dt is to be described as
0
1 0
1
dr¯ df df = av sin (f − vt) − v eˆ1 + av cos (f − vt) + x + vt eˆ2 . dt dt dt
(25)
Note that the foreshortening velocity has been removed because it has historically been considered insignificant; see, for example reference [6]. The potential energy U of the system includes the potential energy in bending and a conservative gas pressure force P, and is expressible as U=
EI 2
g
L 2 vxx dx − Ps.
(26)
0
The system is subject to the following two constraints, F1 and F2 , derivable by forming dot products of a vector loop equation around the mechanism in Figure 3 with unit vectors ıˆ and j :
0 g
1
(27)
vx2 dx +s − a − L = 0. 2
(28)
L
F1 0 a sin (vt) + sin (f) L −
0
0 g
L
F2 0 a cos (vt) + cos (f) L −
0
vx2 dx = 0, 2
1
Equations (22) are substituted into equation (24)–(28), and then into the Lagrange equation (10), to obtain the following equations of motion: 1 2
Lm¯
0 1
d2q1 L 2m¯ d2f p 2q1 sin f p 2q1 cos f df p 4EI l1 + l2 − 12 Lm¯ q1 + q 2 + 2 + dt p dt 2L 2L dt 2L 3 1 2
2m¯av 2 sin L (f − vt) = 0, p
+
(29)
. . . .
950
0 1
L 2m¯ d2q1 df dq + q1 Lm¯q1 1 + dt p dt 2 dt
0
+ L sin f −
0
1 3
L 3m¯ +
1
0
1
Lm¯ 2 d2f p 2q1 cos f q1 l1 2 − L cos f − 2 dt 4L
1
p 2q1 sin f m¯av 2L 2 sin (f − vt) 2m¯av 2Lq1 cos (f − vt) l2 + + = 0, 4L 2 p Ms
d2s − P − l2 = 0, dt 2
(30) (31)
with constraint equations
0
0
F2 0 a cos (vt) + cos (f) L −
1
p 2q12 = 0, 4L
(32)
p 2q12 + s − a − L = 0. 4L
(33)
F1 0 a sin (vt) + sin (f) L −
1
Several non-dimensional parameters are defined to non-dimensionalize equations (29)–(33): g = q1 /L,
R1 = a/L,
V 2 = m¯L 4v 2/p 4EI, L1 = l1 /m¯Lv 2,
u = vt,
M = Ms /m¯L,
F = P/m¯Lv 2,
S = s/L,
L2 = l2 /m¯Lv 2,
(34)
and the non-dimensional equations, now including non-dimensional damping m1 , are 2 1 4 g¨ + m1 g˙ + f + L1 p 2g sin f + L2 p 2g cos f − f 2g + 2 g + R1 sin (f − u) = 0, p V p g¨ + pf gg˙ +
0
1 0
1 0
(35)
1
p p 2 p3 p3 + g f + −p cos f + g cos f L1 + p cos f − g 2 sin f L2 3 2 4 4
p +2gR1 cos (f − u) + R1 sin (f − u) = 0, 2
(36)
MS − F − L2 = 0,
(37)
F1 = R1 sin u + (1 − p 2/4g 2 ) sin f = 0,
(38)
0
F2 = S + 1 −
1
p2 2 g cos f + R1 cos u − (1 + R1 ) = 0. 4
(39)
Equations (38) and (39) are, next, differentiated twice with respect to time. These equations, combined with equations (35)–(37), are now written in a non-dimensional form of equation (14) i.e., eight first order non-dimensionalized equations of the form z˙ = f(z, u),
(40)
= f , S = S , G = g˙ , L 1s = L1 , L 2s = L2 , and f is a where z = (f, S, g, f , S , G , L 1s , L 2s )T, f vector of non-linear functions. Bifurcation analyses are performed based on the algorithm described in section 3. Given flexible link co-ordinates g and G , constraint equation (38) and (39) and their derivations
951
are used to solve for rigid co-ordinates f, f , S and S . Their values are substituted into equation (40) and Gaussian elimination is employed to obtain the values of g˙ and G , which are returned to AUTO through the subroutine FUNC argument list. Two cases— M = 2·0 and M = 9·2—are presented to show the effect of M on stability boundaries in the R1–V plane. A non-dimensional damping ratio m1 = 0·085 was selected in each case. Although M = 9·2 is very large in a practical sense, this case is studied because it models experimental data.
4.1.1. Non-dimensionalized mass M = 2·0 Three resonance curves (with non-dimensional crank lengths of 0·04, 0·08 and 0·2) of non-dimensional amplitude versus non-dimensional crank speed are shown in Figure 4. For the 0·04 crank and increasing speed, the response jumps at V = 0·8034 to another stable period-1 solution, and then eventually period doubles through a flip bifurcation when V = 1·0714. For speeds greater than 1·0714, the curve displays only the unstable period-1 solution amplitude; the actual response is an unspecified sequence of period-doubled solutions, i.e., a period-doubling cascade of unspecified (but large) amplitude and of period 2n (T. In this study non-dimensional amplitudes beyond 0·15 will violate the assumptions of small deformation applicable to the Euler–Bernoulli beam theory used to derive the equations of motion. However, resonance curves will sometimes be presented that are of larger amplitude, simply to demonstrate the mathematical nature of the equations’ response. For the 0·08 crank, the response at V = 0·7322 jumps to an unstable period-double cascade instead of a stable period-1 solution. The 0·2 crank resonance curve presents lower speed bifurcations. A jump (saddle node) bifurcation occurs between V = 0·4203 and V = 0·3641. Like the aforementioned jumps, this one also has a softening spring characteristic. A zone of period doubling occurs between V of 0·4941 and 0·5785. The
Figure 4. Resonance curves with instability zones for M = 2. ———, Stable period-1; r– – – r, jump region; q – – – q, period-doubling cascade region.
952
. . . .
Figure 5. Non-linear equation instability zones in the non-dimensional crank length (R1 ) versus non-dimensional crank speed (V) plane, M = 2.
response of other solutions was not investigated here, the period-1 solution being shown as unstable within this region. Another means of presenting instability zones is shown in Figure 5, in which zones of instability between bifurcation boundaries are delineated in the (R1–V) parameter plane; i.e., a summary of resonance curve bifurcation points. This is similar to defining linear parametric resonance zones of instability, but based upon the non-linear equations. The instability boundaries are designated with the prefix P for period-doubling and J for jump. Bifurcations are numbered in order of occurrence (regardless of the letter prefix designation) as the backbone curve is traced along the period-1 resonance curve; imagine trekking the path of the resonance curve from low speed to high speed, along stable and unstable solution paths, and without jumping. For example, the sequence of bifurcations encountered while trekking the backbone curve at R1 = 0·15 is J1–J2–P3–P4–J5–J6–P7. 4.1.2. Non-dimensionalized mass M = 9·2 This case was the subject of an experimental and analytical study of response presented in references[1] and [2] that showed good correlation between them for a wide range of speeds. Here, those experimental results are compared to the analytically derived resonance curve and its bifurcations. In practice, this non-dimensional mass ratio would be excessive and makes for very complicated non-linear dynamics, but was a consequence of the desire for an inexpensive and easily constructed mechanism. The parameters of that apparatus were measured to be as follows: length of aluminum rod L = 11·5 in; rod cross-section 0·75 in × 0·0603 in; damping ratio = 0·02 (m1 = 0·085 in equations (35)); piston mass Ms = 0·478 lb, rod mass per unit length m¯ = 0·0045 lb/inch (M = 9·2); simply supported rod first natural frequency v = 38 Hz. In Figure 6, experimental resonance curves are shown from reference [1], with the analytical resonance boundaries superposed for comparison. For the 0·125 in crank (R1 = 0·0109) the first and only
953
Figure 6. Resonance curves of mid-span rod deflections from the non-linear equation and experiment, M = 9·2; crank lengths 0·125 in, 0·25 in, 0·5 in and 1 in. ———, Experimental; –––, analytical; r–––r, jump region; q–––q, period-doubling cascade region.
analytical instability is a jump at 1908 rpm; experimentally, the rod yielded just beyond the same speed. For the 0·25 in crank, the first analytical instability is a jump at 1075 rpm, and the experimental jump occurred just below that. Increasing speed caused both analytical and experimental response to decrease, until a period-doubling bifurcation occurred between 1400 rpm and 1467 rpm. In this region, the period-2 solutions exist with a large response for both analytical and experimental results. For the 0·5 in crank, the first analytical instability is a jump at 994 rpm, and the experimental jump occurred just below that speed. Increasing speed caused both the analytical and experimental response to decrease, until a period-doubling bifurcation at 1138 rpm. At 1080 rpm (amplitude not shown) a recorded experimental trace did become transient or chaotic, at which point the experiment was stopped. For the 1 in crank, the first encountered instability is at 690 rpm, where the response, both experimentally and analytically did jump to the top of the backbone curve. Experimental data collection ended at 750 rpm and before 780 rpm, where a period-doubling cascade is expected. 4.2. 2: - The response from the non-linear equations of a four-bar with flexible links has previously been studied by Sadler and Sandor [17] and Sung and Thompson [18]. Numerous parametric resonance studies based upon linear equation approximations have been performed, for example in references [19, 20], but no stability and bifurcation analysis has been attempted when non-linearities are not truncated. The geometry of a four-bar with a flexible coupler and floating frame is shown in Figure 7, with crank and follower assumed rigid. By application of the same aforementioned procedures as were used for the slider–crank, it was found that five equations (three second order dynamic equations in f2 , f3 and v, and two constraint equations formed from a vector loop around the mechanism dotted with ıˆ and j ) also described the system. For brevity, the derivation
954
. . . .
Figure 7. The flexible coupler four-bar linkage with a floating frame.
of equations of motion is not shown here. The interested reader can find the details in reference [21]. AUTO was again used to construct resonance curves and locate bifurcations. In Figure 8 are shown two resonance curves that demonstrate the effect of non-dimensional crank length. In this example, the four-bar has non-dimensional crank length R1 ( = a1 /L) equalling 0·04 and 0·08, R3 ( = a3 /L) = 0·70, R4 ( = a4 /L) = 0·90, where a1 , L, a3 and a4 are the lengths of the crank, undeformed coupler, follower and
Figure 8. Resonance curves of mid-span rod deflections with instability zones, R1 = 0·04 and 0·08, R3 = 0·7, R4 = 0·9, M = 2.
955
Figure 9. Non-linear equation instability zones in the non-dimensional crank length (R1 ) versus non-dimensional crank speed (V) plane.
baseline, respectively. The ratio of the mass of the follower to the coupler (M) is 2·0, and the non-dimensional damping m1 = 0·01. The material and cross-section properties of the flexible coupler are the same as the connecting rod in the slider–crank. The non-dimensional crank speed V is the ratio of the crank speed to a flexible coupler’s first natural frequency of 38 hz. The non-dimensional response g is the coupler mid-span deflection divided by the coupler length. The follower is treated as a rigid and slender rod, the mass moment of inertia of which about its center of mass is 1/12 of the follower mass times a32. The system response g and its stability are totally described through the equation of motion by non-dimensional parameters R1 , R3 , R4 , M, m1 and V. As would be expected, the larger crank size induces a larger bending response in the coupler at a given non-dimensional crank speed (defined as the ratio of crank speed to first natural frequency of the pinned–pinned coupler in bending). The jump instability boundaries with a softening spring, and a period-doubling boundary are shown in Figure 8, and also are compiled in the R1–V parameter plane in Figure 9. In Figure 10 are shown three resonance curves that demonstrate the effect of follower length. R3 has values of 0·5, 0·7 and 0·964286, with R1 = 0·06, R4 = 0·90714 and M = 0·8143. The figure again shows a softening spring character of the jump, and a response that amplifies with followers of larger length. 5. CONCLUSIONS
An easily implemented method to derive equations of motion for flexible linkage mechanisms (such as constant crank speed planar slider–cranks and four-bars) is combined with a numerical scheme to study interesting non-linear effects, such as bifurcations and stability boundaries in the crank length–crank speed parameter plane. The method hinges on both a simple floating frame scheme and a generalized co-ordinate partitioning
956
. . . .
Figure 10. Resonance curves for R3 = 0·5, 0·7 and 0·964286; R1 = 0·06, R4 = 0·90714, M = 0·8143. ———, Stable period-1; r. . . . . r, unstable period-1, jump; q— – — –q, unstable period-1, period-doubling cascade.
algorithm that was specifically designed to interact with software (AUTO) that is used for continuation, bifurcation and construction of resonance curves. The floating frame is pinned on one end, and pinned-and-slotted on the other to permit flexible link length changes which are so important for stability studies. Geometric and dynamic non-linearities are easily included by having sufficient energies as input to the Lagrange equations. These energies are easily derived because of the ease with which they are described with respect to the floating frame. The technique does not require that the system be reduced to a minimum number of differential equations; that is the constraint equations and Lagrange multipliers need not be removed from the analysis. This fact considerably simplifies the determination of instabilities and improves and accuracy of stability boundaries by eliminating the need to apply physical reasoning, intuition and/or assumptions to remove complicating terms that stand in the way of reducing the number and size of equations to a tractable set. Stability boundaries and resonance curves are presented for a slider–crank with a flexible rod and four-bar mechanism with a flexible coupler. In each case a single mode is used to describe the response and stability. Bifurcations included jumps with a softening spring character and period-doubling cascades. In the case of the slider–crank, the results compared favorably with the experiments. Future work will extend this method to situations with more than one flexible link and more than one elastic mode per flexible link. REFERENCES 1. D. H and D. G. B 1995 Nonlinear Dynamics 7, 365–384. Experimental observations of a flexible slider crank mechanism at very high speeds. 2. D. G. B and S. W. L 1993 ASME Design Conference DE-Vol. 56, 123–130. Linear and nonlinear instabilities in a flexible mechanism. 3. D. G. B and S. W. L 1996 Mechanism and Machine Theory 31, 215–227. Nonlinear equation instability boundaries in flexible mechanisms.
957
4. S. W. L and D. G. B 1995 Zeitschrift fur Angewandte Mathematik und Mechanik 75, 883–888. Periodic solutions and stability of a nonlinear flexible slider crank mechanism by HBM/FFT. 5. S. H and S. W. S 1993 Nonlinear Dynamics 4, 573–603. The dynamic stability and nonlinear resonance of a flexible connecting rod: continuous parameter model. 6. S. H and S. W. S 1994 Journal of Sound and Vibration 170, 25–49. The dynamic stability and nonlinear resonance of a flexible connecting rod: single-mode model. 7. E. D and H. B. K 1991 International Journal of Bifurcation and Chaos 1, 493–520. Numerical analysis and control of bifurcation problems (I) bifurcation in finite dimensions. 8. B. V. V and R. S. A 1971 Transactions of the American Society of Mechanical Engineers, Journal of Engineering Industry 93, 251–262. Nonlinear dynamic response of elastic slider–crank mechanism. 9. M. B and W. K 1979 Transactions of the American Society of Mechanical Engineers, Journal of Mechanical Design 101, 149–153. Dynamic stability of elastic mechanisms. 10. H. A, Z. D. M and H. T 1990 Transactions of the American Society of Mechanical Engineers, Journal of Dynamic Systems, Measurement, and Control 112, 177–185. Inverse dynamics of flexible robot arms: modeling and computation for trajectory control. 11. D. G. B and S. W. L 1991 ASME Design Conference DE-Vol. 37, 161–166. Investigation of parametric resonance instability in a flexible rod slider crank mechanism. 12. O. P. A and A. A. S 1985 Computers and Structure 21, 1303–1312. Dynamic analysis of multibody systems using component modes. 13. A. A. S 1989 Dynamics of Multibody Systems. New York: John Wiley. 14. J. S. P 1968 Theory of Matrix Structural Analysis. New York: McGraw-Hill. 15. E. J. H 1989 Computer-aided Kinematics and Dynamics of Mechanical Systems. Boston, Massachussetts: Allyn and Bacon. 16. J. G and P. H 1983 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. 17. J. P. S and G. N. S 1974 Transactions of the American Society of Mechanical Engineers, Journal of Engineering for Industry 96, 411–419. Nonlinear vibration analysis of elastic four-bar linkages. 18. C. K. S, B. S. T, T. M. X and C. H. W 1986 Mechanism and Machine Theory 21, 121–133. An experimental study on the nonlinear elastodynamic response of linkage mechanisms. 19. M. R. S and L. M 1971 Journal of Mechanical Engineering Science 13, 237–242. Stability of a four-bar linkage with flexible coupler. 20. W. L. C, B. T and R. G. F 1984 Mechanism and Machine Theory 19, 307–317. Critical running speeds and stability of high-speed flexible mechanisms. 21. S. W. L 1994 Ph.D. Thesis, Auburn University. Elastic stability and active control of a flexible rod mechanism.