- Email: [email protected]
Perturbation method for the floquet eigenvalues and stability boundary of periodic linear systems
Journal of Sound and Vibration (1995) 182(2), 245–257
PERTURBATION METHOD FOR THE FLOQUET EIGENVALUES AND STABILITY BOUNDARY OF PERIODIC LINEAR SYSTEMS W.-T. W, J. H. G J. A. W Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, U.S.A. (Received 23 April 1993, and in final form 17 February 1994) A multiple parameter perturbation method is developed to determine the Floquet eigenvalues and stability boundary of a linear discrete system that is described by a system of ordinary differential equations with periodic coefficients. In the method, the state of the system is determined by solving a number of decoupled problems and, as a result, implementation of the method is computationally efficient even when the dimension of the system is large. The method is first applied to Mathieu’s equation in order to illustrate its application to a relatively straightforward, standard problem. It is then applied to the torsional analysis of a boat’s drive train to show application to a large scale system of some practical importance.
1. INTRODUCTION
The dynamic behavior of the class of systems that includes reciprocating mechanisms, shafts connected by universal joints, columns subjected to periodic loads, and control systems with periodic feedback is governed by linear ordinary differential equations with periodic coefficients. One important issue in analyzing such systems is their stability, which can be highly sensitive to changes in the amplitude and frequency of the periodic coefficients. A traditional approach in dealing with periodic systems is to identify stable and unstable regions in the amplitude–frequency parameter plane. A variety of methods is available for this purpose; see, for instance, reference [1]. In the present paper, a new method is developed for characterizing stability, and the method is especially useful in the study of large-dimensional systems where the frequency can vary over a wide range, and the amplitudes of the periodic terms are specified to be small. Existing stability methods can be categorized into three groups; Lyapunov-based, numerical and analytical. First, the methods using Lyapunov functions provide conservative bounds on the values of the parameters that yield stable solutions [2–5]. A practical limitation of these methods, however, is that it can be difficult to find a functional which provides tight stability criteria; consequently, the stability estimates are often overly conservative. Second, numerical methods typically solve for the transition matrix through integration of the governing equation. The eigenvalues of the transition matrix determine whether or not the magnitude of the solution increases over an excitation period [6, 7]. In some approaches [8], the dimension of the system can be reduced to simplify the evaluation of the transition matrix. One such method involves finding the derivatives of the transition matrix with respect to the parameters, and using this information directly to trace the 245 0022–460X/95/170245 + 13 $08.00/0
7 1995 Academic Press Limited
.-. .
246
stability boundary [9–11]. The key advantage of the numerical approach is that an accurate prediction of stability can be obtained even when the amplitudes of the periodic coefficients are large. The primary disadvantage, however, is that strictly numerical schemes can be computationally expensive when the dimension of the system is large or when the excitation frequency varies over a wide range. Third, various analytical methods are also available. The harmonic balance method [12, 13] can be combined with a perturbation expansion [14], or with boundary tracing [15, 16], to provide stability information. A combination of averaging and perturbation methods [17–21] has also been used. Another approach [22, 23] is based on perturbation expansions, in which the system’s state is assumed to be the product of periodic and exponential functions in accordance with Floquet theory. In reference [22], for instance, the Floquet eigenvalues were systematically expressed in terms of a perturbation expansion. Unperturbed eigenvalues were classified as being either ‘‘exactly’’ congruent or ‘‘not congruent’’ with other unperturbed eigenvalues. However, when large systems are considered, the existence of many ‘‘almost congruent’’ eigenvalues can cause the perturbation expansions to converge either poorly, or not at all [23]. In the present paper, a perturbation method is developed which takes into account the fact that some of the eigenvalues are almost congruent for large-dimensional systems. The representation used here for the system’s state is similar to that given in reference [22]. However, the periodic and exponential functions are each expressed in terms of perturbation expansions in three variables: the amplitude of the periodic coefficients, a frequency parameter, and a term that quantifies the closeness of the almost congruent eigenvalues. A criterion is given in section 2.2.2 to identify eigenvalues which are members of an almost congruent set. With this criterion, the state space can be divided into a collection of uncoupled subspaces, and perturbation expansions can be determined for the eigenvalues associated with each subspace. Thus, even for large systems, the method can be computationally efficient in that it deals with only a small subgroup of the eigenvalues at any given time.
2. CALCULATION OF FLOQUET EIGENVALUES
2.1. The system of interest is governed by [I + A (t˜ )]
dx = [L + B (t˜ )]x, dt˜
(1)
where x is an N-dimensional state vector, L is a constant N × N matrix, A and B matrices which are periodic functions of time t˜ , and >A > Q >I>. Here, >·> denotes infinity norm. Without loss of generality, L is specified to be diagonal. Because they specified to be periodic, the system matrices A and B can be expressed in terms of combined Fourier and power series P
A (t˜ ) = s
S
isvt˜ s o pA (p) , s e
p = 1 s = −S
P
B (t˜ ) = s
are the are the
S
isvt˜ s o pB (p) , s e
(2)
p = 1 s = −S
where i = z−1. The small parameter e is related to the amplitude of the parametric excitation; v is the fundamental frequency of the periodically varying terms; and P and S are specified integers. For simplicity of the presentation, equation (2) contains only the one parameter, e, but in principle this approach can be extended to multiple parameters as well.
247
In terms of the excitation frequency, consider the change of variable t = vt˜ , where v = v0 (1 + d), in equation (1). Here, v0 is the mid-point of the frequency range of interest, and d is the detuning parameter. For instance, with v1 Q v Q v2 , these become v0 = (v1 + v2 )/2,
=d = Q (v2 − v1 )/(2v0 ).
(3)
With this substitution, equation (1) is written
6
P
I+ s
S
ist s o PA(p) s e
p = 1 s = −S
7
6
76
a dx = 1 + s (−1)kd k dt k=1
P
L+ s
S
7
ist s o pB(p) x, s e
p = 1 s = −S
(4)
(p) where L = L /v0 , A(p) (p) (p) s =A s and Bs = B s /v0 . In this manner, the periodic coefficients of the governing differential equation are represented explicitly as polynomials in the small quantities d and e. There are N independent solutions of equation (4), and each can be expressed in the form
X = F(t) elt .
(5)
Here F is an N-vector that has the same period (2p) as the coefficients in equation (4), and l is a Floquet eigenvalue. Clearly, the stability of the system is determined by the real part of l. 2.2. In subsequently developing the perturbation solution below, the case of repeated eigenvalues is particularly important and requires special treatment. In effect, the l can be repeated if a pair of eigenvalues differs by any integer m multiple of i because from equation (5), terms such as eimt could be grouped with the periodic vector F(t), still having a period of 2p. When this situation occurs, the eigenvalues are said to be congruent. Here it is useful to classify eigenvalues of equation (4) by their unperturbed values, namely, when o = d = 0. If a particular unperturbed eigenvalue lj satisfies lj = mj i + l0 , for j = 1, 2,..., J, it is said to belong to an exactly congruent J is the number of eigenvalues in the set, and l0 is a specified the primary unperturbed eigenvalue. The minimum number of set is one and, consequently, every eigenvalue belongs to at In addition, eigenvalues that are almost congruent satisfy l = mi + l0 + g,
(6) set. The mj are integers, complex number termed elements in a congruent least one congruent set. (7)
where m is an integer, =g = Q K, and K is a specified small number. The asymptotic solution for the Floquent eigenvalues is presented in two steps. First, a procedure is developed to account for exactly congruent eigenvalues. It will be seen that if the individual eigenvalues are treated independently and expanded as polynomials in o and d, then the solvability condition in the perturbation solution cannot be satisfied if l is a member of a congruent set with more than one element. Consequently, a new representation for F and l is constructed in which the eigenvalues for all members of the congruent set are expanded simultaneously. This technique allows the solvability condition to be satisfied. In the second step, the radius of convergence for the solution in the o–d plane is enlarged by transforming the governing equation (4) into a form in which the almost congruent eigenvalues are rendered exactly congruent. The contribution of this work lies primarily in treating congruent, and almost congruent, eigenvalues as distinct subsets that interact when the equations are perturbed,
.-. .
248
and in increasing the radius of convergence of the representation. Congruent and almost congruent eigenvalues often occur when analyzing large-dimensional systems since they can have many eigenvalues which differ primarily in their imaginary parts. When this is the case, a number of the eigenvalues will probably be nearly congruent as implied by equation (7). Thus, the presented method is especially useful for dealing with these larger systems. 2.2.1. Exactly congruent eigenvalues Substitution of equation (5) in equation (4) requires that F and l satisfy
6
P
I+ s
7
S
ist s o pA(p) s e
p = 1 s = −S
6
7
P S dF ist Fl + I + s s o pA(p) s e dt p = 1 s = −S
6
P
=(1 − d + d 2 + · · · ) L + s Assuming that l and F can be expressed as a
a
a
l = s s o pd ql (p,q) ,
a
F(t) = s s
p=0 q=0
S
7
ist s o pB(p) F. s e
p = 1 s = −S
(8)
a
s o pd qFs(p,q) eist ,
(9)
p = 0 q = 0 s = −a
substitution of equation (9) in equation (8), and grouping of like terms, yield {(si + l (0,0) )I − L}Fs(p,q) = −Fs(0,0)l (p,q) + R (p,q) , where p e 0, q e 0 and −a E s E a. The residual term R by p
R (p,q) = − s
(p,q)
(10)
in equation (10) is given
S
(p − j,q) (p,q − 1) s (s − m)A(j) + LFs(p,q − 2) m Fs − m − LFs
j = 1 m = −S
p
−s
S
j = 1 m = −S
p
+s
p−j
q
(n,k) (p − j − n,q − k) s A(j) m s s F s − ml n=0 k=0
S
p
(p − j,q) s B(j) m Fs − m − s
j = 1 m = −S
p
+s
S
(p − j,q − 1) s B(j) m Fs − m
j = 1 m = −S
p−1 q−1
S
(p − j,q − 2) s B(j) + s s Fs(n,k)l (p − n,q − k) m Fs − m
j = 1 m = −S
n=1 k−1
q
q−1
k=1
k=0
+ s Fs(0,k)l (p,q − k) + s Fs(p,k)l (0,q − k) .
(11)
The set of equations described by equation (10) is solved recursively, with the zero order components of l (p,q) and Fs(p,q) being found first, followed by the higher order components. As an example, consider a perturbation expansion when l (0,0) is the nth unperturbed Floquet eigenvalue, i.e., the nth diagonal element of L. Then the Fs(0,0) are each zero vectors, except for F0(0,0) where the nth element alone equals 1. By these assignments, the term (0,0) F0(0,0) el t is precisely the nth solution of the unperturbed system in which the associated eigenvalue is l (0,0) . With the zero order terms so established, and substituted in equation (11), the higher order terms in the l and F expansions can be found recursively through equation (10). In solving for the higher orders terms in equation (10), the nth element of the diagonal matrix {(si + l (0,0) )I − L} vanishes when s = 0. Therefore, a solution for
249
F will exist if and only if the nth element of the vector {−F l + R } also vanishes, which implies that l (p,q) is equal to the nth component of R . This requirement is termed the solvability condition, and it provides the additional constraint which is needed to determine l (p,q) . Once l (p,q) is determined, Fs(p,q) can be found through equation (10). In this case, the form of the perturbation series is uniquely determined by setting the nth component of Fs(p,q) to zero. The above procedure fails when two (or more) eigenvalues of the unperturbed system differ by an integer multiple of i, suggesting that exactly congruent eigenvalues be treated together. For instance, consider the case in which the congruent set contains two eigenvalues. It is necessary to determine the solution for the corresponding modes simultaneously so that the two zero elements in {(si + l (0,0) )I − L} can be accommodated together. This can be done by grouping the eigenvalues and associated vectors in the exactly congruent set. To this end, the term F in equation (5) becomes an N × 2 matrix instead of an N-vector, and l becomes a 2 × 2 matrix. Similarly, equations (9)–(11) hold for the matrices Fs(p,q) and l (p,q) . In a manner analogous to the treatment of a single eigenvalue above, two elements of this matrix are used to accommodate the two zero elements in {(si + l (0,0) )I − L} when solving for one of the eigenvectors, and the other two elements are used when solving for the other eigenvector. The procedure can be extended to J exactly congruent modes. (p,q) 0
(0,0) (p,q) s (p,q)
(p,q)
2.2.2. Almost congruent eigenvalues When two eigenvalues are almost congruent, the norms of Fs(p,q) and l (p,q) substantially increase with perturbation order, weakening convergence of the perturbation expansion. Consequently, for given values of o and d, the summed series in equation (9) might not converge. Consider, for instance, the case in which an expansion is constructed for the nth eigenvalue l (0,0) , and the kth eigenvalue, l (0,0) + mi + g, is almost congruent. In turn, the kth element of {(−mi + l (0,0) )I − L} in equation (10) will be then approximately (p,q) zero, so that an upper bound of the norm of F−m becomes =g =−1>R (p,q) >. Since R (p,q) in equation (10) is a linear combination of the lower order terms of Fs(p,q) and l (p,q) , it is possible (p) that >Fs(p,q) > and =l (p,q) = each grow at the rate {=g =−1 maxp,s (>L>, >A(p) s >, >Bs >)} with respect to perturbation order. If this rate is large, a sufficient condition for convergence requires that (p) max (>A(p) s >, >Bs >)=o = E =g =, p,s,e
(p) max (>L>, >A(p) s >, >Bs >)=d = E =g = p,s,d
(12)
are each satisfied. To improve convergence, the almost congruent eigenvalues are included in the exactly congruent set in the following manner. First, the almost congruent eigenvalues that are included in the set satisfy =l − l* − mi= E K,
(13)
where (p) (p) (p) K = max {max (>A(p) s >, >Bs >)=o =, max (>L>, >As >, >Bs >)=d =}. p,s
p,s
Here, l is an unperturbed eigenvalue which does not belong to the congruent set, the l* are the unperturbed eigenvalues which belong to the congruent set, and again m is any integer. Note that if all eigenvalues are grouped in a common set, then the sufficient condition (12) does not guarantee that the expansion will converge. To include an almost congruent eigenvalue in the set, a new perturbation parameter b is introduced in order to
.-. .
250
redefine the zeroth order eigenvalues. In particular, equation (4) is rewritten as
6
P
I+ s
7
S
ist s o pA(p) s e
p = 1 s = −S
6
76
a dx = 1 − s (−1)kd k dt k=1
P
L0 + bLd + s
S
7
ist s o pB(p) x, s e
p = −P s = −S
(14) where bLd is a constant N × N diagonal matrix with the non-zero elements corresponding to the g values of the associated almost congruent eigenvalues, L0 is a new unperturbed system matrix, and L0 + bLd = L. In this transformation, the original unperturbed system matrix is decomposed through L = L0 + bLd , so as to render the almost congruent eigenvalues in the original unperturbed system exactly congruent in this new formulation. Owing to the term bLd , equation (14) contains the three parameters o, d and b. The J × J matrix Q and the N × J matrix F, defined from x = F(t) eQt , become a
a
a
Q = s s s o pd qb rQ(p,q,r) , p=0 q=0 r=0
a
a
a
F(t) = s s s
a
s o pd qb rF(p,q,r) eist , s
(15)
p = 0 q = 0 r = 0 s = −a
To solve for the terms Q(p,q,r) and F(p,q,r) , equations (11) and (15) are substituted in s equation (14). Grouping of the terms of like order yields the recursion equation {siI + Q(0,0,0) − L0 }F(p,q,r) = −F(0,0,0) Q(p,q,r) + R(p,q,r) , s s
(16)
(p,q,r)
where the N × J residual matrix R is the sum of all remaining terms. The zeroth order term Q(0,0,0) is assigned to be l0 I, where l0 is the primary unperturbed eigenvalue of the congruent set, and then the zeroth order terms F(0,0,0) are found from equation (14) with s (0,0,0) o = b = d = 0. By virtue of this assignment of Q(0,0,0) , the term sa F(0,0,0) eist eQ t s s = −a corresponds to the J solutions of the unperturbed system. F(p,q,r) and Q(p,q,r) are found s (p,q) (p,q) in a manner analogous to that discussed for Fs and l in equation (10). Once Q is determined, the J Floquet eigenvalues are found as the eigenvalues of the matrix Q. It should be noted that since the new perturbation parameter b has been introduced, the sufficient requirement (13) for convergence must also be revised. In the present method, the N modes of the system are divided into a number of smaller sets. Each set is examined separately, unlike a conventional numerical integration scheme. To that extent, the present method can be computationally more efficient than an entirely numerical method. Furthermore, since the present method expresses the eigenvalues as a polynomial, the eigenvalues for different values of the parameters can be obtained with little additional computational effort once the coefficient matrices are determined. 2.3. l0 Convergence of the perturbation expansion can be improved if the quantity >bLd > becomes smaller. Because >bLd > = maxa (=ga =), where the ga are deviations of the almost congruent eigenvalues, and since the ga depend on the choice of the primary unperturbed eigenvalue l0 , it is possible to find a l0 which will minimize >bLd >. A simple optimization procedure follows. First, suppose that the primary unperturbed eigenvalue is temporarily set to the first unperturbed eigenvalue of the set. All the elements of this set can be expressed as la = l 0 + m˜a i + g˜ a = l0 + ma i + ga ,
(17)
where m˜a = ma + m; l0 is the primary unperturbed eigenvalue; and ma and m are integers, with a = 1,..., J. The tilde superscript indicates current values, and the terms without superscripts indicate optimized values. Since >bLd > = maxa (=ga =), the optimum value of l0
251
and the corresponding values of m are found by minimize {maximum (=ga =)}, l 0 ,m
where =ga = = =l 0 + mi + g˜ a − l0 =.
(18)
a
3. STABILITY BOUNDARY
Since the real parts of the Floquet eigenvalues determine stability, the perturbation expansion can be used to establish the transition between stable and unstable behavior in the o–v and o–d planes. Note that the perturbation expansion for Q needs to be calculated only once to be useful over a range of o and v values. The procedure for the stability boundary follows. First, a rectangular, evenly spaced mesh on the parameter plane is generated. Second, the mid-point of the frequency range is used as the frequency of the unperturbed system. Third, the coefficients in the perturbation expansions of F and Q are found. Finally, the Floquet eigenvalues are determined at each mesh point. If the real part of any of the Floquet eigenvalues is greater than zero, then the grid point corresponds to an unstable response. A practical limitation of the method is that if the specified range of v is too wide, then the maximum value of d will also be large and all unperturbed eigenvalues will be included in a congruent set. As a result, the expansion of F and Q may not converge. To avoid this difficulty, the range of v is divided into a number of subdivisions so that the maximum value of d in each division satisfies (p) (p) (p) (19) =dmax = Q =omax =max (>B(p) s >, >As >)/max (>L>, >Bs >, >As >). p,s
p,s
Because of this restriction, the influence of dmax on deciding the value of K in equation (13) is effectively eliminated. 4. EXAMPLES
4.1. The Mathieu equation dx d2x + 2z ˜ + (1 + o cos vt˜ )x = 0 dt dt˜ 2
(20)
is a classical model of parametric excitation, where o characterizes the strength of the excitation, and z is a damping coefficient. Equation (20) is rewritten in the first order form dx = dt˜
0$
i 0
% $
0 −1 +z −i i
% > 0$
−i +o 4 −1
%
$
% 11
1 i e−ivt˜ + −i 1
i 1
1 eivt˜ −i
x.
(21)
The system experiences a primary parametric resonance when v = 2 and, consequently, stability will be investigated for a range of frequencies centered about this critical value. Thus, let d = (v − 2)/2 and t = vt˜ . Equation (21) becomes dx (1 − d + d 2 ) = 2 dt
6$
i 0
% $
0 −1 +z −i 1
% $
−i o i + −i 4 1
%
7
1 (eit + e−it ) x, −i
(22)
where the higher order terms in d have been omitted. In equation (22), if z = 0, the system has the two exactly congruent unperturbed eigenvalues 0·5i and −0·5i. Let Ad = 12
$
−1 i
%
−i −1
and
b = z,
.-. .
252
to employ the approach outlined in section 2.2.2. Then, F and Q for this example will each be 2 × 2 matrices, and each can be expanded: Q = 12
F=
6$
% $
−i 0 −1 +z 0 −i 0
6$ % $ % 7 0 1
+
z 2
0 0 + 0 0
1 it o e + 0 16
6$ % $ % 7 1 0
0 0 + 0 0
% $
0 i +d −1 0
0 it e , 1
% $ %7
0 o 0 + −i 4 1
6$ % $ i 2
0 −it 0 e + 0 0
1 0
% $
,
−2 0 + 0 −2
(23a)
% $
0 it 0 e + 0 0
% 7
2 e2it −i
(23b)
to first order as in equation (15). After restoring the time scale, the Floquet eigenvalues of matrix Q become l = (−i − z 2 z(0·25o)2 − d 2 )(1 + d). (24) In Figure 1, values of the real part of the eigenvalue are plotted as functions of o for a range of v values. The agreement among the result (24) and the eigenvalues calculated numerically from the transition matrix is typically good. For instance, when v = 2·3 (d = 0·15), the difference between the two solutions is within 5%, even though O(o) = 1. Because of the simple form of equation (24), one may directly determine the equation for the stability boundary by requesting Re(l) = 0. The line demarcating the stability zone is given by o 2 = 16(d 2 + z 2 ). (25) In Figure 2 are shown the stability boundaries (25) compared with the numerical results for four different values of the damping coefficient. The error in the perturbation method increases with o, d and z, but in general the agreement between the two methods is good. The damped Mathieu equation is treated by Nayfeh and Mook [24]. There, the damping constant is proportional to o, and the Lindstedt–Poincare´ technique is used to develop a
Figure 1. The real part of the Floquet eigenvalues in Example 1 for v = 2·0, 2·1, 2·2, 2·3 and 2·4: W, numerical integration, ——, perturbation.
883844—JSV 182/2
253
Figure 2. Stability boundaries in Example 1 for z = 0·0, 0·05, 0·1 and 0·2: W, numerical integration; ——, perturbation.
perturbation expansion for the solution. This approach is appropriate if the solution is periodic; namely, o and d lie on the stability boundary. The solution in reference [24] for the stability boundary is, except for differing notation, identical with equation (25). Lion [22] uses a perturbation approach which is similar to that presented here in that it treats a first order system. The perturbation solution for the undamped Mathieu equation, however, does not take into account the existence of congruent eigenvalues. In terms of the current notation, reference [22] predicts the Floquet eigenvalues l = 2i(2/v + o 2/2v(v 2 − 4)),
(26)
and the stability boundary is found by using the criterion that on the stability boundary the solution has the same periodicity as the coefficients. Together with v = 2(1 + d), this implies that the stability boundary is given by o 2 = 2(v − 2)2(v + 2) 1 32d 2.
(27)
In contrast, equation (25) gives the result that o 2 1 16d 2 when z = 0. From Figure 2, the present method yields the correct asymptotic behavior when z vanishes and d is small, and equation (27) is proportionally off by a multiplicative factor of z2. In addition, equation (26) predicts that only the imaginary part of the eigenvalue is affected by the perturbation o. The results depicted in Figure 1 indicate that this is, in fact, not the case and that the real part of the eigenvalue also changes. It should be noted that Lion’s method was not intended for the case of congruent eigenvalues. This is the reason that equation (26) becomes singular as v = 2, and the method does not accurately predict the behavior of the eigenvalues near primary resonance. However, the behavior near resonant points is critical, since this is the region where vibration amplitudes can be the greatest. 4.2. The trend in drive train design for marine applications is towards lightweight, high performance systems. A typical drive train consists of an engine, flywheel, gearbox, shaft and propeller. The power to weight ratio may be increased, for example, by increasing the size of the pistons and by reducing the mass of the flywheel. Since the inertia of the piston changes periodically with time, this design trend tends to increase the ratio of the variable
254
.-. .
inertia to the constant inertia terms in the equation of motion; as a result, instability from parametric excitation becomes a concern. In the general formulation presented here, it is assumed that the inertia, damping and stiffness terms may vary with the angular position of the shaft. In fact, both the stiffness and inertia terms vary periodically if a U-joint is used to offset part of the drive train. Denoting the rotational speed of the drive train by v, and linearizing the governing equations for small amplitude vibration U(t˜ ), the equation of motion becomes [25] (M + oM (t˜ ))U (t˜ ) + (C + voC (t˜ ))U (t˜ ) + (K + v 2oK (t˜ ))U(t˜ ) = E(t˜ ) + v 2oE (t˜ ).
(28)
The parameter o is small and indicates the relative magnitude of the periodically varying terms. In the case of pistons, o = r/(2 + r), where r is the ratio of the varying moment of inertia of the connecting rod and piston to its average value. It is assumed that the matrices M , C , K and the vectors E and E are periodic with a period corresponding to one cycle of the engine, namely 2p/v. The homogeneous form of equation (28) can be transformed into the form of equation (1): [I + oA 1 (t˜ )]x˙(t˜ ) = [L + oB 1 (t˜ ) + ovB 2 (t˜ ) + ov 2(t˜ )B 3 (t˜ )]x(t˜ ),
(29)
with x(t˜ ) = [U(t˜ ), U (t˜ )]T. A 1 , B 1 , B 2 and B 3 are periodic and the elements of the diagonal matrix L are the unperturbed eigenvalues of the system. From equations (3) and (13), a particular perturbation expansion to determine the Floquet eigenvalues holds only over a limited frequency range. A practical consideration which becomes apparent when dealing with large systems is the best choice of the frequency range in the expansion. For example, consider a boat drive with motion characterized by 11 rotational degrees of freedom, and parameters o E 0·5 and 2700 E v E 3000. One approach would be to try to find a single perturbation expansion that would cover the entire frequency range from 2700 to 3000 rad/s. If this approach is used, then from equation (3), v0 = 2850 and =d= E 0·053. Inspection of the matrices involved indicates that =d= determines the value in equation (13) of K = 0·158 and, as a result, all 22 eigenvalues must be grouped together as a congruent set. This adversely affects convergence, since all eigenvalues belong to one congruent set. Instead, if the frequency range is divided into three equal sub-intervals then the term involving =o= in equation (13) dominates, and K is reduced to 0·096; with this choice, the maximum number of eigenvalues that must be grouped together in any of the frequency intervals is eight. As a result, the term b >Ad > is reasonably small and the convergence rate of the perturbation expansion is enhanced. From a practical point of view, the procedure discussed in section 3 is always utilized to restrict the frequency range in accordance with equation (19) so that d is sufficiently small so that o determines K. In this example, the Floquent eigenvalues of the system were calculated for the boat drive train using both the perturbation method, and the numerical transition matrix approach for 2000 E v E 4000 and o E 0·5. The transition matrices were calculated using a fourth order Runge–Kutta algorithm to numerically integrate the equations of motion over one period of excitation. The results are shown in Figures 3 and 4. Both methods indicate similar instability regions. As might be expected, the perturbation method accurately predicts the stability boundary for smaller values of o, but it does not do as well when o is large. The perturbation method provides different expansions in each different frequency range and, as a result, it erroneously predicts that the stability regions change abruptly from one frequency range to the next. The magnitude of the frequency range used in this example was determined by the long computation times for the transition matrices if a lower, more realistic frequency
255
Figure 3. Unstable regions (shaded) in Example 2 for 2000 Q v Q 4000 (perturbation).
range was used. This difficulty occurs because the time step in the Runge–Kutta algorithm must be sufficiently small so that the contribution of the mode with the highest natural frequency can be accurately integrated over the excitation period. Consequently, numerical integration times would increase by almost an order of magnitude if v were decreased to the actual operating range of the drive train, say 50–300 rad/s. In contrast, computer time was not a factor in calculating the perturbation solutions, which took approximately a factor of 200 less time on the computer than the Runge–Kutta method. Furthermore, the time required in calculating the perturbation solution decreases when the system’s stability is analyzed for the actual, lower frequency range. A final observation concerns boundary tracing techniques for large-dimensional systems. When the number of eigenvalues is large, there is a significant number of para-
Figure 4. Unstable regions (shaded) in Example 2 for 2000 Q v Q 4000 (numerical integration).
256
.-. .
metric resonances that can occur in a given frequency range. As a result, there is an abrupt hange in the slope of the stability boundary when the different resonance regions overlap. Furthermore, a number of higher order resonances can result in very narrow stability regions, rendering the boundary tracing technique difficult to implement.
5. SUMMARY
A new perturbation method has been developed for calculating the Floquet eigenvalues and eigenvectors of a system of ordinary differential equations with coefficients that vary periodically in time. The method is useful in dealing with large systems because it directly treats the cases of almost congruent and exactly congruent eigenvalues. This can be an important consideration to the extent that the likelihood of eigenvalues being nearly congruent increases when there are more eigenvalues in the system. The method expresses the solution in terms of a multi-parameter expansion where the parameters typically correspond to a parameter which characterizes the magnitude of the periodic coefficients, a parameter which characterizes the frequency variation, and a third parameter which reflects the closeness of the nearly congruent eigenvalues. The third parameter is selected so as to make all of the nearly congruent eigenvalues exactly congruent and, then, a perturbation expansion is applied which is appropriate for subsets of exactly congruent eigenvalues. It is found that by grouping the nearly congruent eigenvalues together in this manner, the radius of convergence of the expansion is significantly improved. Subsets of exactly congruent eigenvalues should be grouped together. If there are M elements in the subset, then the eigenvalues can be expressed as an expansion in which the coefficients are M × M matrices. Once the coefficient matrices in the perturbation expansions are determined, the eigenvalues can easily be calculated for different values of the parameters. This result in an efficient computational method for larger systems. The method is applied to the damped Mathieu equation, and to torsional drive train vibration. In the case of the damped Mathieu equation, the first order solution gives a single analytical expression for the eigenvalues in terms of the amplitude of the excitation, the frequency and the amount of viscous damping in the system. The results agree reasonably well with those of independently generated numerical solutions. When the real part of the eigenvalue vanishes, the perturbation expansion predicts an equation for the stability boundary which is consistent with other perturbation methods and which also compares favorably with the numerical results. In the case of the boat drive train, the shape of the stable region is much more complicated because of complex interactions that occur when there is a large number of eigenvalues. The perturbation method predicts a stability region which is very similar to that determined numerically. As might be expected, the agreement is very good for small values of the perturbation parameters and worsens as the magnitudes of the perturbation parameters increase. The perturbation method is several orders of magnitude more computationally efficient than the numerical method.
ACKNOWLEDGMENT
The authors acknowledge partial support for this work from the United States Navy under contract N00189-90-C-0331. Mr James Brown is the contract monitor.
257
REFERENCES 1. V. A. Y and V. M. S 1975 Linear Differential Equations with Periodic Coefficients. New York: John Wiley. 2. S. K. S 1981 Journal of Applied Mechanics 48, 174–177. Stability theorems for multidimensional linear systems with variable parameters. 3. S. Z. A-R and G. A 1986 International Journal of Systems Science 17, 1654–1660. Stability analysis of non-autonomous linear systems by a matrix decomposition method. 4. R. K. Y and S. R. K 1988 International Journal of Systems Science 19, 1853–1858. Stability analysis of linear time-varying systems. 5. M. V. M 1990 Avtomatika i Telemekhanika 4, 27–35. An algorithm of analysis of the stability of linear periodic systems and its computer realization. 6. C. S. H 1974 Journal of Mathematical Analysis and Applications 45, 234–251. On approximating a general linear periodic system. 7. P. F, C. E. H and T.-H. W 1977 International Journal for Numerical Methods in Engineering 11, 1117–1136. Efficient numerical treatment of periodic systems with application to stability problems. 8. P. M. G, P. B, G. F and M. L 1987 NASA Conference; 58th Shock and Vibration Symposium 1, 487–494. Dynamic behavior of dissymmetric rotor bearings modelled with a periodic coefficient large system. 9. B. W and H. K 1991 Mechanism and Machine Theory 26(2), 123–131. Direct Floquet method for stability limits determination—I: theory. 10. H. K and B. W 1991 Mechanism and Machine Theory 26(2), 133–144. Direct Floquet method for stability limits determination—II: application and phenomena. 11. J. W. L and I. C 1990 American Institute of Aeronautics and Astronautics Journal 28, 1089–1097. Stability sensitivity analysis of a helicopter rotor. 12. A. U and M. A. B 1981 Journal of Applied Mechanics 48, 948–958. Parametric instability of a rotating shaft due to pulsating torque. 13. I. G. T and C. J. Y 1985 International Journal of Mechanical Science 27, 813–822. Effect of flexibility on the dynamic stability of a four-bar mechanism. 14. D. R. F, J. and H. T. J. L 1988 20th Biennial Mechanisms Conference: Trends and Developments in Mechanisms, Machines, and Robotics, 397–405. Stability analysis by a perturbation method. 15. S. L. L, Y. K. C and S. Y. W 1982 Journal of Applied Mechanics 49, 849–853. A variable parameter incrementation method for dynamic instability of linear and nonlinear elastic systems. 16. B. W 1988 Proceedings of the 2nd International Symposium on Transport Phenomena, Dynamics and Design of Rotating Machinery, Honolulu, 303–314. Rapid stability analysis of parameter-excited large systems. 17. C. S. H 1963 Journal of Applied Mechanics 30, 367–372. On the parametric excitation of a dynamic system having multiple degrees of freedom. 18. C. S. H 1965 Journal of Applied Mechanics 32, 373–377. Further results on parametric excitation of a dynamic system. 19. T. Y and A. S 1968 Bulletin of the Japan Society of Mechanical Engineers 11(43), 92–101. On the oscillations of ‘‘summed and differential types’’ under parametric excitation (vibratory systems with damping). 20. S. T. A and N. S N 1986 Journal of Sound and Vibration 107, 215–230. Dynamic stability of pipes conveying pulsating fluid. 21. S. T. A and S. F. A 1987 Journal of Mechanisms, Transmissions, and Automation in Design 109, 412–418. Dynamic stability of chain drives. 22. P. M. L 1966 Journal of The Franklin Institute 281(1), 27–40. Stability of linear periodic systems. 23. Y. N and P. Y. W 1969 Journal of The Franklin Institute 287(2), 143–157. A method for stability investigation of a periodic dynamic system with many degrees of freedom. 24. A. H. N and D. T. M 1979 Nonlinear Oscillations. New York: John Wiley. 25. W.-T. W 1993 Ph.D. Thesis, Department of Mechanical Engineering, Carnegie Mellon University. Stability analysis and forced response of reciprocating mechanisms.
PERTURBATION METHOD FOR THE FLOQUET EIGENVALUES AND STABILITY BOUNDARY OF PERIODIC LINEAR SYSTEMS W.-T. W, J. H. G J. A. W Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, U.S.A. (Received 23 April 1993, and in final form 17 February 1994) A multiple parameter perturbation method is developed to determine the Floquet eigenvalues and stability boundary of a linear discrete system that is described by a system of ordinary differential equations with periodic coefficients. In the method, the state of the system is determined by solving a number of decoupled problems and, as a result, implementation of the method is computationally efficient even when the dimension of the system is large. The method is first applied to Mathieu’s equation in order to illustrate its application to a relatively straightforward, standard problem. It is then applied to the torsional analysis of a boat’s drive train to show application to a large scale system of some practical importance.
1. INTRODUCTION
The dynamic behavior of the class of systems that includes reciprocating mechanisms, shafts connected by universal joints, columns subjected to periodic loads, and control systems with periodic feedback is governed by linear ordinary differential equations with periodic coefficients. One important issue in analyzing such systems is their stability, which can be highly sensitive to changes in the amplitude and frequency of the periodic coefficients. A traditional approach in dealing with periodic systems is to identify stable and unstable regions in the amplitude–frequency parameter plane. A variety of methods is available for this purpose; see, for instance, reference [1]. In the present paper, a new method is developed for characterizing stability, and the method is especially useful in the study of large-dimensional systems where the frequency can vary over a wide range, and the amplitudes of the periodic terms are specified to be small. Existing stability methods can be categorized into three groups; Lyapunov-based, numerical and analytical. First, the methods using Lyapunov functions provide conservative bounds on the values of the parameters that yield stable solutions [2–5]. A practical limitation of these methods, however, is that it can be difficult to find a functional which provides tight stability criteria; consequently, the stability estimates are often overly conservative. Second, numerical methods typically solve for the transition matrix through integration of the governing equation. The eigenvalues of the transition matrix determine whether or not the magnitude of the solution increases over an excitation period [6, 7]. In some approaches [8], the dimension of the system can be reduced to simplify the evaluation of the transition matrix. One such method involves finding the derivatives of the transition matrix with respect to the parameters, and using this information directly to trace the 245 0022–460X/95/170245 + 13 $08.00/0
7 1995 Academic Press Limited
.-. .
246
stability boundary [9–11]. The key advantage of the numerical approach is that an accurate prediction of stability can be obtained even when the amplitudes of the periodic coefficients are large. The primary disadvantage, however, is that strictly numerical schemes can be computationally expensive when the dimension of the system is large or when the excitation frequency varies over a wide range. Third, various analytical methods are also available. The harmonic balance method [12, 13] can be combined with a perturbation expansion [14], or with boundary tracing [15, 16], to provide stability information. A combination of averaging and perturbation methods [17–21] has also been used. Another approach [22, 23] is based on perturbation expansions, in which the system’s state is assumed to be the product of periodic and exponential functions in accordance with Floquet theory. In reference [22], for instance, the Floquet eigenvalues were systematically expressed in terms of a perturbation expansion. Unperturbed eigenvalues were classified as being either ‘‘exactly’’ congruent or ‘‘not congruent’’ with other unperturbed eigenvalues. However, when large systems are considered, the existence of many ‘‘almost congruent’’ eigenvalues can cause the perturbation expansions to converge either poorly, or not at all [23]. In the present paper, a perturbation method is developed which takes into account the fact that some of the eigenvalues are almost congruent for large-dimensional systems. The representation used here for the system’s state is similar to that given in reference [22]. However, the periodic and exponential functions are each expressed in terms of perturbation expansions in three variables: the amplitude of the periodic coefficients, a frequency parameter, and a term that quantifies the closeness of the almost congruent eigenvalues. A criterion is given in section 2.2.2 to identify eigenvalues which are members of an almost congruent set. With this criterion, the state space can be divided into a collection of uncoupled subspaces, and perturbation expansions can be determined for the eigenvalues associated with each subspace. Thus, even for large systems, the method can be computationally efficient in that it deals with only a small subgroup of the eigenvalues at any given time.
2. CALCULATION OF FLOQUET EIGENVALUES
2.1. The system of interest is governed by [I + A (t˜ )]
dx = [L + B (t˜ )]x, dt˜
(1)
where x is an N-dimensional state vector, L is a constant N × N matrix, A and B matrices which are periodic functions of time t˜ , and >A > Q >I>. Here, >·> denotes infinity norm. Without loss of generality, L is specified to be diagonal. Because they specified to be periodic, the system matrices A and B can be expressed in terms of combined Fourier and power series P
A (t˜ ) = s
S
isvt˜ s o pA (p) , s e
p = 1 s = −S
P
B (t˜ ) = s
are the are the
S
isvt˜ s o pB (p) , s e
(2)
p = 1 s = −S
where i = z−1. The small parameter e is related to the amplitude of the parametric excitation; v is the fundamental frequency of the periodically varying terms; and P and S are specified integers. For simplicity of the presentation, equation (2) contains only the one parameter, e, but in principle this approach can be extended to multiple parameters as well.
247
In terms of the excitation frequency, consider the change of variable t = vt˜ , where v = v0 (1 + d), in equation (1). Here, v0 is the mid-point of the frequency range of interest, and d is the detuning parameter. For instance, with v1 Q v Q v2 , these become v0 = (v1 + v2 )/2,
=d = Q (v2 − v1 )/(2v0 ).
(3)
With this substitution, equation (1) is written
6
P
I+ s
S
ist s o PA(p) s e
p = 1 s = −S
7
6
76
a dx = 1 + s (−1)kd k dt k=1
P
L+ s
S
7
ist s o pB(p) x, s e
p = 1 s = −S
(4)
(p) where L = L /v0 , A(p) (p) (p) s =A s and Bs = B s /v0 . In this manner, the periodic coefficients of the governing differential equation are represented explicitly as polynomials in the small quantities d and e. There are N independent solutions of equation (4), and each can be expressed in the form
X = F(t) elt .
(5)
Here F is an N-vector that has the same period (2p) as the coefficients in equation (4), and l is a Floquet eigenvalue. Clearly, the stability of the system is determined by the real part of l. 2.2. In subsequently developing the perturbation solution below, the case of repeated eigenvalues is particularly important and requires special treatment. In effect, the l can be repeated if a pair of eigenvalues differs by any integer m multiple of i because from equation (5), terms such as eimt could be grouped with the periodic vector F(t), still having a period of 2p. When this situation occurs, the eigenvalues are said to be congruent. Here it is useful to classify eigenvalues of equation (4) by their unperturbed values, namely, when o = d = 0. If a particular unperturbed eigenvalue lj satisfies lj = mj i + l0 , for j = 1, 2,..., J, it is said to belong to an exactly congruent J is the number of eigenvalues in the set, and l0 is a specified the primary unperturbed eigenvalue. The minimum number of set is one and, consequently, every eigenvalue belongs to at In addition, eigenvalues that are almost congruent satisfy l = mi + l0 + g,
(6) set. The mj are integers, complex number termed elements in a congruent least one congruent set. (7)
where m is an integer, =g = Q K, and K is a specified small number. The asymptotic solution for the Floquent eigenvalues is presented in two steps. First, a procedure is developed to account for exactly congruent eigenvalues. It will be seen that if the individual eigenvalues are treated independently and expanded as polynomials in o and d, then the solvability condition in the perturbation solution cannot be satisfied if l is a member of a congruent set with more than one element. Consequently, a new representation for F and l is constructed in which the eigenvalues for all members of the congruent set are expanded simultaneously. This technique allows the solvability condition to be satisfied. In the second step, the radius of convergence for the solution in the o–d plane is enlarged by transforming the governing equation (4) into a form in which the almost congruent eigenvalues are rendered exactly congruent. The contribution of this work lies primarily in treating congruent, and almost congruent, eigenvalues as distinct subsets that interact when the equations are perturbed,
.-. .
248
and in increasing the radius of convergence of the representation. Congruent and almost congruent eigenvalues often occur when analyzing large-dimensional systems since they can have many eigenvalues which differ primarily in their imaginary parts. When this is the case, a number of the eigenvalues will probably be nearly congruent as implied by equation (7). Thus, the presented method is especially useful for dealing with these larger systems. 2.2.1. Exactly congruent eigenvalues Substitution of equation (5) in equation (4) requires that F and l satisfy
6
P
I+ s
7
S
ist s o pA(p) s e
p = 1 s = −S
6
7
P S dF ist Fl + I + s s o pA(p) s e dt p = 1 s = −S
6
P
=(1 − d + d 2 + · · · ) L + s Assuming that l and F can be expressed as a
a
a
l = s s o pd ql (p,q) ,
a
F(t) = s s
p=0 q=0
S
7
ist s o pB(p) F. s e
p = 1 s = −S
(8)
a
s o pd qFs(p,q) eist ,
(9)
p = 0 q = 0 s = −a
substitution of equation (9) in equation (8), and grouping of like terms, yield {(si + l (0,0) )I − L}Fs(p,q) = −Fs(0,0)l (p,q) + R (p,q) , where p e 0, q e 0 and −a E s E a. The residual term R by p
R (p,q) = − s
(p,q)
(10)
in equation (10) is given
S
(p − j,q) (p,q − 1) s (s − m)A(j) + LFs(p,q − 2) m Fs − m − LFs
j = 1 m = −S
p
−s
S
j = 1 m = −S
p
+s
p−j
q
(n,k) (p − j − n,q − k) s A(j) m s s F s − ml n=0 k=0
S
p
(p − j,q) s B(j) m Fs − m − s
j = 1 m = −S
p
+s
S
(p − j,q − 1) s B(j) m Fs − m
j = 1 m = −S
p−1 q−1
S
(p − j,q − 2) s B(j) + s s Fs(n,k)l (p − n,q − k) m Fs − m
j = 1 m = −S
n=1 k−1
q
q−1
k=1
k=0
+ s Fs(0,k)l (p,q − k) + s Fs(p,k)l (0,q − k) .
(11)
The set of equations described by equation (10) is solved recursively, with the zero order components of l (p,q) and Fs(p,q) being found first, followed by the higher order components. As an example, consider a perturbation expansion when l (0,0) is the nth unperturbed Floquet eigenvalue, i.e., the nth diagonal element of L. Then the Fs(0,0) are each zero vectors, except for F0(0,0) where the nth element alone equals 1. By these assignments, the term (0,0) F0(0,0) el t is precisely the nth solution of the unperturbed system in which the associated eigenvalue is l (0,0) . With the zero order terms so established, and substituted in equation (11), the higher order terms in the l and F expansions can be found recursively through equation (10). In solving for the higher orders terms in equation (10), the nth element of the diagonal matrix {(si + l (0,0) )I − L} vanishes when s = 0. Therefore, a solution for
249
F will exist if and only if the nth element of the vector {−F l + R } also vanishes, which implies that l (p,q) is equal to the nth component of R . This requirement is termed the solvability condition, and it provides the additional constraint which is needed to determine l (p,q) . Once l (p,q) is determined, Fs(p,q) can be found through equation (10). In this case, the form of the perturbation series is uniquely determined by setting the nth component of Fs(p,q) to zero. The above procedure fails when two (or more) eigenvalues of the unperturbed system differ by an integer multiple of i, suggesting that exactly congruent eigenvalues be treated together. For instance, consider the case in which the congruent set contains two eigenvalues. It is necessary to determine the solution for the corresponding modes simultaneously so that the two zero elements in {(si + l (0,0) )I − L} can be accommodated together. This can be done by grouping the eigenvalues and associated vectors in the exactly congruent set. To this end, the term F in equation (5) becomes an N × 2 matrix instead of an N-vector, and l becomes a 2 × 2 matrix. Similarly, equations (9)–(11) hold for the matrices Fs(p,q) and l (p,q) . In a manner analogous to the treatment of a single eigenvalue above, two elements of this matrix are used to accommodate the two zero elements in {(si + l (0,0) )I − L} when solving for one of the eigenvectors, and the other two elements are used when solving for the other eigenvector. The procedure can be extended to J exactly congruent modes. (p,q) 0
(0,0) (p,q) s (p,q)
(p,q)
2.2.2. Almost congruent eigenvalues When two eigenvalues are almost congruent, the norms of Fs(p,q) and l (p,q) substantially increase with perturbation order, weakening convergence of the perturbation expansion. Consequently, for given values of o and d, the summed series in equation (9) might not converge. Consider, for instance, the case in which an expansion is constructed for the nth eigenvalue l (0,0) , and the kth eigenvalue, l (0,0) + mi + g, is almost congruent. In turn, the kth element of {(−mi + l (0,0) )I − L} in equation (10) will be then approximately (p,q) zero, so that an upper bound of the norm of F−m becomes =g =−1>R (p,q) >. Since R (p,q) in equation (10) is a linear combination of the lower order terms of Fs(p,q) and l (p,q) , it is possible (p) that >Fs(p,q) > and =l (p,q) = each grow at the rate {=g =−1 maxp,s (>L>, >A(p) s >, >Bs >)} with respect to perturbation order. If this rate is large, a sufficient condition for convergence requires that (p) max (>A(p) s >, >Bs >)=o = E =g =, p,s,e
(p) max (>L>, >A(p) s >, >Bs >)=d = E =g = p,s,d
(12)
are each satisfied. To improve convergence, the almost congruent eigenvalues are included in the exactly congruent set in the following manner. First, the almost congruent eigenvalues that are included in the set satisfy =l − l* − mi= E K,
(13)
where (p) (p) (p) K = max {max (>A(p) s >, >Bs >)=o =, max (>L>, >As >, >Bs >)=d =}. p,s
p,s
Here, l is an unperturbed eigenvalue which does not belong to the congruent set, the l* are the unperturbed eigenvalues which belong to the congruent set, and again m is any integer. Note that if all eigenvalues are grouped in a common set, then the sufficient condition (12) does not guarantee that the expansion will converge. To include an almost congruent eigenvalue in the set, a new perturbation parameter b is introduced in order to
.-. .
250
redefine the zeroth order eigenvalues. In particular, equation (4) is rewritten as
6
P
I+ s
7
S
ist s o pA(p) s e
p = 1 s = −S
6
76
a dx = 1 − s (−1)kd k dt k=1
P
L0 + bLd + s
S
7
ist s o pB(p) x, s e
p = −P s = −S
(14) where bLd is a constant N × N diagonal matrix with the non-zero elements corresponding to the g values of the associated almost congruent eigenvalues, L0 is a new unperturbed system matrix, and L0 + bLd = L. In this transformation, the original unperturbed system matrix is decomposed through L = L0 + bLd , so as to render the almost congruent eigenvalues in the original unperturbed system exactly congruent in this new formulation. Owing to the term bLd , equation (14) contains the three parameters o, d and b. The J × J matrix Q and the N × J matrix F, defined from x = F(t) eQt , become a
a
a
Q = s s s o pd qb rQ(p,q,r) , p=0 q=0 r=0
a
a
a
F(t) = s s s
a
s o pd qb rF(p,q,r) eist , s
(15)
p = 0 q = 0 r = 0 s = −a
To solve for the terms Q(p,q,r) and F(p,q,r) , equations (11) and (15) are substituted in s equation (14). Grouping of the terms of like order yields the recursion equation {siI + Q(0,0,0) − L0 }F(p,q,r) = −F(0,0,0) Q(p,q,r) + R(p,q,r) , s s
(16)
(p,q,r)
where the N × J residual matrix R is the sum of all remaining terms. The zeroth order term Q(0,0,0) is assigned to be l0 I, where l0 is the primary unperturbed eigenvalue of the congruent set, and then the zeroth order terms F(0,0,0) are found from equation (14) with s (0,0,0) o = b = d = 0. By virtue of this assignment of Q(0,0,0) , the term sa F(0,0,0) eist eQ t s s = −a corresponds to the J solutions of the unperturbed system. F(p,q,r) and Q(p,q,r) are found s (p,q) (p,q) in a manner analogous to that discussed for Fs and l in equation (10). Once Q is determined, the J Floquet eigenvalues are found as the eigenvalues of the matrix Q. It should be noted that since the new perturbation parameter b has been introduced, the sufficient requirement (13) for convergence must also be revised. In the present method, the N modes of the system are divided into a number of smaller sets. Each set is examined separately, unlike a conventional numerical integration scheme. To that extent, the present method can be computationally more efficient than an entirely numerical method. Furthermore, since the present method expresses the eigenvalues as a polynomial, the eigenvalues for different values of the parameters can be obtained with little additional computational effort once the coefficient matrices are determined. 2.3. l0 Convergence of the perturbation expansion can be improved if the quantity >bLd > becomes smaller. Because >bLd > = maxa (=ga =), where the ga are deviations of the almost congruent eigenvalues, and since the ga depend on the choice of the primary unperturbed eigenvalue l0 , it is possible to find a l0 which will minimize >bLd >. A simple optimization procedure follows. First, suppose that the primary unperturbed eigenvalue is temporarily set to the first unperturbed eigenvalue of the set. All the elements of this set can be expressed as la = l 0 + m˜a i + g˜ a = l0 + ma i + ga ,
(17)
where m˜a = ma + m; l0 is the primary unperturbed eigenvalue; and ma and m are integers, with a = 1,..., J. The tilde superscript indicates current values, and the terms without superscripts indicate optimized values. Since >bLd > = maxa (=ga =), the optimum value of l0
251
and the corresponding values of m are found by minimize {maximum (=ga =)}, l 0 ,m
where =ga = = =l 0 + mi + g˜ a − l0 =.
(18)
a
3. STABILITY BOUNDARY
Since the real parts of the Floquet eigenvalues determine stability, the perturbation expansion can be used to establish the transition between stable and unstable behavior in the o–v and o–d planes. Note that the perturbation expansion for Q needs to be calculated only once to be useful over a range of o and v values. The procedure for the stability boundary follows. First, a rectangular, evenly spaced mesh on the parameter plane is generated. Second, the mid-point of the frequency range is used as the frequency of the unperturbed system. Third, the coefficients in the perturbation expansions of F and Q are found. Finally, the Floquet eigenvalues are determined at each mesh point. If the real part of any of the Floquet eigenvalues is greater than zero, then the grid point corresponds to an unstable response. A practical limitation of the method is that if the specified range of v is too wide, then the maximum value of d will also be large and all unperturbed eigenvalues will be included in a congruent set. As a result, the expansion of F and Q may not converge. To avoid this difficulty, the range of v is divided into a number of subdivisions so that the maximum value of d in each division satisfies (p) (p) (p) (19) =dmax = Q =omax =max (>B(p) s >, >As >)/max (>L>, >Bs >, >As >). p,s
p,s
Because of this restriction, the influence of dmax on deciding the value of K in equation (13) is effectively eliminated. 4. EXAMPLES
4.1. The Mathieu equation dx d2x + 2z ˜ + (1 + o cos vt˜ )x = 0 dt dt˜ 2
(20)
is a classical model of parametric excitation, where o characterizes the strength of the excitation, and z is a damping coefficient. Equation (20) is rewritten in the first order form dx = dt˜
0$
i 0
% $
0 −1 +z −i i
% > 0$
−i +o 4 −1
%
$
% 11
1 i e−ivt˜ + −i 1
i 1
1 eivt˜ −i
x.
(21)
The system experiences a primary parametric resonance when v = 2 and, consequently, stability will be investigated for a range of frequencies centered about this critical value. Thus, let d = (v − 2)/2 and t = vt˜ . Equation (21) becomes dx (1 − d + d 2 ) = 2 dt
6$
i 0
% $
0 −1 +z −i 1
% $
−i o i + −i 4 1
%
7
1 (eit + e−it ) x, −i
(22)
where the higher order terms in d have been omitted. In equation (22), if z = 0, the system has the two exactly congruent unperturbed eigenvalues 0·5i and −0·5i. Let Ad = 12
$
−1 i
%
−i −1
and
b = z,
.-. .
252
to employ the approach outlined in section 2.2.2. Then, F and Q for this example will each be 2 × 2 matrices, and each can be expanded: Q = 12
F=
6$
% $
−i 0 −1 +z 0 −i 0
6$ % $ % 7 0 1
+
z 2
0 0 + 0 0
1 it o e + 0 16
6$ % $ % 7 1 0
0 0 + 0 0
% $
0 i +d −1 0
0 it e , 1
% $ %7
0 o 0 + −i 4 1
6$ % $ i 2
0 −it 0 e + 0 0
1 0
% $
,
−2 0 + 0 −2
(23a)
% $
0 it 0 e + 0 0
% 7
2 e2it −i
(23b)
to first order as in equation (15). After restoring the time scale, the Floquet eigenvalues of matrix Q become l = (−i − z 2 z(0·25o)2 − d 2 )(1 + d). (24) In Figure 1, values of the real part of the eigenvalue are plotted as functions of o for a range of v values. The agreement among the result (24) and the eigenvalues calculated numerically from the transition matrix is typically good. For instance, when v = 2·3 (d = 0·15), the difference between the two solutions is within 5%, even though O(o) = 1. Because of the simple form of equation (24), one may directly determine the equation for the stability boundary by requesting Re(l) = 0. The line demarcating the stability zone is given by o 2 = 16(d 2 + z 2 ). (25) In Figure 2 are shown the stability boundaries (25) compared with the numerical results for four different values of the damping coefficient. The error in the perturbation method increases with o, d and z, but in general the agreement between the two methods is good. The damped Mathieu equation is treated by Nayfeh and Mook [24]. There, the damping constant is proportional to o, and the Lindstedt–Poincare´ technique is used to develop a
Figure 1. The real part of the Floquet eigenvalues in Example 1 for v = 2·0, 2·1, 2·2, 2·3 and 2·4: W, numerical integration, ——, perturbation.
883844—JSV 182/2
253
Figure 2. Stability boundaries in Example 1 for z = 0·0, 0·05, 0·1 and 0·2: W, numerical integration; ——, perturbation.
perturbation expansion for the solution. This approach is appropriate if the solution is periodic; namely, o and d lie on the stability boundary. The solution in reference [24] for the stability boundary is, except for differing notation, identical with equation (25). Lion [22] uses a perturbation approach which is similar to that presented here in that it treats a first order system. The perturbation solution for the undamped Mathieu equation, however, does not take into account the existence of congruent eigenvalues. In terms of the current notation, reference [22] predicts the Floquet eigenvalues l = 2i(2/v + o 2/2v(v 2 − 4)),
(26)
and the stability boundary is found by using the criterion that on the stability boundary the solution has the same periodicity as the coefficients. Together with v = 2(1 + d), this implies that the stability boundary is given by o 2 = 2(v − 2)2(v + 2) 1 32d 2.
(27)
In contrast, equation (25) gives the result that o 2 1 16d 2 when z = 0. From Figure 2, the present method yields the correct asymptotic behavior when z vanishes and d is small, and equation (27) is proportionally off by a multiplicative factor of z2. In addition, equation (26) predicts that only the imaginary part of the eigenvalue is affected by the perturbation o. The results depicted in Figure 1 indicate that this is, in fact, not the case and that the real part of the eigenvalue also changes. It should be noted that Lion’s method was not intended for the case of congruent eigenvalues. This is the reason that equation (26) becomes singular as v = 2, and the method does not accurately predict the behavior of the eigenvalues near primary resonance. However, the behavior near resonant points is critical, since this is the region where vibration amplitudes can be the greatest. 4.2. The trend in drive train design for marine applications is towards lightweight, high performance systems. A typical drive train consists of an engine, flywheel, gearbox, shaft and propeller. The power to weight ratio may be increased, for example, by increasing the size of the pistons and by reducing the mass of the flywheel. Since the inertia of the piston changes periodically with time, this design trend tends to increase the ratio of the variable
254
.-. .
inertia to the constant inertia terms in the equation of motion; as a result, instability from parametric excitation becomes a concern. In the general formulation presented here, it is assumed that the inertia, damping and stiffness terms may vary with the angular position of the shaft. In fact, both the stiffness and inertia terms vary periodically if a U-joint is used to offset part of the drive train. Denoting the rotational speed of the drive train by v, and linearizing the governing equations for small amplitude vibration U(t˜ ), the equation of motion becomes [25] (M + oM (t˜ ))U (t˜ ) + (C + voC (t˜ ))U (t˜ ) + (K + v 2oK (t˜ ))U(t˜ ) = E(t˜ ) + v 2oE (t˜ ).
(28)
The parameter o is small and indicates the relative magnitude of the periodically varying terms. In the case of pistons, o = r/(2 + r), where r is the ratio of the varying moment of inertia of the connecting rod and piston to its average value. It is assumed that the matrices M , C , K and the vectors E and E are periodic with a period corresponding to one cycle of the engine, namely 2p/v. The homogeneous form of equation (28) can be transformed into the form of equation (1): [I + oA 1 (t˜ )]x˙(t˜ ) = [L + oB 1 (t˜ ) + ovB 2 (t˜ ) + ov 2(t˜ )B 3 (t˜ )]x(t˜ ),
(29)
with x(t˜ ) = [U(t˜ ), U (t˜ )]T. A 1 , B 1 , B 2 and B 3 are periodic and the elements of the diagonal matrix L are the unperturbed eigenvalues of the system. From equations (3) and (13), a particular perturbation expansion to determine the Floquet eigenvalues holds only over a limited frequency range. A practical consideration which becomes apparent when dealing with large systems is the best choice of the frequency range in the expansion. For example, consider a boat drive with motion characterized by 11 rotational degrees of freedom, and parameters o E 0·5 and 2700 E v E 3000. One approach would be to try to find a single perturbation expansion that would cover the entire frequency range from 2700 to 3000 rad/s. If this approach is used, then from equation (3), v0 = 2850 and =d= E 0·053. Inspection of the matrices involved indicates that =d= determines the value in equation (13) of K = 0·158 and, as a result, all 22 eigenvalues must be grouped together as a congruent set. This adversely affects convergence, since all eigenvalues belong to one congruent set. Instead, if the frequency range is divided into three equal sub-intervals then the term involving =o= in equation (13) dominates, and K is reduced to 0·096; with this choice, the maximum number of eigenvalues that must be grouped together in any of the frequency intervals is eight. As a result, the term b >Ad > is reasonably small and the convergence rate of the perturbation expansion is enhanced. From a practical point of view, the procedure discussed in section 3 is always utilized to restrict the frequency range in accordance with equation (19) so that d is sufficiently small so that o determines K. In this example, the Floquent eigenvalues of the system were calculated for the boat drive train using both the perturbation method, and the numerical transition matrix approach for 2000 E v E 4000 and o E 0·5. The transition matrices were calculated using a fourth order Runge–Kutta algorithm to numerically integrate the equations of motion over one period of excitation. The results are shown in Figures 3 and 4. Both methods indicate similar instability regions. As might be expected, the perturbation method accurately predicts the stability boundary for smaller values of o, but it does not do as well when o is large. The perturbation method provides different expansions in each different frequency range and, as a result, it erroneously predicts that the stability regions change abruptly from one frequency range to the next. The magnitude of the frequency range used in this example was determined by the long computation times for the transition matrices if a lower, more realistic frequency
255
Figure 3. Unstable regions (shaded) in Example 2 for 2000 Q v Q 4000 (perturbation).
range was used. This difficulty occurs because the time step in the Runge–Kutta algorithm must be sufficiently small so that the contribution of the mode with the highest natural frequency can be accurately integrated over the excitation period. Consequently, numerical integration times would increase by almost an order of magnitude if v were decreased to the actual operating range of the drive train, say 50–300 rad/s. In contrast, computer time was not a factor in calculating the perturbation solutions, which took approximately a factor of 200 less time on the computer than the Runge–Kutta method. Furthermore, the time required in calculating the perturbation solution decreases when the system’s stability is analyzed for the actual, lower frequency range. A final observation concerns boundary tracing techniques for large-dimensional systems. When the number of eigenvalues is large, there is a significant number of para-
Figure 4. Unstable regions (shaded) in Example 2 for 2000 Q v Q 4000 (numerical integration).
256
.-. .
metric resonances that can occur in a given frequency range. As a result, there is an abrupt hange in the slope of the stability boundary when the different resonance regions overlap. Furthermore, a number of higher order resonances can result in very narrow stability regions, rendering the boundary tracing technique difficult to implement.
5. SUMMARY
A new perturbation method has been developed for calculating the Floquet eigenvalues and eigenvectors of a system of ordinary differential equations with coefficients that vary periodically in time. The method is useful in dealing with large systems because it directly treats the cases of almost congruent and exactly congruent eigenvalues. This can be an important consideration to the extent that the likelihood of eigenvalues being nearly congruent increases when there are more eigenvalues in the system. The method expresses the solution in terms of a multi-parameter expansion where the parameters typically correspond to a parameter which characterizes the magnitude of the periodic coefficients, a parameter which characterizes the frequency variation, and a third parameter which reflects the closeness of the nearly congruent eigenvalues. The third parameter is selected so as to make all of the nearly congruent eigenvalues exactly congruent and, then, a perturbation expansion is applied which is appropriate for subsets of exactly congruent eigenvalues. It is found that by grouping the nearly congruent eigenvalues together in this manner, the radius of convergence of the expansion is significantly improved. Subsets of exactly congruent eigenvalues should be grouped together. If there are M elements in the subset, then the eigenvalues can be expressed as an expansion in which the coefficients are M × M matrices. Once the coefficient matrices in the perturbation expansions are determined, the eigenvalues can easily be calculated for different values of the parameters. This result in an efficient computational method for larger systems. The method is applied to the damped Mathieu equation, and to torsional drive train vibration. In the case of the damped Mathieu equation, the first order solution gives a single analytical expression for the eigenvalues in terms of the amplitude of the excitation, the frequency and the amount of viscous damping in the system. The results agree reasonably well with those of independently generated numerical solutions. When the real part of the eigenvalue vanishes, the perturbation expansion predicts an equation for the stability boundary which is consistent with other perturbation methods and which also compares favorably with the numerical results. In the case of the boat drive train, the shape of the stable region is much more complicated because of complex interactions that occur when there is a large number of eigenvalues. The perturbation method predicts a stability region which is very similar to that determined numerically. As might be expected, the agreement is very good for small values of the perturbation parameters and worsens as the magnitudes of the perturbation parameters increase. The perturbation method is several orders of magnitude more computationally efficient than the numerical method.
ACKNOWLEDGMENT
The authors acknowledge partial support for this work from the United States Navy under contract N00189-90-C-0331. Mr James Brown is the contract monitor.
257
REFERENCES 1. V. A. Y and V. M. S 1975 Linear Differential Equations with Periodic Coefficients. New York: John Wiley. 2. S. K. S 1981 Journal of Applied Mechanics 48, 174–177. Stability theorems for multidimensional linear systems with variable parameters. 3. S. Z. A-R and G. A 1986 International Journal of Systems Science 17, 1654–1660. Stability analysis of non-autonomous linear systems by a matrix decomposition method. 4. R. K. Y and S. R. K 1988 International Journal of Systems Science 19, 1853–1858. Stability analysis of linear time-varying systems. 5. M. V. M 1990 Avtomatika i Telemekhanika 4, 27–35. An algorithm of analysis of the stability of linear periodic systems and its computer realization. 6. C. S. H 1974 Journal of Mathematical Analysis and Applications 45, 234–251. On approximating a general linear periodic system. 7. P. F, C. E. H and T.-H. W 1977 International Journal for Numerical Methods in Engineering 11, 1117–1136. Efficient numerical treatment of periodic systems with application to stability problems. 8. P. M. G, P. B, G. F and M. L 1987 NASA Conference; 58th Shock and Vibration Symposium 1, 487–494. Dynamic behavior of dissymmetric rotor bearings modelled with a periodic coefficient large system. 9. B. W and H. K 1991 Mechanism and Machine Theory 26(2), 123–131. Direct Floquet method for stability limits determination—I: theory. 10. H. K and B. W 1991 Mechanism and Machine Theory 26(2), 133–144. Direct Floquet method for stability limits determination—II: application and phenomena. 11. J. W. L and I. C 1990 American Institute of Aeronautics and Astronautics Journal 28, 1089–1097. Stability sensitivity analysis of a helicopter rotor. 12. A. U and M. A. B 1981 Journal of Applied Mechanics 48, 948–958. Parametric instability of a rotating shaft due to pulsating torque. 13. I. G. T and C. J. Y 1985 International Journal of Mechanical Science 27, 813–822. Effect of flexibility on the dynamic stability of a four-bar mechanism. 14. D. R. F, J. and H. T. J. L 1988 20th Biennial Mechanisms Conference: Trends and Developments in Mechanisms, Machines, and Robotics, 397–405. Stability analysis by a perturbation method. 15. S. L. L, Y. K. C and S. Y. W 1982 Journal of Applied Mechanics 49, 849–853. A variable parameter incrementation method for dynamic instability of linear and nonlinear elastic systems. 16. B. W 1988 Proceedings of the 2nd International Symposium on Transport Phenomena, Dynamics and Design of Rotating Machinery, Honolulu, 303–314. Rapid stability analysis of parameter-excited large systems. 17. C. S. H 1963 Journal of Applied Mechanics 30, 367–372. On the parametric excitation of a dynamic system having multiple degrees of freedom. 18. C. S. H 1965 Journal of Applied Mechanics 32, 373–377. Further results on parametric excitation of a dynamic system. 19. T. Y and A. S 1968 Bulletin of the Japan Society of Mechanical Engineers 11(43), 92–101. On the oscillations of ‘‘summed and differential types’’ under parametric excitation (vibratory systems with damping). 20. S. T. A and N. S N 1986 Journal of Sound and Vibration 107, 215–230. Dynamic stability of pipes conveying pulsating fluid. 21. S. T. A and S. F. A 1987 Journal of Mechanisms, Transmissions, and Automation in Design 109, 412–418. Dynamic stability of chain drives. 22. P. M. L 1966 Journal of The Franklin Institute 281(1), 27–40. Stability of linear periodic systems. 23. Y. N and P. Y. W 1969 Journal of The Franklin Institute 287(2), 143–157. A method for stability investigation of a periodic dynamic system with many degrees of freedom. 24. A. H. N and D. T. M 1979 Nonlinear Oscillations. New York: John Wiley. 25. W.-T. W 1993 Ph.D. Thesis, Department of Mechanical Engineering, Carnegie Mellon University. Stability analysis and forced response of reciprocating mechanisms.