- Email: [email protected]
Dynamic stability of periodic solutions of large scale nonlinear systems
COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
48 (1985) 191-202
DYNAMIC STABILITY OF PERIODIC SOLUTIONS OF LARGE SCALE NONLINEAR SYSTEMS
H.G. MATTHIES
and Chr. NATH
Germanischer Lloyd AG, D-2000 Hamburg 11, West Germany Received 6 October 1983 Revised manuscript received 21 May 1984
A method suitable for the stability analysis of periodic solutions of large nonlinear systems is presented. The stability behaviour is inferred from a Galerkin approximation of the associated Poincare map. The practical value of the method is shown by applying it to the nonlinear finite element model of a large scale wind energy convertor.
1. Introduction
For systems with periodic solutions, such as mechanical systems with rotating parts, it is often important to establish the conditions under which this periodic solution is stable. The term dynamic stability appears to be redundant, since in any system loss of stability is a dynamical phenomenon, but in many cases it can be investigated by considering the statics of the system only, e.g., buckling analysis of an elastic structure under conservative loading. In the study of periodic solutions, this is not possible and the dynamical behaviour of the system becomes important. The method for dynamic stability analysis is given in Section 2 and is based on the theory of periodic differential equations. A vast amount of literature, both from the mathematical and an engineering viewpoint, exists on this subject, particularly for linear equations (e.g. [l-8] and the references therein). To be specific, we shall phrase the results in terms of a mechanical system, but the theory and numerical approximation are equally applicable to any other dynamical system. The stability of periodic solutions, which we shall investigate, has to be distinguished from and is more subtle than the phenomenon of resonance, in which case most probably a periodic solution will not exist at all. Since the systems we want to consider are described by a large number of equations, a full stability analysis seems infeasible and an approximate procedure must be introduced. We propose a Galerkin approximation in Section 3 with the natural vibration modes as basis vectors, and the stability analysis is performed for the reduced system. The last section contains the application of the theory to the nonlinear finite element model of the wind energy convertor GROWIAN (GroBe Windenergieanlage). Due to a lack of experience in the construction and operation of such very large wind energy conversion systems the global dynamic stability analysis emerges as one of the main issues besides the stress and vibration analyses. 00457825/S/$3.30
@ 1985, Elsevier Science Publishers B.V. (North-Holland)
192
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
2. Stability of periodic solutions
Let us take as the starting point the equations of motion of a mechanical system with a finite number of degrees of freedom. When dealing with a continuous system, we shall assume that an appropriate discretisation has been introduced, e.g. with finite elements. Our definition of stability, which is Liapunov’s, is formulated with reference to the norm topology on a finite dimensional space. Since, in contrast to the finite dimensional case, different norms on an infinite dimensional space are in general not equivalent, there is no unique extension of Liapunov’s definition to continuous systems. The concept of stability is norm dependent and even definitions involving two different norms are being proposed 171. In order to circumvent these di~culties, we shall simply assume that the stability of the continuous system can be judged by considering the discretised system. Consider now the equations of motion of a mechanical system, Mti + S(w, I+) = F(w, t) .
(1)
Here M is the mass matrix, S the vector of internal forces, generally dependent on the displacement w and velocity k, and F the vector of external forces, possibly displacement dependent. Let the external force be periodic with period 77 and assume that (1) has a periodic solution w(t), i.e., w(T)= w(O)= If y(t) is an arbitrary equation
u*
and
solution
G(T) = lit(O)= t1’. to (l), then the difference
Mi’+S(z+w,i+~+)-S(w,G)=F(z+w,t)-F(w,t),
(2)
z(t) = y(t)-
w(t) satisfies the
(3)
which has the origin z(t) = o as a solution. The periodic solution w(t) of (1) is called stable or asymptotically stable respectively according to whether the origin is stable or asymptotically stable for (3) (cf. [3]). Stability and asymptotic stability are defined following Liapunov (cf. [l-3]), and we recall the definitions for the sake of completeness: The origin is ~fu~Ze for (3) if for any bound on the solution we can find corresponding bounds for the initial conditions, such that the solution will stay within the given bounds, i.e., i ( t )I]< r for all t > 0 if llz(0)Jl-t Ili(O)l]< s. for all r > 0 there is a s > 0 such that flz(t)ll+ 11 The origin is asymptoticdly stable for (3) if it is stable and additionally the solution tends to zero if the initial condition is small enough, i.e., there is a s > 0 such that /z(t)ll+ Ili(t)ll=+0 for t + 02 if llz(O)ll + Ili(O)ll< s. The customary way of dealing with the stability of (1) or (3) is to consider the variational equation (the linearisation or first approximation) Mi+Ci+Kz=o,
(4
where C(w, *it)= as/&& is the tangential damping matrix and K(w, *it)= ~S/~w - ~F/~w is the tangential stiffness matrix. Both are periodic with period T (since w and 3 are), and hence this
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
193
is a linear equation with periodic coefficients. According to the Floquet-Liapunov theory (cf. [3, S]) the knowledge of the monodromy or transition matrix over one period is sufficient to determine stability or instability of (4). The monodromy matrix B is defined by the relation
for arbitrary initial conditions (z(O), i(0)) of the solution z(t) to (4). Equation (5) defines a linear time-discrete dynamical system, and the origin is asymptotically stable, stable or unstable for (4) if and only if the origin is asymptotically stable, stable or unstable for (5). The definitions of stability for (5) are essentially a rewording of the above stability definitions for (3). Let p(B) = max{(Ail: Ai eigenvalue of B} denote the spectral radius of B; then it is well known that (cf. [l-3, 81: (i) The origin is asymptotically stable for (5) if and only if p(B) < 1. (ii) If p(B) > 1, the origin is unstable (not stable) for (5). (iii) If p(B) = 1, the origin is unstable for (5) unless for all eigenvalues on the unit circle the corresponding Jordan block is of size 1 in the Jordan decomposition of B, in this case the origin is stable but not asymptotically stable. Since stability is a local property and the nonlinear contributions of higher order in (4) vanish, the above conditions (i) and (ii) carry over to (3) and hence to the periodic solution w(t) of (1) (cf. [3]). Condition (iii) itself is not stable under perturbations; the stability behaviour is determined by the nonlinear terms. In the case of an autonomous system (i.e. F = o), certain refinements are necessary (cf. [l, 3]), since always p(B) = 1. Here we have the desired stability if A = 1 is a simple eigenvalue, all other eigenvalues lying strictly inside the unit circle. Summarizing, (1) was ‘differentiated’ to arrive at (4) which in turn led to the time-discrete system (5) used to determine stability. Since we want to organize the numerical approximation somewhat differently, we consider another approach which essentially interchanges the two steps ‘differentiation’ and ‘time-discretisation’. Let us define the Poincare map associated to (1) by
where u(t) is the solution of (1) with initial conditions (u(O), zi(0)). In other words, the Poincare-map transforms a given initial state into the final state after one period. Under reasonable assumptions (cf. [2, 31) the solution u(t) of (1) is unique and cannot blow up in finite time. Hence, the Poincare map P in (6) is well defined for all initial conditions in phase space. As is evident from (2) (1) h as a periodic solution w(t) with initial conditions (u’, 0’) if and only if (u’, v”) is a fixed point for the map P in (6) i.e., initial and final state after one period are identical. Again under reasonable assumptions (cf. [2, 31) a solution u(t) of (1) and its derivative ri(t) are differentiable w.r.t. the initial conditions (u’, v”). As is well known from standard results of the theory of ordinary differential equations (differentiation w.r.t. initial conditions), the quantities &/au”, W&O are matrix solutions of the variational equation (4). This may be seen by formal differentiation of (1) (justified under suitable hypotheses on S and F). It is easily verified that
194
H.G. Matthies,
a(u(T), rip)) ~(~(O), i(0)) =
Chr. Nath, Dynamic
stability
of periodic solutions of large systems
(7)
B,
due to the linearity of (3), but on the other hand from (6) we may infer Q(T), Ii(T)) a(a(0), zi(0)) =
iq$llti(O)) ’
This shows that B is the derivative (Jacobian) of P at the point of linearisation (u’, DO),and the linear time-discrete system (5) is the first approximation or linearisation of the nonlinear time-discrete system (6). The stability de~nitions for (6) are analogous to those fur (3) (only the continuous time has to be replaced by discrete time points), as well as the stability conditions (i) and (ii) in terms of the spectral radius of B (cf. [I]). These conditions can now be seen in a different light; they indicate whether the fixed point (st’, u”) is attractive or repellent and are well known from the theory of successive approximations~ Attractive fixed points are indicated by a spectral radius of B less than unity, and they correspond to asymptoticafiy stable solutions. If on the other hand there is an eigenvalue outside the unit circle, the fixed point is repellent and corresponds to an unstable solution of (1).
3. Numerical approximation In systems with many degrees of freedom, as they are encountered in the discretisation of partial di~erential equations for example, the mere calculation of the monodr~my matrix B would require a tremendous numerical effort, not to speak of the subsequent computation of the spectral radius of B, Additiona~Iy, one often wants further information on the stability behaviour, and several of the largest eigenvalues and corresponding eigenvectors of B are required, thus further increasing the computational effort. One has to recall that B has a dimension of twice the number of kinematic degrees of freedoms and is generally fully populated and not symmetric. A way to come around this obstacle is to introduce a Galerkin approximation: (9) As in any Calerkin approximation, the accuracy will depend on how wet1 the function u(t) can be approximated from the subspace spanned by the basis-vectors (pi]. For mechanical systems, a natural choice is the basis of the M-orthonormal eigenvectors of the generalized eigenproblem
is the symmetric part of 1%
195
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
In many cases most of the response is contained in the first few eigenmodes, hence the subspace spanned by these {vi} may be chosen of much smaller dimension than the original number of kinematic degrees of freedom, thereby greatly reducing the size of the problem. Denoting by
(11)
@k= [Vl, . . .? Vkl a matrix with the first k eigenvectors as columns, and projecting on the subspace spanned by the {Cpi},
the variational
@:MfDkti+ @:CQi,ti + @:K@ka = o .
equation
(4)
(12)
with (Y= [al,. . . , cq.]‘, one may compute a reduced monodromy matrix for the linear periodic system (12). There are many variations on this approach; often it is used in conjunction with component mode synthesis where part of the structure is fixed and another part performs a periodic motion and the periodic matrices M, C and K are hereby usually approximated by their first harmonic. In any case the response is probably well approximated by span{qi} at t = 0, but this may not be the case for other times t # 0. In order to avoid this inaccuracy and the numerical effort involved in the computation of the projected matrices and the fundamental solution to (12) we shall instead adopt the view that the monodromy matrix B is the derivative of the Poincare map P at the point (u’, u’), and we shall project this derivative onto the subspace spanned by the {vi}. In this way nothing more but the assumption that the response is well approximated by the subspace at t = 0 is used. To this end, we compute solutions ui(t) of (1) with initial conditions Ui(0)=
U” +
E#i
t&(O)=
7
i = 1,. . . , k,
v’,
(13)
and &+k(O) = TV’+ &i+kqi
&+/c(O) = 24’7
,
i = I, . . . , VI ,
(14)
to obtain the values of P at the points (Uj(O),tij(O)) (j = 1, . . . , m + k). The derivative of P restricted to the subspace span{(qp,, o), . . . , (&, o), (0, ql), . . . , (0, q,,,)} of phase space is then computed by numerical differentiation. Setting
z
ET’[u~(T)- u', t&(T)-
0’1’ ,
i = 1, . . . , k
(1%
and #i+k :=
B[O,
pi]’
z
&F+!;lk[Ui+k(T)-
u”,
hi+k(T)-
vO]‘,
i = 1, . . . , m ,
(16)
the reduced monodromy matrix B k,,, is computed as the Galerkin approximation to B by projecting onto the above subspace. The vectors ((pi, o) and (0, vi) form an orthonormal basis if we take the norm in phase space as
lb, v)ll’ =
U'MU
+ u'Mv .
196
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
Hence, Bkmis given by (17) Since Bkm is typically of much smaller dimension than B, the computation of all of the eigenvalues is no major numerical effort. If additionally the eigenvalues are computed for an increasing sequence of k and m, convergence of the eigenvalues, which is usually very quick, may be monitored. Usually very few qpi are necessary to obtain a good approximation. The projection of B onto the subspace spanned by the {vi} can be regarded as a change of basis, so if the original system is finite dimensional we will certainly obtain a good approximation by choosing k and m large enough. If on the other hand (1) is a discretisation of a continuous system, the difference between Bkmand B becomes immaterial, since the difference between Bkmand the monodromy operator of the continuous system is now important. Any error estimates available for the monodromy matrix B in (5) may in this case be used to yield error estimates for Bkm.
4. Time discretisation
and application
of the method to a large wind energy convertor
The methods to determine stability or instability as presented in the previous section assume that it is possible to compute the final state (w(T), G(T)) for a given initial state (u’, u’). Analytical computation of the final state will not be possible for most applications and one has to consider a numerical approximation in time of the differential equation (1). It is clear that for a dynamic stability analysis the response has to be computed very accurately, otherwise the numerical differentiation in (1.5) may not yield any useful results. Thus, the choice of time-stepping algorithm and time step size becomes very important. There is a vast body of literature on the merits of various time stepping algorithms (e.g., cf. [9-121 and the references therein). We have to recall that the considerations about the stability properties and approximation order are statements about the asymptotic behaviour, i.e., when the time step At + ~0 or At + 0. But in an actual computation one has to employ a finite step size At, large enough to make the solution feasible at all and small enough to ascertain a desired accuracy. Certainly any numerical scheme used should be convergent as At + 0, but the asymptotic approximation order alone is not a good criterion for the selection of an algorithm. In theory, the advantage of a higher-order method should be to allow a larger step size than a lower-order method for a comparable accuracy, but this is often not the case. In fact, the stability properties may become the decisive aspect. For linear problems, the trapezoidal rule (or Newmark average acceleration method) is most popular, as it has no numerical damping, is unconditionally stable and is of second order. It has been known for some time, however, that the conclusions drawn from linear problems do not necessarily carry over to the nonlinear regime, (cf. [lo, 121). Th e main difference is that for an implicit method the linear equations may be solved almost to machine accuracy, whereas in nonlinear problems one has to employ some iterative scheme and, hence, there will be an equilibrium imbalance due to the termination criterion for the iteration. These error forces may cause
H.G. Matthies, Chr. Nath, Dynamic stabilityof periodic solutions of large systems
197
severe numerical problems for the trapezoidal rule {cf. [IO, 121). Backward differentiation formulas are not plagued by these spurious oscillations (cf. 111, 121) but the one-step backward formula has only first-order accuracy. We have conducted a series of tests to establish the properties of different methods for the problem at hand, and they show that for reasonable time step sizes the approximation order is irrelevant, whereas the stability is decisive. The test system was a very coarse model of a wind energy convertor (5 nodes, 4 beam elements, 12 degrees of freedom). Initial conditions were chosen corresponding to a steady rotation of the rotor. Since damping was not present, the system should stay in a state of steady rotation indefinitely. We compared the computed values for the rotational speed and axial extension of the rotor blades with the exact values, which are constant in the case, in Figs. 1 and 2. The methods which were tested are the backward Euler method (first-order backward di~erentiation), the trapezoidal rule, Wilson’s 0 method (cf. [14]) with B = 1.4, Newmark’s method with y = 0.5, /3 = 0.4 and with y = 1.0, j3 = $(y + 4)“. The backward Euler method and the trapezoidal rule are from the following family of methods: z&+~= zi” + At(aii,+~ + (1 un+l = u, +
At(c&+l+
(1
a)&) ,
- ~l)ti,,)
(18)
= U, + A& + At2(02ii,+l + (1 - CY)(Y&) . Here u, is the value of u at the nth step. For rr = 1 we obtain the Euler backward method and for a! = OS the trapezoidal rule. The well-known Newmark family of methods is defined by I-backwordEuler I-trapezoidal rule 3-WllSOrl 0 (0=1.4 L-Newmark (Y=lO.P=O.561 I-Newmork lV=O.5.0 =0 LOI
-3 -4.-5.. -6
1
Fig. 1. Relative error in rotational speed (340 steps per revolution).
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
198
--3-Wilson 9, > .-
Ei
z :
1..
-
____----
___---
-\-
-'
/
.' '
x.
-VP.
-
-
-
-
-
.
-1.~
.5-Newmorklv=0.5,P=~LI I'
'\
,--*
\
1'
.
' '1
.
-2 I
\ .
\
.,
1
2
3
4
5
/
\ . .\
-51
- - --2-trapezoidal rule
\
'.
@=1.4)
\
/
\
0.
_‘_ /
__--
e
6
7
)-bockwardEuler
8
9 time
10
step
Fig. 2. Relative error in axial extension (360 steps per revolution).
(in+1= zi,+ At(yii,+l + (1 - y)iin) , un+l = u,, + At&, + At*(Pz&+l+ ($+ @ii,,).
(19)
For y = 0.5, /3 = 0.25 we again have the trapezoidal rule. Figs. 1 and 2 show that at this step size only the backward Euler method gives acceptable results, whereas those from the other methods are contaminated by parasitic oscillations and consequent unacceptable errors. The error in the backward Euler method has the same order of magnitude as the truncation error in the initial data. The step size is such that 360 steps are needed for one full revolution. By numerical experiment it was found that with the trapezoidal rule one would need approximately one thousand times as many steps to achieve the accuracy of the backward Euler method. For these reasons the backward Euler method was chosen for the stability computations with 360 steps per revolution. The finite element model of the wind energy converter GROWIAN and the main dimensions are shown in Fig. 3. GROWIAN has been erected at the mouth of the river Elbe in northern Germany and is presently undergoing final tests before operation. With a tower height of 100 m and a rotor diameter of also 100 m it is the largest of its kind. The finite element model consists of nonlinear beam and truss elements for the tower, the rotor and the parallel wire strands and the linear beam elements to model the stiffness of the foundation. To fulfill vibrational requirements, i.e., lower eigenfrequencies not to be resonant with multiples of rotor shaft revolutions at service speeds, each of the wire strands are given an initial tension force of 5300 kN. To reduce the vibrational input from the rotor the convertor is constructed with a teetering hub. This additional degree of freedom increases the risk of dynamic instability.
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
199
wind
generated height diameter mass
f e-model
model
power
of the large wind energy
3000
kW
100 m
of rotor
of nacelle
rotating
Fig. 3. Finite element
electric
of tower
100 m and
rotor
mass 160 degrees
converter
370 Mg 150 Mg
of freedom-13
elements
GROWIAN.
Equal attention as to the time-stepping algorithm has to be paid to the initial conditions. Transients from inaccurate initial data could propagate through a whole period (revolution) and hence invalidate the stability analysis. To avoid this, the initial conditions are taken from a quasi-static analysis using the same loads as for the dynamic case but now including explicitly the forces due to rotating masses. The deformations calculated in this manner correspond to the initial displacements of a stable rotating system at that instant (dynamic load step zero). The corresponding initial velocities are also easily computed under the assumption of a steadily rotating system. In general it will not always be easy to obtain initial conditions. A systematic approach would be to exploit the fixed point property of the Poincare map and employ some suitable fixed-point algorithm. The structure with applied initial conditions is given in Fig. 4. These initial conditions are now perturbed by the mode shapes of the lowest eigenfrequencies as described in Section 3. The 15 lowest eigenfrequencies in the state of rotation are shown in Fig. 5. The six vibration modes due to the wire strands were not deemed to be critical for the global stability, since they may be damped very easily, and all but one were therefore excluded from the stability investigation. The one retained (mode #lo: 0.95 Hz) was thought to be most critical in conjunction with the first bending mode of the tower (mode #2: 0.47 Hz). In Table 1 the moduli of the eigenvalues of the monodromy matrix are given for an increasing number of eigenmodes used in the perturbation. Obviously mode f4 excites a strong dynamic instability. Also the inclusion of the wire strand mode (the fifth mode in the stability analysis) raises the modulus of the eigenvalue corresponding with the tower bending slightly, but the eigenvalue stays inside the unit circle. The instability excited by mode R4 was also detected by parallel computations [15], but along the lines of the more conventional approach as described in Section 2. This instability means physically that the rotational speed of the rotor will not be stable. The rotor is driven by the wind forces and the power is extracted by the generator. In the finite element model used the generator is modeled very crudely as a simple damper of the rotational degree of freedom. In reality there is an elaborate feedback control system to keep the rotor at a constant speed. Computations of the control system show, however, that the rotation will be stable. Therefore, the mode f4 was excluded from the stability analysis, as well as the wire strand mode #lo, which did not
200
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
displacement scale 20 :l
Fig. 4. Structure
with applied
initial conditions.
contribute much. The moduli of the eigenvalues of the monodromy matrix for this perturbation pattern are given in Table 2 for an increasing sequence of modes. The eigenvalues are all strictly well inside the unit circle and the system is stable. In both cases the quick convergence of eigenvalues is remarkable.
0.88Hz P
1.71Hz
Fig. 5. Mode shapes
of the 15 lowest eigenfrequencies.
\
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
201
Table 1 Moduli of eigenvalues of I&,, (rotational and wire strand mode included) k (m =O)
1
2
3
4
5.
6
7
8
9
10
0.327
0.457 0.803
0.469 0.826 0.762
1.961 0.882 0.779 0.206
1.961 0.907 0.777 0.207 0.253
1.961 0.907 0.779 0.207 0.287 0.287
1.961 0.907 0.779 0.207 0.276 0.276 0.070
1.961 0.908 0.779 0.207 0.276 0.238 0.238 0.026
1.962 0.906 0.782 0.207 0.199 0.199 0.195 0.195 0.011
1.961 0.906 0.784 0.207 0.205 0.205 0.202 0.202 0.201 0.201
Table 2 Moduli of eigenvalues of I&,,, (rotational and wire strand mode excluded) k (m =0)
1
2
3
4
5
6
7
8
0.327
0.457 0.803
0.469 0.826 0.762
0.469 0.825 0.761 0.322
0.469 0.826 0.761 0.296 0.074
0.469 0.826 0.761 0.241 0.241 0.023
0.469 0.826 0.762 0.241 0.241 0.027 0.027
0.469 0.826 0.762 0.254 0.254 0.197 0.197 0.012
5. Conclusion The method presented for the dynamic stability analysis differs from the more conventional approach in that the two steps ‘linearisation’ and ‘time-discretisation’ are interchanged. In this form it is better suited for the stability investigation of large systems, along with some other advantages as described in Section 3. This is particularly true if the linearisation is only performed in a subspace spanned by a well-chosen set of basis vectors. Here, the vibration models of the lowest eigenfrequencies were used. This can be regarded as a Galerkin approximation of the original system. The choice of basis vectors for the subspace may be influenced, in that modes which are either not important or not well represented in the model may be excluded from the stability analysis. The method was applied to study the stability of the large wind energy convertor GROWIAN. The eigenvalues of the monodromy matrix numerically show a very quick convergence for an increasing sequence of subspaces.
202
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
References [l] J.P. LaSalle, The Stability of Dynamical Systems (SIAM, Philadelphia, PA, 1976). [2] M.W. Hirsch and S. Smale, Differential Equations, Dynamical Systems and Linear Algebra (Academic Press, New York, 1974). [3] N. Rouche and J. Mawhin, Ordinary Differential Equations-Stability and Periodic Solutions (Pitman, London, 1980) (translated from the French 1973 edition). [4] H. Ziegler, Principles of Structural Stability (Birkhauser, Base], 2nd ed., 1977). [S] H. Leipholz, Stability of Elastic Systems (Sijthoff & Noordhoff, Alphen aan den Rijn, 1980). [6] P. Friedmann, C.E. Hammond and Tze-Hsin Woo, Efficient numerical treatment of periodic systems with application to stability problems, Internat. J. Numer. Meths. Engrg. 11 (1977) 1117-1136. [7] H. Brauchli, On the norm dependence of the concept of stability, in: P. Germain and B. Nayroles, eds., Applications of Methods of Functional Analysis to Problems in Mechanics, Lecture Notes in Mathematics 503 (Springer, Berlin, 1976). [8] L. Collatz, Differettialgleichungen (Teubner, Stuttgart, 1970). [9] G. Dahlquist and A. Bjiirck, Numerical Methods (Prentice Hall, Englewood Cliffs, NJ, 1974). [lo] T.J.R. Hughes, Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics, Comput. & Structures 6 (1976) 313-324. [ll] Z. Zlatev and P.G. Thomsen, Application of backward differentiation methods to the finite element solution of time-dependent problems, Internat. J. Numer. Meths. Engrg. 14 (1979) 1051-1061. [12] A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations, Math. Comput. 28 (125) (1974) 145-162. [13] O.C. Zienkiewicz, The Finite Element Method (McGraw Hill, London, 3rd ed., 1977). [14] K.J. Bathe, Finite Element Procedures in Engineering Analysis (Prentice-Hall, Englewood Cliffs, NJ, 1982). [15] K. Kehl, W. Keim, F. KieBling and A. Rippl, Zur Dynamik von grol3en Windkraftanlagen, VDI-Berichte 381 (1980) 113-120. [16] J.H. Argyris, K.A. Braun, B. KirchgaBner and R. Walther, Statische und dynamische Untersuchungen an einem Windrotormodell, ISD-Bericht No. 259, Universitat Stuttgart, 1979.
IN APPLIED
MECHANICS
AND ENGINEERING
48 (1985) 191-202
DYNAMIC STABILITY OF PERIODIC SOLUTIONS OF LARGE SCALE NONLINEAR SYSTEMS
H.G. MATTHIES
and Chr. NATH
Germanischer Lloyd AG, D-2000 Hamburg 11, West Germany Received 6 October 1983 Revised manuscript received 21 May 1984
A method suitable for the stability analysis of periodic solutions of large nonlinear systems is presented. The stability behaviour is inferred from a Galerkin approximation of the associated Poincare map. The practical value of the method is shown by applying it to the nonlinear finite element model of a large scale wind energy convertor.
1. Introduction
For systems with periodic solutions, such as mechanical systems with rotating parts, it is often important to establish the conditions under which this periodic solution is stable. The term dynamic stability appears to be redundant, since in any system loss of stability is a dynamical phenomenon, but in many cases it can be investigated by considering the statics of the system only, e.g., buckling analysis of an elastic structure under conservative loading. In the study of periodic solutions, this is not possible and the dynamical behaviour of the system becomes important. The method for dynamic stability analysis is given in Section 2 and is based on the theory of periodic differential equations. A vast amount of literature, both from the mathematical and an engineering viewpoint, exists on this subject, particularly for linear equations (e.g. [l-8] and the references therein). To be specific, we shall phrase the results in terms of a mechanical system, but the theory and numerical approximation are equally applicable to any other dynamical system. The stability of periodic solutions, which we shall investigate, has to be distinguished from and is more subtle than the phenomenon of resonance, in which case most probably a periodic solution will not exist at all. Since the systems we want to consider are described by a large number of equations, a full stability analysis seems infeasible and an approximate procedure must be introduced. We propose a Galerkin approximation in Section 3 with the natural vibration modes as basis vectors, and the stability analysis is performed for the reduced system. The last section contains the application of the theory to the nonlinear finite element model of the wind energy convertor GROWIAN (GroBe Windenergieanlage). Due to a lack of experience in the construction and operation of such very large wind energy conversion systems the global dynamic stability analysis emerges as one of the main issues besides the stress and vibration analyses. 00457825/S/$3.30
@ 1985, Elsevier Science Publishers B.V. (North-Holland)
192
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
2. Stability of periodic solutions
Let us take as the starting point the equations of motion of a mechanical system with a finite number of degrees of freedom. When dealing with a continuous system, we shall assume that an appropriate discretisation has been introduced, e.g. with finite elements. Our definition of stability, which is Liapunov’s, is formulated with reference to the norm topology on a finite dimensional space. Since, in contrast to the finite dimensional case, different norms on an infinite dimensional space are in general not equivalent, there is no unique extension of Liapunov’s definition to continuous systems. The concept of stability is norm dependent and even definitions involving two different norms are being proposed 171. In order to circumvent these di~culties, we shall simply assume that the stability of the continuous system can be judged by considering the discretised system. Consider now the equations of motion of a mechanical system, Mti + S(w, I+) = F(w, t) .
(1)
Here M is the mass matrix, S the vector of internal forces, generally dependent on the displacement w and velocity k, and F the vector of external forces, possibly displacement dependent. Let the external force be periodic with period 77 and assume that (1) has a periodic solution w(t), i.e., w(T)= w(O)= If y(t) is an arbitrary equation
u*
and
solution
G(T) = lit(O)= t1’. to (l), then the difference
Mi’+S(z+w,i+~+)-S(w,G)=F(z+w,t)-F(w,t),
(2)
z(t) = y(t)-
w(t) satisfies the
(3)
which has the origin z(t) = o as a solution. The periodic solution w(t) of (1) is called stable or asymptotically stable respectively according to whether the origin is stable or asymptotically stable for (3) (cf. [3]). Stability and asymptotic stability are defined following Liapunov (cf. [l-3]), and we recall the definitions for the sake of completeness: The origin is ~fu~Ze for (3) if for any bound on the solution we can find corresponding bounds for the initial conditions, such that the solution will stay within the given bounds, i.e., i ( t )I]< r for all t > 0 if llz(0)Jl-t Ili(O)l]< s. for all r > 0 there is a s > 0 such that flz(t)ll+ 11 The origin is asymptoticdly stable for (3) if it is stable and additionally the solution tends to zero if the initial condition is small enough, i.e., there is a s > 0 such that /z(t)ll+ Ili(t)ll=+0 for t + 02 if llz(O)ll + Ili(O)ll< s. The customary way of dealing with the stability of (1) or (3) is to consider the variational equation (the linearisation or first approximation) Mi+Ci+Kz=o,
(4
where C(w, *it)= as/&& is the tangential damping matrix and K(w, *it)= ~S/~w - ~F/~w is the tangential stiffness matrix. Both are periodic with period T (since w and 3 are), and hence this
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
193
is a linear equation with periodic coefficients. According to the Floquet-Liapunov theory (cf. [3, S]) the knowledge of the monodromy or transition matrix over one period is sufficient to determine stability or instability of (4). The monodromy matrix B is defined by the relation
for arbitrary initial conditions (z(O), i(0)) of the solution z(t) to (4). Equation (5) defines a linear time-discrete dynamical system, and the origin is asymptotically stable, stable or unstable for (4) if and only if the origin is asymptotically stable, stable or unstable for (5). The definitions of stability for (5) are essentially a rewording of the above stability definitions for (3). Let p(B) = max{(Ail: Ai eigenvalue of B} denote the spectral radius of B; then it is well known that (cf. [l-3, 81: (i) The origin is asymptotically stable for (5) if and only if p(B) < 1. (ii) If p(B) > 1, the origin is unstable (not stable) for (5). (iii) If p(B) = 1, the origin is unstable for (5) unless for all eigenvalues on the unit circle the corresponding Jordan block is of size 1 in the Jordan decomposition of B, in this case the origin is stable but not asymptotically stable. Since stability is a local property and the nonlinear contributions of higher order in (4) vanish, the above conditions (i) and (ii) carry over to (3) and hence to the periodic solution w(t) of (1) (cf. [3]). Condition (iii) itself is not stable under perturbations; the stability behaviour is determined by the nonlinear terms. In the case of an autonomous system (i.e. F = o), certain refinements are necessary (cf. [l, 3]), since always p(B) = 1. Here we have the desired stability if A = 1 is a simple eigenvalue, all other eigenvalues lying strictly inside the unit circle. Summarizing, (1) was ‘differentiated’ to arrive at (4) which in turn led to the time-discrete system (5) used to determine stability. Since we want to organize the numerical approximation somewhat differently, we consider another approach which essentially interchanges the two steps ‘differentiation’ and ‘time-discretisation’. Let us define the Poincare map associated to (1) by
where u(t) is the solution of (1) with initial conditions (u(O), zi(0)). In other words, the Poincare-map transforms a given initial state into the final state after one period. Under reasonable assumptions (cf. [2, 31) the solution u(t) of (1) is unique and cannot blow up in finite time. Hence, the Poincare map P in (6) is well defined for all initial conditions in phase space. As is evident from (2) (1) h as a periodic solution w(t) with initial conditions (u’, 0’) if and only if (u’, v”) is a fixed point for the map P in (6) i.e., initial and final state after one period are identical. Again under reasonable assumptions (cf. [2, 31) a solution u(t) of (1) and its derivative ri(t) are differentiable w.r.t. the initial conditions (u’, v”). As is well known from standard results of the theory of ordinary differential equations (differentiation w.r.t. initial conditions), the quantities &/au”, W&O are matrix solutions of the variational equation (4). This may be seen by formal differentiation of (1) (justified under suitable hypotheses on S and F). It is easily verified that
194
H.G. Matthies,
a(u(T), rip)) ~(~(O), i(0)) =
Chr. Nath, Dynamic
stability
of periodic solutions of large systems
(7)
B,
due to the linearity of (3), but on the other hand from (6) we may infer Q(T), Ii(T)) a(a(0), zi(0)) =
iq$llti(O)) ’
This shows that B is the derivative (Jacobian) of P at the point of linearisation (u’, DO),and the linear time-discrete system (5) is the first approximation or linearisation of the nonlinear time-discrete system (6). The stability de~nitions for (6) are analogous to those fur (3) (only the continuous time has to be replaced by discrete time points), as well as the stability conditions (i) and (ii) in terms of the spectral radius of B (cf. [I]). These conditions can now be seen in a different light; they indicate whether the fixed point (st’, u”) is attractive or repellent and are well known from the theory of successive approximations~ Attractive fixed points are indicated by a spectral radius of B less than unity, and they correspond to asymptoticafiy stable solutions. If on the other hand there is an eigenvalue outside the unit circle, the fixed point is repellent and corresponds to an unstable solution of (1).
3. Numerical approximation In systems with many degrees of freedom, as they are encountered in the discretisation of partial di~erential equations for example, the mere calculation of the monodr~my matrix B would require a tremendous numerical effort, not to speak of the subsequent computation of the spectral radius of B, Additiona~Iy, one often wants further information on the stability behaviour, and several of the largest eigenvalues and corresponding eigenvectors of B are required, thus further increasing the computational effort. One has to recall that B has a dimension of twice the number of kinematic degrees of freedoms and is generally fully populated and not symmetric. A way to come around this obstacle is to introduce a Galerkin approximation: (9) As in any Calerkin approximation, the accuracy will depend on how wet1 the function u(t) can be approximated from the subspace spanned by the basis-vectors (pi]. For mechanical systems, a natural choice is the basis of the M-orthonormal eigenvectors of the generalized eigenproblem
is the symmetric part of 1%
195
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
In many cases most of the response is contained in the first few eigenmodes, hence the subspace spanned by these {vi} may be chosen of much smaller dimension than the original number of kinematic degrees of freedom, thereby greatly reducing the size of the problem. Denoting by
(11)
@k= [Vl, . . .? Vkl a matrix with the first k eigenvectors as columns, and projecting on the subspace spanned by the {Cpi},
the variational
@:MfDkti+ @:CQi,ti + @:K@ka = o .
equation
(4)
(12)
with (Y= [al,. . . , cq.]‘, one may compute a reduced monodromy matrix for the linear periodic system (12). There are many variations on this approach; often it is used in conjunction with component mode synthesis where part of the structure is fixed and another part performs a periodic motion and the periodic matrices M, C and K are hereby usually approximated by their first harmonic. In any case the response is probably well approximated by span{qi} at t = 0, but this may not be the case for other times t # 0. In order to avoid this inaccuracy and the numerical effort involved in the computation of the projected matrices and the fundamental solution to (12) we shall instead adopt the view that the monodromy matrix B is the derivative of the Poincare map P at the point (u’, u’), and we shall project this derivative onto the subspace spanned by the {vi}. In this way nothing more but the assumption that the response is well approximated by the subspace at t = 0 is used. To this end, we compute solutions ui(t) of (1) with initial conditions Ui(0)=
U” +
E#i
t&(O)=
7
i = 1,. . . , k,
v’,
(13)
and &+k(O) = TV’+ &i+kqi
&+/c(O) = 24’7
,
i = I, . . . , VI ,
(14)
to obtain the values of P at the points (Uj(O),tij(O)) (j = 1, . . . , m + k). The derivative of P restricted to the subspace span{(qp,, o), . . . , (&, o), (0, ql), . . . , (0, q,,,)} of phase space is then computed by numerical differentiation. Setting
z
ET’[u~(T)- u', t&(T)-
0’1’ ,
i = 1, . . . , k
(1%
and #i+k :=
B[O,
pi]’
z
&F+!;lk[Ui+k(T)-
u”,
hi+k(T)-
vO]‘,
i = 1, . . . , m ,
(16)
the reduced monodromy matrix B k,,, is computed as the Galerkin approximation to B by projecting onto the above subspace. The vectors ((pi, o) and (0, vi) form an orthonormal basis if we take the norm in phase space as
lb, v)ll’ =
U'MU
+ u'Mv .
196
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
Hence, Bkmis given by (17) Since Bkm is typically of much smaller dimension than B, the computation of all of the eigenvalues is no major numerical effort. If additionally the eigenvalues are computed for an increasing sequence of k and m, convergence of the eigenvalues, which is usually very quick, may be monitored. Usually very few qpi are necessary to obtain a good approximation. The projection of B onto the subspace spanned by the {vi} can be regarded as a change of basis, so if the original system is finite dimensional we will certainly obtain a good approximation by choosing k and m large enough. If on the other hand (1) is a discretisation of a continuous system, the difference between Bkmand B becomes immaterial, since the difference between Bkmand the monodromy operator of the continuous system is now important. Any error estimates available for the monodromy matrix B in (5) may in this case be used to yield error estimates for Bkm.
4. Time discretisation
and application
of the method to a large wind energy convertor
The methods to determine stability or instability as presented in the previous section assume that it is possible to compute the final state (w(T), G(T)) for a given initial state (u’, u’). Analytical computation of the final state will not be possible for most applications and one has to consider a numerical approximation in time of the differential equation (1). It is clear that for a dynamic stability analysis the response has to be computed very accurately, otherwise the numerical differentiation in (1.5) may not yield any useful results. Thus, the choice of time-stepping algorithm and time step size becomes very important. There is a vast body of literature on the merits of various time stepping algorithms (e.g., cf. [9-121 and the references therein). We have to recall that the considerations about the stability properties and approximation order are statements about the asymptotic behaviour, i.e., when the time step At + ~0 or At + 0. But in an actual computation one has to employ a finite step size At, large enough to make the solution feasible at all and small enough to ascertain a desired accuracy. Certainly any numerical scheme used should be convergent as At + 0, but the asymptotic approximation order alone is not a good criterion for the selection of an algorithm. In theory, the advantage of a higher-order method should be to allow a larger step size than a lower-order method for a comparable accuracy, but this is often not the case. In fact, the stability properties may become the decisive aspect. For linear problems, the trapezoidal rule (or Newmark average acceleration method) is most popular, as it has no numerical damping, is unconditionally stable and is of second order. It has been known for some time, however, that the conclusions drawn from linear problems do not necessarily carry over to the nonlinear regime, (cf. [lo, 121). Th e main difference is that for an implicit method the linear equations may be solved almost to machine accuracy, whereas in nonlinear problems one has to employ some iterative scheme and, hence, there will be an equilibrium imbalance due to the termination criterion for the iteration. These error forces may cause
H.G. Matthies, Chr. Nath, Dynamic stabilityof periodic solutions of large systems
197
severe numerical problems for the trapezoidal rule {cf. [IO, 121). Backward differentiation formulas are not plagued by these spurious oscillations (cf. 111, 121) but the one-step backward formula has only first-order accuracy. We have conducted a series of tests to establish the properties of different methods for the problem at hand, and they show that for reasonable time step sizes the approximation order is irrelevant, whereas the stability is decisive. The test system was a very coarse model of a wind energy convertor (5 nodes, 4 beam elements, 12 degrees of freedom). Initial conditions were chosen corresponding to a steady rotation of the rotor. Since damping was not present, the system should stay in a state of steady rotation indefinitely. We compared the computed values for the rotational speed and axial extension of the rotor blades with the exact values, which are constant in the case, in Figs. 1 and 2. The methods which were tested are the backward Euler method (first-order backward di~erentiation), the trapezoidal rule, Wilson’s 0 method (cf. [14]) with B = 1.4, Newmark’s method with y = 0.5, /3 = 0.4 and with y = 1.0, j3 = $(y + 4)“. The backward Euler method and the trapezoidal rule are from the following family of methods: z&+~= zi” + At(aii,+~ + (1 un+l = u, +
At(c&+l+
(1
a)&) ,
- ~l)ti,,)
(18)
= U, + A& + At2(02ii,+l + (1 - CY)(Y&) . Here u, is the value of u at the nth step. For rr = 1 we obtain the Euler backward method and for a! = OS the trapezoidal rule. The well-known Newmark family of methods is defined by I-backwordEuler I-trapezoidal rule 3-WllSOrl 0 (0=1.4 L-Newmark (Y=lO.P=O.561 I-Newmork lV=O.5.0 =0 LOI
-3 -4.-5.. -6
1
Fig. 1. Relative error in rotational speed (340 steps per revolution).
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
198
--3-Wilson 9, > .-
Ei
z :
1..
-
____----
___---
-\-
-'
/
.' '
x.
-VP.
-
-
-
-
-
.
-1.~
.5-Newmorklv=0.5,P=~LI I'
'\
,--*
\
1'
.
' '1
.
-2 I
\ .
\
.,
1
2
3
4
5
/
\ . .\
-51
- - --2-trapezoidal rule
\
'.
@=1.4)
\
/
\
0.
_‘_ /
__--
e
6
7
)-bockwardEuler
8
9 time
10
step
Fig. 2. Relative error in axial extension (360 steps per revolution).
(in+1= zi,+ At(yii,+l + (1 - y)iin) , un+l = u,, + At&, + At*(Pz&+l+ ($+ @ii,,).
(19)
For y = 0.5, /3 = 0.25 we again have the trapezoidal rule. Figs. 1 and 2 show that at this step size only the backward Euler method gives acceptable results, whereas those from the other methods are contaminated by parasitic oscillations and consequent unacceptable errors. The error in the backward Euler method has the same order of magnitude as the truncation error in the initial data. The step size is such that 360 steps are needed for one full revolution. By numerical experiment it was found that with the trapezoidal rule one would need approximately one thousand times as many steps to achieve the accuracy of the backward Euler method. For these reasons the backward Euler method was chosen for the stability computations with 360 steps per revolution. The finite element model of the wind energy converter GROWIAN and the main dimensions are shown in Fig. 3. GROWIAN has been erected at the mouth of the river Elbe in northern Germany and is presently undergoing final tests before operation. With a tower height of 100 m and a rotor diameter of also 100 m it is the largest of its kind. The finite element model consists of nonlinear beam and truss elements for the tower, the rotor and the parallel wire strands and the linear beam elements to model the stiffness of the foundation. To fulfill vibrational requirements, i.e., lower eigenfrequencies not to be resonant with multiples of rotor shaft revolutions at service speeds, each of the wire strands are given an initial tension force of 5300 kN. To reduce the vibrational input from the rotor the convertor is constructed with a teetering hub. This additional degree of freedom increases the risk of dynamic instability.
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
199
wind
generated height diameter mass
f e-model
model
power
of the large wind energy
3000
kW
100 m
of rotor
of nacelle
rotating
Fig. 3. Finite element
electric
of tower
100 m and
rotor
mass 160 degrees
converter
370 Mg 150 Mg
of freedom-13
elements
GROWIAN.
Equal attention as to the time-stepping algorithm has to be paid to the initial conditions. Transients from inaccurate initial data could propagate through a whole period (revolution) and hence invalidate the stability analysis. To avoid this, the initial conditions are taken from a quasi-static analysis using the same loads as for the dynamic case but now including explicitly the forces due to rotating masses. The deformations calculated in this manner correspond to the initial displacements of a stable rotating system at that instant (dynamic load step zero). The corresponding initial velocities are also easily computed under the assumption of a steadily rotating system. In general it will not always be easy to obtain initial conditions. A systematic approach would be to exploit the fixed point property of the Poincare map and employ some suitable fixed-point algorithm. The structure with applied initial conditions is given in Fig. 4. These initial conditions are now perturbed by the mode shapes of the lowest eigenfrequencies as described in Section 3. The 15 lowest eigenfrequencies in the state of rotation are shown in Fig. 5. The six vibration modes due to the wire strands were not deemed to be critical for the global stability, since they may be damped very easily, and all but one were therefore excluded from the stability investigation. The one retained (mode #lo: 0.95 Hz) was thought to be most critical in conjunction with the first bending mode of the tower (mode #2: 0.47 Hz). In Table 1 the moduli of the eigenvalues of the monodromy matrix are given for an increasing number of eigenmodes used in the perturbation. Obviously mode f4 excites a strong dynamic instability. Also the inclusion of the wire strand mode (the fifth mode in the stability analysis) raises the modulus of the eigenvalue corresponding with the tower bending slightly, but the eigenvalue stays inside the unit circle. The instability excited by mode R4 was also detected by parallel computations [15], but along the lines of the more conventional approach as described in Section 2. This instability means physically that the rotational speed of the rotor will not be stable. The rotor is driven by the wind forces and the power is extracted by the generator. In the finite element model used the generator is modeled very crudely as a simple damper of the rotational degree of freedom. In reality there is an elaborate feedback control system to keep the rotor at a constant speed. Computations of the control system show, however, that the rotation will be stable. Therefore, the mode f4 was excluded from the stability analysis, as well as the wire strand mode #lo, which did not
200
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
displacement scale 20 :l
Fig. 4. Structure
with applied
initial conditions.
contribute much. The moduli of the eigenvalues of the monodromy matrix for this perturbation pattern are given in Table 2 for an increasing sequence of modes. The eigenvalues are all strictly well inside the unit circle and the system is stable. In both cases the quick convergence of eigenvalues is remarkable.
0.88Hz P
1.71Hz
Fig. 5. Mode shapes
of the 15 lowest eigenfrequencies.
\
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
201
Table 1 Moduli of eigenvalues of I&,, (rotational and wire strand mode included) k (m =O)
1
2
3
4
5.
6
7
8
9
10
0.327
0.457 0.803
0.469 0.826 0.762
1.961 0.882 0.779 0.206
1.961 0.907 0.777 0.207 0.253
1.961 0.907 0.779 0.207 0.287 0.287
1.961 0.907 0.779 0.207 0.276 0.276 0.070
1.961 0.908 0.779 0.207 0.276 0.238 0.238 0.026
1.962 0.906 0.782 0.207 0.199 0.199 0.195 0.195 0.011
1.961 0.906 0.784 0.207 0.205 0.205 0.202 0.202 0.201 0.201
Table 2 Moduli of eigenvalues of I&,,, (rotational and wire strand mode excluded) k (m =0)
1
2
3
4
5
6
7
8
0.327
0.457 0.803
0.469 0.826 0.762
0.469 0.825 0.761 0.322
0.469 0.826 0.761 0.296 0.074
0.469 0.826 0.761 0.241 0.241 0.023
0.469 0.826 0.762 0.241 0.241 0.027 0.027
0.469 0.826 0.762 0.254 0.254 0.197 0.197 0.012
5. Conclusion The method presented for the dynamic stability analysis differs from the more conventional approach in that the two steps ‘linearisation’ and ‘time-discretisation’ are interchanged. In this form it is better suited for the stability investigation of large systems, along with some other advantages as described in Section 3. This is particularly true if the linearisation is only performed in a subspace spanned by a well-chosen set of basis vectors. Here, the vibration models of the lowest eigenfrequencies were used. This can be regarded as a Galerkin approximation of the original system. The choice of basis vectors for the subspace may be influenced, in that modes which are either not important or not well represented in the model may be excluded from the stability analysis. The method was applied to study the stability of the large wind energy convertor GROWIAN. The eigenvalues of the monodromy matrix numerically show a very quick convergence for an increasing sequence of subspaces.
202
H.G. Matthies, Chr. Nath, Dynamic stability of periodic solutions of large systems
References [l] J.P. LaSalle, The Stability of Dynamical Systems (SIAM, Philadelphia, PA, 1976). [2] M.W. Hirsch and S. Smale, Differential Equations, Dynamical Systems and Linear Algebra (Academic Press, New York, 1974). [3] N. Rouche and J. Mawhin, Ordinary Differential Equations-Stability and Periodic Solutions (Pitman, London, 1980) (translated from the French 1973 edition). [4] H. Ziegler, Principles of Structural Stability (Birkhauser, Base], 2nd ed., 1977). [S] H. Leipholz, Stability of Elastic Systems (Sijthoff & Noordhoff, Alphen aan den Rijn, 1980). [6] P. Friedmann, C.E. Hammond and Tze-Hsin Woo, Efficient numerical treatment of periodic systems with application to stability problems, Internat. J. Numer. Meths. Engrg. 11 (1977) 1117-1136. [7] H. Brauchli, On the norm dependence of the concept of stability, in: P. Germain and B. Nayroles, eds., Applications of Methods of Functional Analysis to Problems in Mechanics, Lecture Notes in Mathematics 503 (Springer, Berlin, 1976). [8] L. Collatz, Differettialgleichungen (Teubner, Stuttgart, 1970). [9] G. Dahlquist and A. Bjiirck, Numerical Methods (Prentice Hall, Englewood Cliffs, NJ, 1974). [lo] T.J.R. Hughes, Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics, Comput. & Structures 6 (1976) 313-324. [ll] Z. Zlatev and P.G. Thomsen, Application of backward differentiation methods to the finite element solution of time-dependent problems, Internat. J. Numer. Meths. Engrg. 14 (1979) 1051-1061. [12] A. Prothero and A. Robinson, On the stability and accuracy of one-step methods for solving stiff systems of ordinary differential equations, Math. Comput. 28 (125) (1974) 145-162. [13] O.C. Zienkiewicz, The Finite Element Method (McGraw Hill, London, 3rd ed., 1977). [14] K.J. Bathe, Finite Element Procedures in Engineering Analysis (Prentice-Hall, Englewood Cliffs, NJ, 1982). [15] K. Kehl, W. Keim, F. KieBling and A. Rippl, Zur Dynamik von grol3en Windkraftanlagen, VDI-Berichte 381 (1980) 113-120. [16] J.H. Argyris, K.A. Braun, B. KirchgaBner and R. Walther, Statische und dynamische Untersuchungen an einem Windrotormodell, ISD-Bericht No. 259, Universitat Stuttgart, 1979.