ARTICLE IN PRESS

JOURNAL OF
SOUND AND
VIBRATION
Journal of Sound and Vibration 271 (2004) 1015–1038
www.elsevier.com/locate/jsvi

Assessing the local stability of periodic motions for

large multibody non-linear systems using proper
orthogonal decomposition
G. Quaranta*, P. Mantegazza, P. Masarati
Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, via La Masa 34, 20156 Milano, Italy
Received 25 October 2002; accepted 14 March 2003

Abstract

The eigenvalues of the monodromy matrix, known as Floquet characteristic multipliers, are used to study
the local stability of periodic motions of a non-linear system of differential-algebraic equations (DAE).
When the size of the underlying system is large, the cost of computing the monodromy matrix and its
eigenvalues may be too high. In addition, for non-minimal set equations, such as those of a DAE system,
there is a certain number of spurious eigenvalues associated with the algebraic constraint equations, which
are meaningless for the assessment of the stability of motions. An approach to extract the dominant
eigenvalues of the transition matrix without explicitly computing it is presented. The selection of the
eigenvalues is based on a proper orthogonal decomposition (POD), which extracts the minimal set of
dominant local modes of the transient dynamics on an ‘‘energy’’ contents basis. To make the procedure
applicable to both numerical and experimental tests, a unifying experimental philosophy is pursued for the
analysis of complex multidisciplinary multibody models. Some applications are outlined, and comparisons
with other techniques are presented to demonstrate the accuracy of the proposed procedure.
r 2003 Elsevier Ltd. All rights reserved.

1. Introduction

Non-linear systems may possess a large variety of periodic motions, denominated limit cycles
[1]. The local stability of these periodic orbits is assessed by computing the rate of attraction of
their transients. These transients may be caused by imposing initial conditions which do not
belong to the periodic orbit, or by perturbations of a steady state. Their behaviour can be
investigated using Floquet’s theory, originally developed for the analysis of the solution of linear

*Corresponding author. Tel.: +39-02-2399-8387; fax: +39-02-2399-8334.

E-mail address: [email protected] (G. Quaranta).

0022-460X/$ - see front matter r 2003 Elsevier Ltd. All rights reserved.
doi:10.1016/j.jsv.2003.03.004
 ARTICLE IN PRESS

1016 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

ordinary differential equations with periodic coefficients [1–3]. The stability of the system can be
inferred from the spectral radius of the monodromy matrix, which is the transition matrix that
relates to two periodic solution state vectors separated by a time period. In classical applications
of this method, the monodromy matrix is computed first; then its eigenvalues are evaluated. There
are different methods to compute the monodromy matrix [4], but all of them are somewhat
impractical to pursue for large numerical models, so this approach has been limited to systems
with relatively few states.
Another convenient way to study the stability of periodic orbits is based on the construction of
stroboscopic Poincare! maps. The idea is to reduce the study of continuous time systems to the
analysis of associated discrete time systems. A Poincare! map can be constructed by defining a
cross-section surface in the multidimensional phase space and measuring the transition of the
system orbits through this surface [5]. In this way a local description of the transients is obtained.
Of course, the periodic orbit in the phase space becomes a fixed point for the map. It can be shown
that the eigenvalues of the Poincare! map are perfectly equivalent to the Floquet characteristic
multipliers [5]. The stability condition, analogous to that of maps, states that the magnitude of the
characteristic multipliers must be less than unity.
As a consequence of the increase in computer power, complex dynamical systems can be
simulated using large-scale models. A typical example is the design of complex deformable
aeroservomechanical systems that require sophisticated analyses, which are often conducted by
means of general-purpose modelling codes based on a multibody/multidisciplinary methodology
[6]. By means of these codes the designer can progress from simple rigid body models up to fully
detailed ones, featuring accurate and fully non-linear description of constraints, deformable
elements, servohydraulic circuits, aerodynamic forces and control system components [7]. These
general-purpose modelling tools allow the investigation of a wide range of interactional
phenomena between different subsystems. The size of the resulting model, pushed by the
increasing demand for details in deformable bodies modelling, may soon become very large. The
computation of the monodromy matrix and of its eigenvalues for these problems can be extremely
demanding, both in terms of computational power and memory requirements.
An empirical method to reconstruct a local Jacobian of the Poincare! map from data obtained
either by experiments or by simulations has been presented by Murphy et al. [8]. Despite being
very interesting, this method, based on a least-squares identification of the Jacobian matrix, might
become quickly unmanageable as the number of degrees-of-freedom used to represent the system
increase.
Among the attempts to simplify the estimate of the stability boundaries for large-scale models,
the work of Matthies and Nath [9] deserves a special mention; they conducted the stability
analysis of an ordinary differential equation (ODE) system, resulting from the discretization of a
partial differential equation (PDE) finite element model, by means of a coarser discretization,
trying to reduce the effective system dimensions. Subramanian and Gaonkar proposed a parallel
algorithm to speed up the research of the Floquet multipliers for large models, achieving only a
limited overall computational time reduction [10]. A more computationally efficient approach has
been proposed by Bauchau and Nikishkov [11,12]. It allows the evaluation of a limited number of
eigenvalues of the monodromy matrix with the largest modulus using Arnoldi’s iterative
algorithm, without requiring the explicit computation of the matrix. Even though this method
correctly assesses stability boundaries when the system is unstable, it may give misleading
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1017

information in evaluating the stability margin of periodic orbits since, for models based on DAE,
spurious eigenvalues associated with algebraic constraints may arise. Furthermore, it cannot be
applied to experimental data, as opposed to the approaches based on the explicit reconstruction of
the Jacobian matrix.
The huge amount of information generated by numerical simulations, represented by the state
vectors at each integration time step, must be synthesized to obtain few quantities of interest. In
fact, providing a viable and efficient approach to condense the available data, while maintaining
the capability of reliably finding the stability properties of periodic orbits, would make the
technique presented in Ref. [8] feasible also for large systems, with the additional advantage of
having at hand a unifying approach for numerical analyses and experimental tests. A powerful
and elegant method of data analysis to obtain a low-dimensional approximate description of high-
dimensional processes is the proper orthogonal decomposition (POD), also known as Karhunen–
Lo!eve decomposition. There are many applications of POD. It has been extensively used in fluid
mechanics, to study turbulent flows and to reduce the degrees-of-freedom of the numerical models
[13], or in the reduction of complex viscous transonic aerodynamic fields inside turbomachinery
[14]. Recently, it has been applied in dynamic studies of structural vibrations [15,16], and for
damage detection [17]. Data analysis using POD is conducted to extract a set of basis functions,
called proper orthogonal modes (POMs), from experimental data or from numerical simulations,
to be subsequently used in a Galerkin projection yielding low-dimensional dynamical models.
These functions are optimal, in the sense that fewer POD modes are needed to account for the
same amount of ‘‘signal energy’’, compared to any other orthogonal basis [13]. Thus, the POMs
are a minimal set of output signals that can be used to identify the dominant eigenvalues of the
Jacobian matrix of the Poincare! map. After evaluating the multi-dimensional Poincare! maps
associated with these signals, pieces of information regarding the dominant eigenvalues can be
extracted by means of standard system identification procedures.
This paper is structured as follows: Section 2 briefly reviews the basics of the stability theory for
periodic orbits and how it relates to Poincare! maps. Section 3 presents POD, and how it is applied
to the extraction of the leading information regarding the stability of the motions under
investigation. Section 4 summarizes all the steps required to apply the proposed method to assess
the stability of motions. In Section 5 a brief introduction to DAE systems and how they are used
to build multibody models is presented. By means of a simple example it is shown how spurious
eigenvalues may arise. Finally, Section 6 presents two numerical examples to validate the
proposed method and illustrate its accuracy. The first one addresses the stability of the classical
problem of a beam subjected to a harmonic compressive follower force [18]. The second one
shows the application of the presented technique to a multibody/multidisciplinary model, to study
the ground resonance of a helicopter [19], a typical mechanical periodic instability.

2. Periodic orbits and their stability

Consider an autonomous N-dimensional dynamic system

x’ ¼ fðxÞ; ð1Þ
 ARTICLE IN PRESS

1018 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

* ¼ xðt
which possesses a periodic solution xðtÞ * þ TÞ; characterized in the phase space by the orbit
g: The stability of g can be assessed by studying the way g-neighbouring trajectories behave. By
linearizing the differential equation about g one obtains
n’ ¼ DfðxðtÞÞn;
* ð2Þ

*
where DfðxðtÞÞ is a periodic matrix. The monodromy matrix is the transition matrix computed on
a period T; so it is defined as
nðTÞ ¼ UðT; 0Þnð0Þ: ð3Þ

Rearranging Eqs. (2) and (3), the transition matrix can be written as the solution of the following
differential problem (see Seydel [4]):
’ ¼ DfðxðtÞÞU;
U * Uð0; 0Þ ¼ I: ð4Þ

Resorting to Floquet’s theory [3,5], it can be shown how any fundamental solution matrix of this
periodic system has the form
Uðt; 0Þ ¼ ZðtÞeBt ð5Þ

where Z is a T-periodic matrix. Since Uð0; 0Þ ¼ Zð0Þ ¼ I; the monodromy matrix is

UT ¼ UðT; 0Þ ¼ eBT : ð6Þ

The eigenvalues of this matrix are the characteristic multipliers of the periodic orbit g: The
multiplier associated with perturbations along g in the phase space is always equal to one, so the
modules of the remaining ðN  1Þ multipliers orthogonal to g determine the stability properties of
the periodic orbit g:
The study of the stability can also be accomplished by reducing the continuous system (1) to the
ðN  1Þ-dimensional associated discrete map ynþ1 ¼ Fðyn Þ; called Poincare! map, where y is an
ðN  1Þ state vector. This map is obtained by sampling the ðN  1Þ state variables, triggered by the
orbit crossing an hyperplane transverse to the vector field fðxÞ; usually called a Poincar!e section
(for further details see Ref. [5]). By definition a periodic orbit is one that revisits a point in the
phase space once every period. So, the stability of a periodic orbit is reflected by the properties of
the fixed points for the ðN  1Þ-dimensional map. In fact, the eigenvalues of the linearized
Poincare! map DFðyg Þ correspond to the ðN  1Þ characteristic multipliers which determine the
stability properties of the periodic orbit.
Non-autonomous periodically forced systems, x’ ¼ fðx; tÞ; may be treated in the same way by
simply noting that they can be viewed as autonomous, at the expense of an increase of their
dimension by one, i.e., including time as a state variable
x’ ¼ fðx; yÞ;
y’ ¼ o;

where o is the circular frequency of the periodic forcing signal. In this case the Poincare! section
will be simply defined by any y ¼ y0 ; with y0 A½0; 2p:
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1019

2.1. Evaluation of the monodromy matrix

The previous paragraph provides a sequence of steps required to estimate the stability of a
periodic motion. First of all it is essential to evaluate the monodromy matrix UT ; or equivalently
the Jacobian of the associated Poincare! map. As soon as this matrix is available, the stability of a
periodic trajectory may be assessed by means of its eigenvalues. A fixed point for a discrete system
is stable if and only if the spectral radius of the Jacobian, computed around the fixed point, is less
or equal to unity. The classical method to generate the monodromy matrix is based on Eq. (4):
integrating these N 2 differential equations for 0ptpT yields the monodromy matrix but, since
*
the Jacobian DfðxðtÞÞ varies along the orbit g and hence with time, it is necessary for the general
non-linear case to integrate until t ¼ T the ðN þ N 2 Þ-dimension initial value problem
*
x’ ¼ fðxÞ; xð0Þ ¼ xð0Þ;
U’ ¼ DfðxðtÞÞU;
* Uð0Þ ¼ I:

In case of linear systems only the last N 2 equations are required. Numerous methods have been
proposed to reduce the overwhelming computational effort required to solve this problem. For an
introductory description the reader is referred to the methods suggested by Seydel in Ref. [4]. In
the linear case the problem can be simplified by using the method proposed by Friedmann et al.
[20]. In any case, when the dimension N of the problem is large, the calculation of UT and all its
eigenvalues may be too costly to be effectively used in most applications. Moreover, in practical
analyses UT is seldom known analytically, so it must be obtained by a numerical linearization; as
a consequence, the quality of its approximation must often be verified by means of repeated
computations with different numerical parameters to assess the discretization errors. Finally, in
any case most of the eigenvalues of UT do not give any useful insight for the investigation of
stability problems. Thus, to simplify the stability analysis, it is necessary to find a method capable
of estimating and selecting the dominant eigenvalues of UT ; which characterize the phenomenon
under investigation, while reducing the computational burden as much as possible.

3. Extracting dominant transients from simulation

The approach followed here stems from the method proposed by Murphy et al. [8] and is
strictly connected to the idea of numerical simulations viewed as ‘‘numerical experiments’’, whose
results are analyzed by techniques that are suitable for experimental data as well. In essence, the
monodromy matrix allows one to write the following expression:
xðkþ1Þ ¼ UT xðkÞ ; ð7Þ

where the index k is related to samples one period T apart. In the recursive equation written above
the matrix UT must be evaluated starting from appropriate vectors xðkÞ ; that are computed
numerically by integrating Eq. (1) after applying a suitable perturbation. The perturbation must
guarantee an adequate linearization while keeping the solution within the attraction basin of the
periodic orbit under investigation. The different approaches for the determination of UT can be
related to different interpretations of Eq. (7).
 ARTICLE IN PRESS

1020 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

On one hand, it may be viewed as an auto-regressive (AR) discrete time model

Xðkþ1Þ ¼ UT XðkÞ ; ð8Þ
where
Xðkþ1Þ ¼ ½yxðkþ1Þ xðkÞ y;
XðkÞ ¼ ½yxðkÞ xðk1Þ y;
which is a generalization of the least-squares linear fitting of the system response to an initial
perturbation, observed either in experiments or simulations, to compute an approximation of the
linearized Poincare! map presented in Ref. [8].
On the other hand, it can also be viewed as a power iteration, in which case the components
ratio xðkþ1Þ =xðkÞ is known to converge to the largest eigenvalue of UT ; if the initial vector is not
defective in the direction of the related eigenvector [21]. The power iteration (7), in principle,
suffices to determine the stability boundaries. However, since the largest eigenvalue may not be
well separated from neighbouring eigenvalues, the method must be extended to look the
eigenvalues in the vicinity of the largest modulus one by using an iterative subspace method [22]
where, at each step, the algorithm chooses the best initial conditions for Eq. (7).
Bauchau and Nikishkov propose an implicit matrix method [11] exploiting the properties of the
Arnoldi’s algorithm: a subspace method, which needs only a matrix-vector multiplication to extract
the highest modulus eigenvalues of large and sparse matrices, and, at the same time, does not
require the solution of the transpose problem, as other subspace methods for non-symmetric
problems do. Therefore, the eigenvalues of the transition matrix Uðt; 0Þ can be obtained by
performing a computation of the response after a certain time interval to an initial condition vector
chosen by the algorithm. If a sufficiently small time interval is chosen, it is possible to evaluate the
eigenvalues of the linearized system near an equilibrium condition. By choosing the period T as time
interval, the eigenvalues of the monodromy matrix can be computed (if a periodic response exists).
Even though this method gives the correct estimate of the stability boundary, it can be misleading in
giving information about the stability margin of DAE systems. In fact, as it will be seen in Section 5,
there are always unit modulus eigenvalues associated with the algebraic constraints.
A way must be devised to select a reduced set of signals synthesizing the information related to the
dominant transient, to be able to apply any method for the evaluation of stability properties. These
signals can be selected by resorting to a technique capable of extracting spatial coherence in an
oscillating system, when the time history of its state variables is known from either numerical
simulations or from the output of several sensors measured in real-life experiments. This method is
represented by the POD used for the analysis of multi-dimensional data. It provides a way to find the
best approximating subspace to a given set of data in a least-squares sense. The POD allows one to
obtain a modal decomposition that is only data dependent and does not assume any prior knowledge of
the system. Consequently, it perfectly fits the analysis needs of multidisciplinary models. A reduced
order model of a dynamical system can be obtained by simply using a Galerkin projection procedure
where the POD modes are used as basis functions. The same ideas can be applied here for the generation
of a small subset of significant signals to be used for the identification of the leading eigenvalues
characterizing system behaviour and stability. For the analysis of mechanical systems, it must be further
noticed that these base POMs have an interesting interpretation, since they can be viewed as the result of
a least-squares error optimization for a linear representation of the non-linear normal modes [15].
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1021

3.1. Computation of the proper orthogonal decomposition

Consider a system where all the N state variables are measured at n time steps. Their time
averages are usually subtracted from the signals and data are arranged in a ðN  nÞ matrix

X ¼ ½xð1Þ ; xð2Þ ; y; xðnÞ :

An approximation of the system dynamics is obtained by projecting the original N-dimensional

state space onto an m-dimensional subspace S: The main purpose of POD is to find a projection
operator Q mapping RN onto S; which minimizes the Euclidean distance of the sampled points
from the m-dimensional hyperplane
X
n
HðQÞ ¼ jjxðiÞ  QxðiÞ jj:
i¼1

It can be shown (see Ref. [13]) that, given l1 Xl2 X?Xln ; the eigenvalues of the data correlation
matrix R ¼ XXT ; the optimal m-dimensional projection operator is represented by the ðm  NÞ
matrix whose rows are the first m eigenvectors of R: To obtain these eigenvectors, the singular
value decomposition (SVD) of matrix XT can be computed:

XT ¼ URVT ;

where R is the diagonal matrix of the singular values si ; and the columns of V are the POMs. By
sorting the singular values si of X in a descending order, it can be shown that the matrix
XTm ¼ Um Rm VTm ; where Um and Vm are the rectangular matrices obtained by retaining the first m
columns, and Rm is the ðm  mÞ principal minor of R; is the closest rank m matrix to XT in the
Frobenius norm. To choose the dimension m of the approximate subspace that will contain all the
significant information one must look at the singular values, since their value expresses the ‘‘signal
energy’’ related to the associated POM. The following examples will show that these values reach
a constant plateau, usually called ‘‘noise floor’’, that characterizes the modes which do not contain
any significant information. In fact, the use of some form of SVD as a tool to compute the order
of a model is a common practice in system identification as well.
If Nbn; then it is more efficient to use the so-called method of snapshots. It consists of first
computing the matrix U as the eigenvectors matrix of XT X: Once U is known, since
UT XT ¼ RVT ;
it is clear that the norm of the first m rows of RVT are the singular values, and that these rows,
after normalization, are the POMs.

4. Assessing local stability

The proposed method to assess local stability is summarized in the following steps:
(1) Find a periodic orbit. This is usually done by integrating the system for a few periods, possibly
with some artificial damping to accelerate convergence to the periodic solution.
 ARTICLE IN PRESS

1022 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

(2) Perturb the system and analyze the induced transient to obtain the elements of Eq. (8). When
dealing with systems with a large number of degrees-of-freedom the choice of an effective set
of perturbations may not be trivial.
(3) Choose a Poincare! section in the phase space and sample the related time histories at periods
T to obtain the corresponding Poincare! map.
(4) Analyze the sampled signals, using the snapshot method to compute the POD. The dominant
POMs will be selected, based on the relative magnitude of the singular values.
(5) Using the time histories of the POM, obtain a minimal set of signals to be processed. By
means of any standard system identification method for multiple output systems [23], e.g. AR
based on the least-squares method, an approximate Jacobian matrix is determined, which
possesses the dominant eigenvalues of the phenomenon under investigation. The more
sophisticated algorithm N4SID [24] can be used for the analysis of experimental data, since it
also contains an optimal treatment of noisy signals.
(6) Extract the Floquet multipliers by using any numerical method for the determination of
eigenvalues.

Often a prior knowledge of the physical phenomena under investigation will help in choosing
the type and level of perturbations. In the extensive experiences of the authors, some of which are
reported in the following test cases, no problems were encountered for the selection of correct
perturbations. An analysis of the effects of different perturbations and a possible method to
choose appropriate signals in complex cases is presented in the test section of this paper. Steps 4
and 5 can be applied recursively to verify the convergence of the solution by increasing the
number of samples evaluated. If the problem size is small, step 4 is not necessary. Instead, it
becomes nearly mandatory when the dimension of the system is large or when it is necessary to
distinguish between physical and spurious eigenvalues, as explained in the following section.

5. Multibody systems and related eigenvalues

The large-scale models that are investigated in this work fall into the class of multidisciplinary
multibody systems [6]. They are represented by special DAE systems. To clarify the peculiarities
of these models, a brief presentation of the basic formulation is needed.
The motion of a system of rigid/deformable bodies can be described using the principles of
classical mechanics. They may lead to a redundant co-ordinate set approach, in which the
constraints between inertial bodies are explicitly expressed, or to a minimal co-ordinate set
approach where all the constraints are eliminated. The numerical reduction to minimal set
equations usually results in methods which are computationally less effective if complicated
mechanisms are investigated. The redundant approach instead can be used to easily and
automatically generate the equations that model complex systems by means of the multibody
formalism. A multibody system can be defined as a collection of bodies in arbitrary motion with
respect to each other. These bodies can be connected by algebraic constraints and flexible
elements; they can also be subjected to interactional forces that represent other multidisciplinary
elements. The advantage of this technique is a simpler formulation, resulting in an easy and
versatile implementation, since it is possible to write the free-body dynamics equations for each
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1023

body in the global reference frame, with the addition of the constraint equations. One drawback is
a larger amount of required degrees-of-freedom, compared to that of the reduced co-ordinate set
formulation, which can nonetheless be managed effectively by sparse matrix numerical methods.
The basic equations of motion for the bodies can be described recurring to the Lagrangian
multipliers technique to introduce algebraic constraints. Such multipliers are related to the
reaction forces exchanged by the bodies, so they are of great significance for engineers. Since their
calculation is often required, their direct availability further enhances the overall efficiency of the
redundant co-ordinate set against minimal set formulations. The dynamics of the bodies are
written in the form of a fully implicit first order differential system by coupling the definition of
the momentum, b; and of the momenta moment, C; of each body to the force and moment
equilibrium equations [7]. Given the static moments vector S and the inertia moments matrix J;
the resulting set of equations for each body is
b ¼ mx’ þ x  S;
C ¼ S  x’ þ Jx;
b’ ¼ Fðx; x;
’ R; x; a; a’ ; y; tÞ þ VF ;
C’ þ x’  b ¼ Cðx; x;
’ R; x; a; a’ ; y; tÞ þ VC ; ð9Þ

where the forces F and couples C may depend on: the configuration represented by the position
vector x; the orientation matrix R and the velocities, translational x’ and angular x; the internal
states a that can be generally associated with the multidisciplinary fields that are simultaneously
considered; and the time t: The forces VF and VC are related to the Lagrangian multipliers that
represent the actual reaction forces and couples generated by the constraints. Kinematic
constraints, both holonomic and non-holonomic, are added as algebraic/differential equations:
Wh ðx; R; y; tÞ ¼ 0;
’ R; x; y; tÞ ¼ 0;
Wnh ðx; x; ð10Þ

The resulting system of non-linear DAE, obtained joining Eqs. (9) and (10), requires a special
treatment [25,26] due to the singularity of the algebraic equations if the problem is treated as
differential. By calling y the kinematic unknowns, z the momentum unknowns and k the algebraic
Lagrangian multipliers unknowns, the system can be cast in the following general form
Mðy; tÞ’y ¼ z;
z’ ¼ Qðy; y’ ; tÞ  GT k;
Wðy; tÞ ¼ 0: ð11Þ

In these equations, M is a configuration dependent inertia matrix, Q are arbitrary external forces
and couples and G ¼ Dy W is the Jacobian of the holonomic constraints with respect to the
kinematic unknowns. The final system is DAE of index three, meaning that three differentiations
with respect to time are required to obtain y’ as a continuous function of ðy; tÞ [25]. Non-
holonomic constraints require a slightly different treatment and result in a lower index DAE
system, whose solution is less critical. Numerous techniques have been proposed to solve this kind
of problem [25–28]; all the examples presented in the following are solved directly in DAE form
 ARTICLE IN PRESS

1024 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

resorting to a fully implicit A/L-stable, second order accurate predictor–corrector integration

scheme [7,28].
For index three DAE systems a constraint represents a ð2n  mÞ-dimensional manifold
M ¼ fðy; zÞ j WðyÞ ¼ 0; Dy WðyÞM1 ðyÞz ¼ 0g

on which the solution must lie, where n is the dimension of the y vector, and m is the number of
constraint equations WðyÞ ¼ 0: The dynamic behaviour of the system represented by Eq. (11) will
be locally dominated by the eigenvalues of the linearized vector field that lies in the tangent space
of M: To obtain this information one should ideally compute the eigenvalues of the problem
stated in terms of the minimal co-ordinate set, which represents the system of differential
equations projected on the manifold M: The evaluation of the eigenvalues directly on the DAE
systems will result in ð2n  mÞ correct eigenvalues, since m is also the number of Lagrange
multipliers, and 2m spurious eigenvalues that will not give any useful insight into the dynamics of
the mechanism. Consequently, a method to discern among the computed eigenvalues is needed.
Usually the m eigenvalues associated to the Lagrangian multipliers assume a minus infinite value,
to indicate that the constraint equations are instantaneously satisfied (infinitely fast dynamics).
The eigenvalues associated with the redundant co-ordinates instead, assume zero value, since the
constrained degrees-of-freedom have no dynamics. So, if the entire spectrum is computed, the
significant eigenvalues must be selected in a post-processing phase. However, if only a subset of
the eigenvalues is required, an algorithm to select the interesting roots is mandatory. The
problems that can be encountered while selecting the significant eigenvalues are illustrated in the
following.

5.1. Eigenvalues of a three-masses system

Consider a chain of three masses m ¼ 1 shown in Fig. 1. The first and the second mass are
connected in series through a spring of stiffness k ¼ 1 and a damper with d ¼ 0:1: The second and
the third mass are rigidly connected. The first mass is constrained to the ground so that the whole
system represents a simple oscillator by means of linear DAEs. The resulting equations read:
mx. 1 þ d x’ 1  d x’ 2 þ kx1  kx2  l1 ¼ 0;
mx. 2  d x’ 1 þ d x’ 2  kx1 þ kx2  l2 ¼ 0;
mx. 3 þ l2 ¼ 0;
x1 ¼ 0;
x2  x3 ¼ 0: ð12Þ

x1 x2 x

Fig. 1. Three-masses model layout.

 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1025

Table 1
Eigenvalues computed by means of the implicit matrix approach
L l
1.000 0.0000
1.000 0.0000
1.000 0.0000
1.000 0.0000
0.9997+j 7.0665e4 0.0250+j 0.7067
0.9997j 7.0665e4 0.0250j 0.7067
1.698e17 38614
2.600e22 49701

Table 2
First four singular values computed for the three-masses problem
s
6.0745
2.8217
2:2470e  16
9:2296e  17

These equations, cast in state space form, yield an eight-dimensional DAE system. Clearly, the
same problem can be described in a more compact form by eliminating the constraints
2mx. 2 þ d x’ 2 þ kx2 ¼ 0: ð13Þ
Eq. (13) is characterized by the two complex conjugated eigenvalues l1j2 ¼ 0:0257j 0:7067; that
contain all the useful information to be extracted from system (12).
To evaluate these eigenvalues from the DAE system, the method of the implicit matrix
proposed by Bauchau and Nikishkov [11] can be used. The eigenvalues of the linearized system
matrix l; and that of the transition matrix Uðh; 0Þ; L; are related by
1
l ¼ ln L; ð14Þ
h
where h is the time step used to compute the eigenvalues of the transition matrix, and ‘‘ln’’ is the
principal natural logarithm. Using a time step of h ¼ 0:001 for the present case leads to the eight
eigenvalues listed in Table 1.
The correct eigenvalues that represent the dynamics of the system are well approximated, but,
as can clearly be seen, they are not the largest modulus ones (it is worth recalling that the
eigenvalues on the left, related to the transition matrix, are actually computed). So, seeking the
first two eigenvalues with the Arnoldi’s method would lead to erroneous information about the
stability margin of the system.
The same problem can be solved by applying the method proposed in this paper. In this case the
step regarding the sampling process to obtain the Poincare! map is not required since the problem
is not periodic. The system is perturbed by giving an initial unit displacement for both the second
and third mass. Table 2 lists the amplitude of the first four higher singular values. It clearly shows
 ARTICLE IN PRESS

1026 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

Table 3
Eigenvalues computed by means of the POD modes
L l
0:9997 þ j 7:0665e  4 0:0250 þ j 0:7067
0:9997  j 7:0665e  4 0:0250  j 0:7067

that there are only two dominant dynamics in the system. A simple least-square fitting of the time
histories of the associated POMs gives a Jacobian matrix which possesses the two eigenvalues
listed in Table 3. No spurious eigenvalues appear.

6. Numerical tests

The numerical method presented in the previous sections is best illustrated by applying it to two
problems. The first example is a classical test case that will also be used to validate the
methodology. The second example shows the application of the method to a real multidisciplinary
case, a complex helicopter model, including flexible non-linear elements and landing gear
components, such as tyres and shock absorbers, are modelled.

6.1. Beam under compressive follower force

A convenient example, useful for the purpose of validating the methodology, is the study of the
instability of a simply supported beam under the influence of a periodic compressive follower axial
load applied at one of its ends [18]. Consider the continuous beam of length L with mass per unit
length m and uniform bending stiffness EJ; subjected to an axial compressive load of the form
PðtÞ ¼ P0 þ P1 cos oP t; shown in Fig. 2. Call s the co-ordinate along the beam axis and wðs; tÞ the
beam transverse deflection. A distributed linear viscous structural damping d is assumed. The
equations of motion can be derived by using Lagrange’s formalism. The kinetic energy associated
with the transverse beam vibration is
Z L
1
T¼ ’ tÞ2 ds:
mwðs;
0 2
The potential energy is given by the bending strain energy
Z L
1 w00 ðs; tÞ2
V¼ EJ ds:
0 2 ð1  w0 ðs; tÞ2 Þ
where ð Þ0  @=@s: Assuming as solution
X ips
wðs; tÞ ¼ qi ðtÞ sin ; i ¼ 1; 2; y; ð15Þ
i
L
which satisfies all the boundary conditions, it is possible to write the usual Lagrange’s
equations in terms of the qi degrees-of-freedom. The generalized forces corresponding to the
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1027

w(x)

x P

Fig. 2. Simply supported beam under compressive follower load.

follower force PðtÞ is

Z L qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 
@ 0 2
Qi ¼ PðtÞ 1  w ðs; tÞ  1 ds;
@qi 0

while the terms associated with the damping forces are equal to
Z L
@
Di ¼ ’ tÞ ds:
dwðs;
@qi 0
The resulting Lagrange’s equations are
 
d @T @V
 þ Di ¼ Qi ; i ¼ 1; 2; y : ð16Þ
dt @q’ i @qi
In the non-linear system of equations (16) the different modal degrees-of-freedom are coupled,
and each single equation resembles Duffing’s equation [5]. The corresponding linearized
governing equation is the following:
@4 w @2 w
mw. þ d w’ þ EJ 4 þ ðP0 þ P1 cos oP tÞ 2 ¼ 0: ð17Þ
@s @s
Using the modal expansion of Eq. (15), the PDE above is transformed in a series of Mathieu
uncoupled equations for each modal amplitude
q. i þ eq’ i þ o2i ð1  mi cos oP tÞqi ¼ 0; ð18Þ
where the following parameters are defined
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 ffi
d i2 p2 EJ i2 p2 EJ P0 P1
e¼ ; Pi ¼ ; oi ¼ 2 1 ; mi ¼ :
m L2 L m Pi ðPi  P0 Þ
It is observed that in many situations the lowest mode is dominant, so that a lumped single mode
linearized approximation represents a good simplified model, useful enough to validate stability
results.
By considering a solution made of only the first buckling mode of Eq. (15), the stability regions
in the m–O plane, O ¼ oP =2o1 being a non-dimensional excitation frequency, can be found using
the classical Hill’s infinite determinant method [29]. A simpler way to obtain the value of the
spectral radius for the different conditions is given in Ref. [5], where it is noticed that the
determinant of the Jacobian matrix of the Poincare! map associated with Eq. (18) is equal to
detðDF Þ ¼ l1 l2 ¼ e2pe=op : If the eigenvalues are complex this product is also equal to the square
 ARTICLE IN PRESS

1028 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

0.50003 1000
0.50002
500
0.50001
0.5 0
x

w
0.49999
-500
0.49998
0.49997 -1000

0.49996 -1500
0 5 10 15 20 25 30 0 5 10 15 20 25 30
t/T
(a) (b) t/T

Fig. 3. Unstable periodic motion: m ¼ 0:6 and O ¼ 0:9; - -, time history; —, sampled signal; (a) beam mid-point axial
displacement; (b) beam mid-point transversal displacement.

of the spectral radius:

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
r¼ e2pe=op : ð19Þ
In any case the eigenvalues can easily be found by first computing the ð2  2Þ monodromy matrix
following the classical approach.
The multibody model built to represent the system under investigation is made of 10 non-linear
finite volume three-node beam elements [30]. Such a model includes also the exact representation
of the two kinematic constraints, resulting in 264 states. An extremely small static load is applied
transversally at the mid-point to model the effect of a beam imperfection. In this case the
perturbation load can be easily chosen. It is known that the instability will be related to the first
mode, so by applying an impulsive transverse load at the mid-point a perturbation is obtained
which possesses a component in the direction of the dominant eigenvector.
Fig. 3 shows two different signals measured for an unstable combination of the parameter m and
O; the first one is the axial displacement at the mid-point, while the second is the transverse
deflection. It is clear how, if the wrong signals are chosen, the instability may not be tangibly
captured. So, an analysis capable of synthesizing all the dominant information from a large model
is necessary, especially when it is not clear which is the main instability mechanism. In the latter
case the choice of the perturbation signal can be made starting from the identification of the
unstable mode running the POD analysis in an unstable condition, while the system is subjected to
an impulsive excitation on all of its states. The obtained results will be used in the subsequent
tests, where the excitation load will be applied to several antinodes of the first few POMs identified
in the instability condition.
After the simulation is run, sampled data are processed to extract the singular values and the
POMs. All the sampled time histories of the system states are scaled first, to avoid misleading
results in all the cases where mixed physics signals are analyzed together (a situation that can often
happen for multidisciplinary models).
Figs. 4 and 5 compare the spectral radius computed from the single mode linear Mathieu model
with the results obtained by applying the proposed method to the non-linear multibody model, for
two different values of the parameter m: The data reported on top of the figures show the number
of POMs to which the 99.99% of the energy of the time response is associated. Notice how the
number of significant modes decreases in the instability areas where the unstable mode dominates
the time response.
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1029

N POM
5
3
1
0.125 0.25 0.5 1
(a)

1.6
1.5
1.4
1.3

1.2
1.1
1
0.9
0.8
0.125 0.25 0.5 1
(b) Ω

Fig. 4. Comparison of the spectral radius of the monodromy matrix at m ¼ 0:6 with variable O; - -, linear single mode
results; &; POD non-linear identification; (a) POD dimensions; (b) spectral radius.

7
N POM

5
3
1
0.125 0.25 0.5 1
(a)

2.2
2
1.8
1.6

1.4
1.2
1
0.8
0.125 0.25 0.5 1
(b) Ω

Fig. 5. Comparison of the spectral radius of the monodromy matrix at m ¼ 0:9 with variable O; - -, linear single mode
results; &; POD non-linear identification; (a) POD dimensions; (b) spectral radius.

Fig. 6 shows several time histories and phase space projections of the response and the
corresponding sampled signal of the beam mid-point transversal displacement, for different values
of the parameters. The signals have been scaled to make the diagrams more readable. The
agreement between the single mode model and the complete multibody non-linear analysis is very
good for stability boundaries computed around the value of OC1; while for low values of O the
non-linear model gives more damped results. This difference can be explained by looking at the
diagrams reported in Fig. 6. For O ¼ 0:33–0.842 (Figs. 6e–h), the responses of the system are
mainly due to a single mode, as the emerging single turn closed path on the phase space projection
shows. Instead, for O ¼ 0:22–0.24 (Figs. 6a–d), the responses are characterized by a visible
 ARTICLE IN PRESS

1030 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

0.6 0.6
0.4 0.4
0.2 0.2
0 0
-0.2 -0.2
-0.4 -0.4
w

w
-0.6 -0.6
-0.8 -0.8
-1 -1
-1.2 -1.2
-1.4 -1.4
0 5 10 15 20 25 30 35 40 -8 -6 -4 -2 0 2 4 6 8
.
t/T w
(a) Ω = 0.22 max= 0.2533 + j 0.8792 (b)

8 8
6 6
4 4
2 2
0 0
w

w
-2 -2
-4 -4
-6 -6
-8 -8
-10 -10
0 5 10 15 20 25 30 35 40 45 -80 -60 -40 -20 0 20 40 60 80
.
t/T w
(c) Ω = 0.24 max= 0.9182 + j 0.0917 (d)
0.4 0.4

0.2 0.2

0 0

-0.2 -0.2
w
w

-0.4 -0.4

-0.6 -0.6

-0.8 -0.8
0 10 20 30 40 50 60 70 80 -6 -4 -2 0 2 4 6
.
t/T w
(f)
(e) Ω = 0.33 max= −0.8602 + j 0.4612

0.1 0.1
0.05 0.05
0 0
-0.05 -0.05
-0.1 -0.1
-0.15 -0.15
w

-0.2 -0.2
-0.25 -0.25
-0.3 -0.3
-0.35 -0.35
-0.4 -0.4
0 5 10 15 20 25 30 35 40 -3 -2 -1 0 1 2 3
.
t/T w
(g) Ω = 0.8420 max= −0.9376 + j 0.0936 (h)

Fig. 6. Stable orbits: transversal displacement w at beam mid-point for different values of O; m ¼ 0:6: On the right
phase space projection onto the w–w’ plane; - -, time history; —, sampled signal; .; Poincar!e map sample.
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1031

30000 200000
150000
20000 100000
10000 50000
0
0

w
w

-50000
-10000 -100000
-150000
-20000
-200000
-30000 -250000
0 10 20 30 40 50 60 70 80 -2e+06-1e+06 0 1e+06 2e+06
.
t/T w
(a) Ω = 0.4650 max= 1.0919 (b)

4e+11 4e+11
3e+11 3e+11
2e+11 2e+11
1e+11 1e+11
0 0
w

w
-1e+11 -1e+11
-2e+11 -2e+11
-3e+11 -3e+11
-4e+11 -4e+11
0 20 40 60 80 100 120 140 160 -5e+12 -2e+12 0 2e+12 5e+12
.
t/T w
(c) Ω = 0.90 max= −1.4596 (d)

Fig. 7. Unstable orbits: transversal displacement at beam mid-point for different values of O; m ¼ 0:6: On the right
phase space projection onto the w–w’ plane; - -, time history; —, sampled signal; .; Poincare! map sample.

activation of the second mode, as the emerging double turn closed path on the phase space
projection shows. In the latter cases discrepancies are to be expected since the model cannot be
correctly represented by a lumped single mode beam, as in the analytical simplified case.
The Poincare! section of Fig. 7 at m ¼ 0:6 shows that there are two possible routes to instability.
The first is a simple jump bifurcation due to a turning point, while the second drives the system
through a period doubling instability. By sampling the same time history with a double period it is
possible to evaluate the stability of the new periodic solution.
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi In this case the spectral radius of the
single mode linear model will be equal to r ¼ e4pe=op ¼ 0:9634: The spectral radius obtained by
means of the simulation is equal to r ¼ 0:9411; in this case 19 POMs are necessary to cover
99.99% of the energy of the time response, consequently the single mode linear model represents
only a rough approximation of the more complex non-linear system as captured by the more
precise multibody simulation.
Table 4 presents the first three computed eigenvalues for different numbers of time samples
processed. It shows how the first one is correctly approximated with 10 samples only, which
means that it is sufficient to run a 10-period simulation to correctly identify the dominant l: For
the same problem Bauchau and Nikishkov [12] declare the need for 20 simulation steps to
correctly compute the first eigenvalue with the Arnoldi’s implicit matrix approach. Considering
that, for all the cases when the system is stable, the first significant eigenvalue is not the one with
the highest modulus (in this specific case it is the eighth), and that in both cases the most
 ARTICLE IN PRESS

1032 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

Table 4
Convergence of eigenvalues for increasing number of samples: m ¼ 0:6; O ¼ 0:8
n First Eig. Second Eig. Third Eig.
Value Error % Value Error % Value Error %
10 Real 0.8574 0.27 0.5943 2.47 0.0393 22.43
Imag 0.4627 0.22 0.3485 0.91 0.2133 45.20

30 Real 0.8597 0.00 0.5858 1.00 0.0347 8.10

Imag 0.4616 0.02 0.3517 0.00 0.1865 26.96

60 Real 0.8595 0.02 0.5844 0.76 0.0342 6.54

Imag 0.4616 0.02 0.3516 0.03 0.1754 19.40

100 Real 0.8597 0.00 0.5800 0.00 0.0321 0.00

Imag 0.4617 0.00 0.3517 0.00 0.1469 0.00

1
1

0.5
0.5
w
w

0 0

-0.5
-0.5

-1
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
(a) x (b) x

Fig. 8. Computed POM forms for two different parameter sets; (a) m ¼ 0:6 O ¼ 0:3; (b) m ¼ 0:9 O ¼ 0:75; POMs: &;
1; 3; 2; W; 3; }; 4; \; 5; ; 6; m; 7.

computationally intensive phase is the system simulation, the presented method seems to be more
effective for these specific multibody applications. Anyway, the efficiency strongly depends on the
size of the subspace of unit value eigenvectors related to the number of constraints present in the
model.
Fig. 8 shows the POMs computed in two different conditions. The first POM belonging to the
‘‘noise floor’’ is plotted. Its scattered form shows why these modes are negligible. Other
perturbations have been tested, all giving correct results, except the case of an axial perturbation
applied at the mid-point, because it is represented by a vector tangential to the periodic orbit.
A final test has been conducted using the POM forms computed at m ¼ 0:6 and O ¼ 0:3
(Fig. 8(a)) as a basis to obtain the signals subsequently used for the identifications in different
operating conditions. This allows one to save all the operations necessary for the determination of
POM forms at the end of each simulation, which may be, for large models, a computationally
intensive phase. While the optimality of the POD modal base is lost, for this simple case the small
differences do not affect the quality of the results in terms of values of the eigenvalues. These
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1033

results cannot be taken as a validation test of a general rule, since there are no great differences
between the modal bases computed in different test conditions for this particular problem (Fig. 8).
Anyway, it must be noticed that the possibility to use previously computed POM bases, or to
partially re-use them, has already been exploited in a completely different field, such as the
determination of reduced order model for transonic aerodynamic loads [31]. Further
investigations are required to assess the potential of this technique.

6.2. Ground resonance of a helicopter

This example has been chosen to show how this technique can be effectively applied to large
and complex models. A brief analysis of the problem will be presented; a more extensive study on
this particular phenomenon by means of this novel methodology will be the subject of a future
work. The ground resonance is a mechanical instability that may occur when a helicopter stands
on the ground [19,32]. The instability is caused by the interaction between the lead-lag in-plane
movement of the rotor blades and the dynamics of the rotor supports, in this specific case
represented by the aircraft body and its landing gear. In all the articulated helicopter rotors the
instability is avoided by using lead-lag dampers on the rotor blades. In this case the periodicity is a
result of the rotational movement imposed on the rotor mast, consequently the problem can be
treated again as an assessment of stability of a periodically forced non-autonomous system.
The investigation is usually tackled by means of reduced order models in which the structural
properties are represented by lumped parameters. The aerodynamic forces of the rotor are usually
neglected because they show a limited dependency on the lead-lag motion of the blade. However,
this last simplification may not be true since the helicopter rolling and pitching causes a change in
blade pitch and the presence of pitch-lag and pitch-flap couplings could affect the rotor
aerodynamic loads. To avoid any undue simplification, a detailed multibody model rotorcraft has
been built to analyze the ground resonance [33]. By means of the mentioned multibody/
multidisciplinary formulation it is possible to represent all the non-linear effects due to large
displacements and rotations, aerodynamic loads, non-linear damping and constitutive effects that
are usually neglected. The model is representative of a generic medium weight helicopter, and is
composed of two main parts: the rotor and the landing gear. For a detailed description of both the
reader is referred to Ref. [33]. A short description will be reported here to give an idea of the
complexity of the multibody system under investigation. The essential properties of the rotor are
detailed in Table 5. The blades are modelled by four parabolic beam elements. The aerodynamic
forces are obtained resorting to the classical strip theory, with tabulated aerodynamic coefficients,
which account for unsteady flow, uniform inflow and radial flow drag. A dynamic inflow model is
considered as well. The blades are attached to the hub by means of elastomeric bearings and a
viscoelastic damper. Flexible pitch links and a complete exact kinematic representation of the
swashplate are added. The airframe in this case is a single rigid body, even though more complex
flexible superelements can easily be added. The tail rotor is instead modelled by a single
concentrated force. The three legs of the landing gear are basically modelled by the structural
cylinder, the non-linear oleopneumatic shock absorber, the fork, the wheel and the attachment
braces, all connected by exact kinematic constraints. The resulting model requires 795 degrees-of-
freedom. Analyses made at the nominal angular velocity of 25 rad=s; at different values of the
lead-lag damping coefficient are presented. In this case 25 POMs are necessary to obtain a basis
 ARTICLE IN PRESS

1034 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

Table 5
Helicopter and rotor properties
Number of blades, nb 5
Radius, R 8:0 m
Cut-out, R0 1:66 m
Chord, c 0:5 m
Angular velocity, O 25:0 rad=s
Hinge offset, e 0:2 m
Root pitch stiffness, Ky 10:0 N=rad
Root flap and lag stiffness, Kb ; Kx 100:0 N=rad
Pitch link stiffness, Mc=y 1:6e5 N m=rad

Blade mass per unit length, m 10:0 Kg=m

Blade inertia per unit length, Jp 1:28 Kg m
Blade twist, y 12:0 deg
Twist stiffness, GJ 1:8e6 N m2
Beam stiffness, EJy 1:5e5 N m2
Chord stiffness, EJz 1:0e7 N m2

Overall mass, M 9500:0 kg

Overall length, L 20:0 m
Roll Inertia, Ix 3000:0 kg m2
Pitch and Yaw Inertia, Iy ; Iz 15000:0 kg m2

25

20
Singular value

15

10

0
0 5 10 15 20 25
Index of singular value

Fig. 9. Singular values for the helicopter model; lead-lag damper coefficient d ¼ 7e3 ½N m=rad=s:

which covers 99.99% of the energy. Looking at the singular values associated with each POM,
shown in Fig. 9, it is clear that the periodic response is dominated by the first 12 eigenvalues.
The periodic orbit is perturbed by means of an impulsive load applied at the tip of blade 1.
Fig. 10 shows the identified dominant multipliers. The same multipliers can be transformed by
means of Eq. (14) to obtain the eigenvalues of the linearized problem, shown in Fig. 11.
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1035

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-1 - 0.5 0 0.5 1

Fig. 10. Floquet multipliers for the helicopter for different lead-lag damper coefficient d values; ; 8:0e3 N m=rad=s;
; 7:5e3 N m=rad=s; ; 7:0e3 N m=rad=s; &; 6:8 e3 N m=rad=s; }; 6:5e3 N m=rad=s:

0.2 1.4
0 1.2
-0.2
-0.4 1
-0.6 0.8
-0.8 0.6
-1
0.4
-1.2
-1.4 0.2
-1.6 0
6.4e+3 6.8e+3 7.2e+3 7.6e+3 8e+3 6.4e+3 6.8e+3 7.2e+3 7.6e+3 8e+3
(a) lead-lag damping coefficient [Nm/rad/s] (b) lead-lag damper coefficient [Nm/rad/s]

Fig. 11. System eigenvalues obtained transforming the Floquet multipliers; (a) real part; (b) frequency.

The path of the estimated multipliers shows that the system is going through a Naimark–Saker
bifurcation, the equivalent for maps of the Hopf bifurcation [5], crossing the value of
6:9e3 N m=rad=s for the lag damper coefficient. The instability is caused by an interaction
between the helicopter roll motion, at about 1 Hz; and the retreating lag motion. The resulting
orbit followed by the system after the bifurcation belongs to a multi-dimensional two-torus T 2 :
Since the two frequencies are not simply commensurate, the resulting regime orbit is quasi-
periodic, at least from a numerical point of view. Fig. 12 shows two phase space projections of the
system orbit before and after the bifurcation. In the first picture it is clear how the system is going
toward a fixed point (left bottom corner), as it is expected for a Poincare! map of a periodic orbit.
The second instead shows how the Poincare! map seems to form a closed curve made of discrete
points, a typical behaviour of a quasi-periodic solution, where the map represents a section of the
T 2 torus.
After this first bifurcation the system can still be considered globally stable, since the amplitudes
of the resulting oscillations are small. A second bifurcation will lead the system to a major
instability that cannot be recovered, as shown in Ref. [33]. However, the stability analysis of the
new quasi-periodic orbit is beyond the scope of this work so it is not presented here.
 ARTICLE IN PRESS

1036 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

2.8935 2.915

2.91
2.893
2.905
lead-lag angle [deg]

lead-lag angle [deg]

2.8925 2.9

2.895
2.892
2.89

2.8915 2.885

2.88
2.891
2.875

2.8905 2.87
0.3544 0.3548 0.3552 0.3556 0.33 0.34 0.35 0.36 0.37 0.38
(a) flap angle [deg] (b) flap angle [deg]

Fig. 12. Phase space projection of the system transient at two different values of lead-lag damping coefficient d: .;
Poincar"e sample; (a) d ¼ 8e3 ½N m=rad=s periodic behaviour; (b) d ¼ 6:5e3 ½N m=rad=s quasi-periodic behaviour.

7. Concluding remarks

The multibody approach allows the use of modular elements to build a complete system with all
the details required for each subpart. In this way it is possible to obtain a model that can be used
for different analyses at each development stage, avoiding risky physical oversimplifications.
When periodic systems are analyzed, the assessment of their stability boundaries for different sets
of operating parameters is often needed before going through any further investigation on the
system behaviour. Unfortunately, dealing with such large models, a very high computational
burden is required to gain this information by means of classical techniques based on the
reconstruction of the entire monodromy matrix. By means of the multibody/multidisciplinary
modelling paradigm the designer may deal with a single model useful in a wide range of operating
conditions. However, to make this tool really efficient, it must necessarily be coupled to a
methodology capable of synthesizing the significant pieces of information in a small set of results.
This paper presented an effective method to overcome these hurdles on the base of a proper
orthogonal decomposition. It is shown how, by means of this technique, it is possible to select the
dominant dynamics which govern the system response in the form of a small set of combined
states/signals. By using classical identification methods, it is straightforward to obtain the
dominant Floquet multipliers from these signals. By means of POD it is possible to achieve two
additional goals: first of all it gives information about the real dimensions of the subspace on
which the dynamic evolution for a particular set of parameters takes place; second, it allows one
to rule out all the spurious eigenvalues that may eventually arise when a non-minimal set system
of equations is analyzed. The application of the proposed method to the test case of the simply
supported beam subjected to an oscillating axial load gave results that have been matched
successfully with those that can be obtained for the same problem with more classical
methodologies. Few results obtained to assess the ground resonance stability of a large helicopter
model show the capabilities of the proposed method in dealing with large multidisciplinary
models.
 ARTICLE IN PRESS

G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038 1037

Further tests are required to verify the possibility of adopting a single set of POMs in the
identification of the Floquet multipliers for different operating conditions. This set does not need
to be the result of a single analysis made at a specific operating condition but can be obtained by
merging POMs computed with different parameter values. In fact, if the modal forms are already
determined, it will be possible to identify the Floquet multipliers while the simulation is running,
stopping it only when a sufficiently accurate result is achieved.

References

[1] G. Sansone, R. Conti, Non-linear Differential Equations, Pergamon Press, Oxford, UK, 1964 (translated from the
Italian by A. H. Diamond).
[2] A. Lyapunov, The General Problem of the Stability of Motion, Princeton University Press, Princeton, NJ, 1947.
[3] T.A. Burton, Stability and Periodic Solutions of Ordinary and Functional Differential Equations, Academic Press,
Orlando, FL, 1985.
[4] R. Seydel, Practical Bifurcation and Stability Analysis. From Equilibrium to Chaos, Springer, New York, 1994.
[5] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Field, 3rd
Edition, Springer, New York, 1990.
[6] W. Schiehlen, Multibody Systems Handbook, Springer, Berlin, 1990.
[7] G. Quaranta, P. Masarati, P. Mantegazza, Multibody analysis of controlled aeroelastic systems on parallel
computers, Multibody Systems Dynamics 8 (1) (2002) 71–102.
[8] K. Murphy, P. Bayly, L. Virgin, J. Gottwald, Measuring the stability of periodic attractors using perturbation-
induced transients: applications to two nonlinear oscillators, Journal of Sound and Vibration 172 (1) (1994) 85–102.
[9] H. Matthies, C. Nath, Dynamic stability of periodic solution of large scale nonlinear systems, Computational
Methods in Applied Mechanics and Engineering 48 (1985) 191–202.
[10] S. Subramanian, G.H. Gaonkar, Parallel fast-floquet analysis of trim and stability for large helicopter models, in:
22nd European Rotorcraft Forum, Brighton, United Kingdom, September 17–19, 1996, pp. 94.1–94.14.
[11] O.A. Bauchau, Y.G. Nikishkov, An implicit transition matrix approach to stability analysis of flexible multi-body
systems, Multibody Systems Dynamics 5 (3) (2001) 279–301.
[12] O.A. Bauchau, Y.G. Nikishkov, An implicit floquet analysis for rotorcraft stability evaluation, Journal of the
American Helicopter Society 46 (2001) 200–209.
[13] P. Holmes, J. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry,
Cambridge University Press, Cambridge, 1996.
[14] B. Epureanu, E. Dowell, K. Hall, Reduced-order models of unsteady viscous flows in turbomachinery using
viscous-inviscid coupling, Journal of Fluids and Structures 15 (2001) 255–273.
[15] B. Feeny, R. Kappagantu, On the physical interpretation of proper orthogonal modes in vibrations, Journal of
Sound and Vibration 211 (4) (1998) 607–616.
[16] M. Azeez, A. Vakakis, Proper orthogonal decomposition (POD) of a class of vibroimpact oscillations, Journal of
Sound and Vibration 240 (5) (2001) 859–889.
[17] R. Ruotolo, C. Surace, Using SVD to detect damage in structures with different operational conditions, Journal of
Sound and Vibration 226 (3) (1999) 425–439.
[18] V.V. Bolotin, Nonconservative Problems of the Theory of Elastic Stability, Pergamon Press, Oxford, 1963.
[19] R.P. Coleman, A.M. Feingold, Theory of self-excited mechanical oscillations of helicopter rotors with hinged
blades, REPORT 1351, NACA, 1958.
[20] P. Friedmann, C. Hammond, T.-H. Woo, Efficient numerical treatment of periodic systems with application to
stability problems, International Journal for Numerical Methods in Engineering 11 (1977) 1117–1136.
[21] G.H. Golub, C.F. van Loan, Matrix Computation, 2nd Edition, The John Hopkins University Press, Baltimore,
1991.
[22] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, Boston, MA, 1996.
[23] L. Lijung, System Identification—Theory for the User, 2nd Edition, Prentice-Hall, Upper Saddle River, NJ, 1999.
 ARTICLE IN PRESS

1038 G. Quaranta et al. / Journal of Sound and Vibration 271 (2004) 1015–1038

[24] P. Van Overschee, B. De Moor, N4SID: subspace algorithms for the identification of combined deterministic–
stochastic systems, Automatica 30 (1) (1994) 75–93.
[25] K.E. Brenan, S.L.V. Campbell, L.R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic
Equations, North-Holland, New York, 1989.
[26] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential Algebraic Problems, 2nd
revised edition, Springer, Berlin, 1996.
[27] C.L. Bottasso, M. Borri, L. Trainelli, Integration of elastic multibody systems by invariant conserving dissipating
algorithms. II. Numerical schemes and applications, Computer Methods in Applied Mechanics and Engineering 190
(2001) 3701–3733.
[28] P. Masarati, M. Lanz, P. Mantegazza, Multistep integration of ordinary, stiff and differential-algebraic problems
for multibody dynamics applications, in: XVI Congresso Nazionale AIDAA, Palermo, 24–28 Settembre, 2001.
[29] P. Hagedorn, Non-Linear Oscillations, Clarendon Press, Oxford, UK, 1981.
[30] G.L. Ghiringhelli, P. Masarati, P. Mantegazza, A multi-body implementation of finite volume beams, American
Institute of Aeronautics and Astronautics Journal 38 (2000) 131–138.
[31] J.P. Thomas, E.H. Dowell, K.C. Hall, Three-dimensional transonic aeroelasticity using proper orthogonal
decomposition based reduced order models, in: Proceedings of the 42nd AIAA/ASME/ASCE/AHS/ASC Structures,
Structural Dynamics, and Materials Conference and Exhibit, Seattle, WA, 16–19 April 2001.
[32] R.L. Bielawa, Rotary Wing Structural Dynamics and Aeroelasticity, AIAA, Washington, DC, 1992.
[33] S. Gualdi, P. Masarati, M. Morandini, G.L. Ghiringhelli, A multibody approach to the analysis of helicopter-
terrain interaction, in: 28th European Rotorcraft Forum, Bristol, UK, 2002.