Model reduction of nonlinear systems with external periodic excitations via construction of invariant manifolds

Journal of Sound and Vibration 330 (2011) 2596–2607

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Model reduction of nonlinear systems with external periodic excitations via construction of invariant manifolds Amit P. Gabale, S.C. Sinha n Nonlinear Systems Research Laboratory, Department of Mechanical Engineering, Auburn University, Auburn, AL 36849, USA

a r t i c l e in f o

abstract

Article history: Received 28 April 2010 Received in revised form 8 December 2010 Accepted 10 December 2010 Handling Editor: L.N. Virgin Available online 17 January 2011

A methodology for determining reduced order models of periodically excited nonlinear systems with constant as well as periodic coefficients is presented. The approach is based on the construction of an invariant manifold such that the projected dynamics is governed by a fewer number of ordinary differential equations. Due to the existence of external and parametric periodic excitations, however, the geometry of the manifold varies with time. As a result, the manifold is constructed in terms of temporal and dominant state variables. The governing partial differential equation (PDE) for the manifold is nonlinear and contains time-varying coefficients. An approximate technique to find solution of this PDE using a multivariable Taylor–Fourier series is suggested. It is shown that, in certain cases, it is possible to obtain various reducibility conditions in a closed form. The case of time-periodic systems is handled through the use of Lyapunov–Floquet (L–F) transformation. Application of the L–F transformation produces a dynamically equivalent system in which the linear part of the system is time-invariant; however, the nonlinear terms get multiplied by a truncated Fourier series containing multiple parametric excitation frequencies. This warrants some structural changes in the proposed manifold, but the solution procedure remains the same. Two examples; namely, a 2-dof mass–spring–damper system and an inverted pendulum with periodic loads, are used to illustrate applications of the technique for systems with constant and periodic coefficients, respectively. Results show that the dynamics of these 2-dof systems can be accurately approximated by equivalent 1-dof systems using the proposed methodology. & 2010 Elsevier Ltd. All rights reserved.

1. Introduction Dynamic analysis of many engineering systems often requires solution of a set of large number of nonlinear ordinary differential equations subjected to external excitations. Such large dimensional nonlinear systems are costly to solve in terms of computer time and storage. Further, dynamic characterization of such systems also poses serious challenges due to complex interrelationships that may exist between system parameters and excitation variables; which is true even for moderate sized systems. As a result, constant efforts are being made by various researchers to develop new techniques for reduced dimensional modeling of nonlinear dynamic systems. These reduced order models are comparatively easy to simulate and analyze as they consists of relatively small number of states. In this article we present one such technique for

n

Corresponding author. Tel.: + 1 334 844 3325; fax: + 1 334 844 3307. E-mail address: [email protected] (S.C. Sinha).

0022-460X/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.jsv.2010.12.013

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

2597

a class of problems which arises from systems subjected to external periodic excitation. The approach is based on a nonlinear projection of the full state-space onto a low dimensional state-space through the construction of an invariant manifold. The discussion encompasses systems with constant as well as the periodic coefficients. In general, nonlinear projection based order reduction techniques can be broadly classified into three categories, namely, center manifold, nonlinear normal modes (NNMs) and singular perturbation. Extensive literature is available on the theory and applications of these techniques. Here we mention only noteworthy references and some salient features for brevity. First technique, the ‘center manifold’, is a useful order reduction tool when the linear system matrix has a special spectral decomposition; wherein, some of the eigenvalues of the system lie on the imaginary axis (critical modes) while the remaining are on left half of the complex plane (stable modes). In such systems, the evolution of dynamics is governed by the critical modes alone and damped modes follow this ‘dominant dynamics’ passively. As a result, it is possible to construct a minimal sized model which depends only on the critical modes of the system. A detailed discussion of the center manifold technique can be found in Carr [1] and Guckenheimer and Holmes [2]. The second technique, NNMs, was originally introduced by Rosenberg [3]. He defined NNMs for autonomous, conservative nonlinear systems as a vibration in unison with fixed relations between generalized coordinates. Later, Shaw and Pierre [4] generalized the concept of NNMs to non-conservative systems. They used invariant manifold approach in which NNMs are assumed to be invariant subspaces, which are non-planer and tangential to corresponding normal modes of the linearized system at the equilibrium. This nonlinear extension of the modal analysis allows one to generate reduced order models by projecting non-essential modes on the invariant manifold which is parameterized by a set of desired modes. NNMs have gained popularity among structural/mechanical engineers and many references could be found (e.g., [5,6]) where this technique has been used to reduce the order of real engineering problems. The third technique, singular perturbation, is more popular among the control engineering community where it has been used, effectively, to design reduced order controllers for various nonlinear systems. Systems which come under the premise of singular perturbation admit (or can be transformed into) a special form in the state-space. In this form some of the derivatives of states are multiplied by a small positive parameter e. Such systems typically display a two time-scale phenomenon, by which, the fast subsystem approaches a steady state more quickly than the slow subsystem. Consequently the long term behavior of such systems can be approximated by the slow subsystem, alone, on a lower dimensional subspace. The singular perturbation provides a theoretical frame work to decompose an original large system into two subsystems; the slow subsystem representing the reduced order model. A detailed discussion on the subject can be found in books by Kokotovic´ et al. [7] and Khalil [8]. Among the techniques discussed above, the center manifold technique has been extended for systems with periodic coefficients [9]. An interesting account of the application of time-periodic center manifold theory for the bifurcation analysis of parametrically exited systems can be found in Pandiyan and Sinha [10] and Da´vid and Sinha [11]. Further, Sinha et al. [12] have also generalized the NNM technique, introduced by Shaw and Pierre [4], for systems with periodic coefficients. The treatment for singular perturbation technique also accommodates time-varying parameters in some of the literature. However, to authors’ best knowledge, the practical implementation of this technique to time-varying nonlinear systems has not been reported. Nonlinear dynamic systems subjected to external excitation, however, require special attention as none of the above techniques can directly be applied to this class of problems. Following the general idea of introducing a new state variable to represent the forcing term, such that the augmented system is homogenous, Shaw et al. [13] first applied the concept of NNMs to reduce the order of periodically excited autonomous systems. They proposed polynomial expansion to approximate a solution to a set of PDEs that governs the geometry of invariant manifold. Later, Jiang et al. [14] used a similar approach and employed Galerkin based numerical solution procedure to solve these PDEs. Their objective was to improve the validity of reduced order models in the large amplitude region. Further, Agnes and Inman [5] used the method of multiple scales to construct NNMs for such systems; however, they also employed augmented state-space to deal with the periodic forcing. The idea of introducing new state variables to transform the system into a homogenous form, however, is not very appealing for order reduction of systems which are subjected to multifrequency external and/or parametric excitations. Since every excitation frequency, present in the system, increases the dof of the equivalent homogenous system by one. This may make the construction of the invariant manifold a cumbersome process and may outweigh the advantages of order reduction all together. To overcome this drawback, recently, Redkar and Sinha [15] proposed an innovative approach where they directly constructed a time-varying invariant manifold in system’s modal coordinates. They assumed the invariant manifold as time modulated nonlinear functions of the dominant states which have a certain structure in terms of spatial and temporal terms. Their approach is quite general, in the sense that it can also be used for systems with periodic coefficients after the application of L–F transformation [16]. However, they observed that the formulation does not yield promising results when the unmodeled dynamics of the systems is near secondary resonance (sub- or super-harmonic) with external excitation frequencies. The work presented here draws on the idea of Redkar and Sinha [15,16] but significantly differs in the way we construct the time-varying manifolds. It is expected that reduced order models obtained with the proposed methodology will accurately estimate the response of the original system even if the unmodeled dynamics of the system is near primary or secondary resonance with external excitation. The remaining paper is organized as follows. Sections 2 and 3 outline the order reduction methodology for nonlinear systems with constant and periodic coefficients, respectively. Examples are presented in Section 4 which demonstrate applicability of the technique for these two types of systems. Discussion and conclusions are presented in Section 5.

2598

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

2. Order reduction of nonlinear systems with constant coefficients In this case, we are interested in approximating the forced response of a general nonlinear system of the form x_ ¼ Ax þ ef r ðxÞ þ FðtÞ

(1)

by a system of differential equations with smaller dimensions. In the above equation e is a small positive number, A is an n  n time-invariant matrix, f r ðxÞ is an n dimensional vector of nonlinear terms containing monomials in x up to order r and Fðt Þ is a T periodic force vector of commensurate dimensions. The proposed technique follows a two step procedure where we first transform the system into a state-space realization in which the new state variables can be grouped according to some measure of importance. In this work the important states variables are referred as the dominant states and the non-important as the non-dominant states. Second step involves finding an expression (invariant manifold) that relates the non-dominant states to the dominant states such that the resulting system of differential equations has lower order. Following this idea, we apply a similarity transformation x ¼ Mz to Eq. (1) which produces z_ ¼ Jz þ ef r ðzÞ þ FðtÞ 1

where J is in Jordan canonical form, f r ðdÞ ¼ M partitioned as ( ) " Jp z_ p ¼ 0 z_ q

1

f r ðdÞ and FðtÞ ¼ M 0 Jq

#(

zp zq

)

( þe

(2)

FðtÞ. Without any loss of generality, Eq. (2) can be

f p r ðzp ,zq Þ f q r ðzp ,zq Þ

(

) þ

Fp ðtÞ Fq ðtÞ

) (3)

where Jp is a p  pðp5nÞ Jordan block associated with the dominant states (zp ). These are the states we want to keep in the reduced order model. Jq is q  qðq ¼ npÞ Jordan block associated with the non-dominant states (zq ). These are the states we want to eliminate from the reduced order model. The proposed technique does not put any restrictions on the selection of the dominant and non-dominate states provided that they satisfy certain reducibility conditions. As a result, one can choose any set of states, which are of practical interest, as the dominant states. For systems which are subjected to periodic excitation, the response is usually dominated by the evolution of states which are in or near resonance (primary or secondary) with the external excitation or those which corresponds to eigenvalues of J with large negative real parts. Therefore, these states usually become natural choice of the dominate states. Once the dominant and non-dominant dynamics of the system is identified in Eq. (3), the non-dominant states can be eliminated from the dominant dynamics by expressing them in terms of the dominant states. We propose a nonlinear, time-varying constraint relationship between the dominant and non-dominant states in the following form: zq ¼ Hðzp ,tÞ ¼ h0 ðtÞ þ eðhr1 ðzp ,tÞ þ hr2 ðzp ,tÞ þ    þ hrr1 ðzp ,tÞ þ hrr ðzp ÞÞ

(4)

where the unknown h0 ðtÞ is assumed to be a function of time alone, hrs ðzp ,tÞ constitutes of monomials in zp of degree sðs ¼ 1,. . .,r1Þ with unknown time-periodic coefficients, while hrr ðzp Þ are monomials in zp of degree r with unknown constant coefficients. The constraint relationship represents a p dimensional manifold. For this manifold to be invariant (i.e., given zq ð0Þ ¼ Hðzp ð0Þ,0Þ, the solution ðzp ðtÞ,zq ðtÞÞ lies on the manifold 8 t Z 0), it must satisfy Eq. (3). Substituting Eq. (4) into Eq. (3) yields z_ p ¼ Jp zp þ ef pr ðzp ,Hðzp ,tÞÞ þ Fp ðtÞ

(5a)

@Hðzp ,tÞ @Hðzp ,tÞ þ z_ p ¼ Jq zq þ ef qr ðzp ,Hðzp ,tÞÞ þ Fq ðtÞ @t @zp

(5b)

Dynamics of the system on the invariant manifold is described by pðp o nÞ dimensional differential system (5a), which is a desired reduced order model of the original system (1), while function Hðzp ,tÞ must satisfy the PDE (5b). Thus, the solution of Eq. (5b) is essential for a successful order reduction. Eq. (5b) is a non-homogenous, nonlinear PDE with time-varying coefficients. To authors’ best knowledge this PDE is not tractable analytically and one has to use numerical techniques or series expansion to solve this equation approximately. In this work we propose the series expansion for this purpose. Since the invariant manifold is assumed in a specific temporal and spatial form, a combined Fourier and multivariable Taylor–Fourier series is used. This technique not only be easily automated for a nonlinearity of any degree, it also yields various reducibility conditions in a closed form in certain cases. However, it should be noted that such approximation is only valid for systems with weak nonlinearities and small amplitudes of motion. Accordingly we expand the nonlinear terms in Eq. (5b) and collect terms of like powers in zp which yields the following set of equations: z0p : z1p :

@hr1 ðzp ,h0 ðtÞÞ @h0 ðtÞ Jq h0 ðtÞef qr0 ðh0 ðtÞÞ þ e Fp ðtÞ ¼ Fq ðtÞ @zp @t

@hr1 ðzp ,tÞ @hr1 ðzp ,tÞ @hr2 ðzp ,tÞ þ Jp zp Jq hr1 ðzp ,tÞf qr1 ðzp ,h0 ðtÞÞ þ Fp ðtÞ ¼ 0 @t @zp @zp ^

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

zr1 : p

@hrr1 ðzp ,tÞ @hrr1 ðzp ,tÞ @hrr ðzp ,tÞ þ Jp zp Jq hrr1 ðzp ,tÞf qrr1 ðzp ,h0 ðtÞÞ þ Fp ðtÞ ¼ 0 @t @zp @zp @hrr ðzp Þ Jp zp Jq hrr ðzp Þ ¼ f qrr ðzp Þ zrp : @zp

2599

(6)

where terms up to the order of e are retained in the expansion. f qrs represents sth degree terms in zp arising from the nondominant dynamics in Eq. (3). The elements of invariant manifold are expanded as h0 ðtÞ ¼

q k X X

hj,l eilot ej

j ¼ 1 l ¼ k

hrs ðzp ,tÞ ¼

q X X k X

 m hr,j,ms ,l eilot zp  s ej

j ¼ 1 ms l ¼ k q X X

 m hr,j,mr zp  r ej

hrr ðzp Þ ¼

(7)

j ¼ 1 mr

pffiffiffiffiffiffiffi  m P mp 1 m2 1,o ¼ 2p=T and ej is the where s ¼ 1,2,. . .,r1, ma ¼ ðm1 ,m2 ,. . .,mp Þ, pi¼ 1 mi ¼ aða ¼ 1,2,. . .,r Þ, zp  a ¼ zm 1 z2    zp , i ¼ jth member of the natural basis. Substituting Eq. (7) into equation set (6) and comparing coefficients of similar terms of the expansion one gets a set of coupled nonlinear algebraic equations. These equations then can be solved numerically to obtain the Fourier and Taylor–Fourier coefficients of (7). However, depending upon the dimension of the original system, type of external excitation (single frequency/ multifrequency) and the degree of the nonlinearity, finding a solution to this set of coupled; nonlinear, algebraic equations may become a difficult task. Further, these algebraic equations typically have multiple solutions which give rise to multiple invariant manifolds. As a result, an investigator faces a challenge of choosing one manifold, out of several possible manifolds that yields the most accurate reduced order model of the system. Moreover, this technique does not reveal the reducibility conditions, which, if not satisfied can lead to singularities in solution of the nonlinear algebraic equations. In the case, however, when the external periodic excitation projected on the dominant dynamics is small ðFp ðtÞ  e Fp ðtÞÞ, the above drawbacks can be overcome as described next. Assuming Fp ðtÞ  e Fp ðtÞ, segregation of the terms in Eq. (5.b) according to their degree in zp yields z0p :

z1p :

zr1 p

:

@hr1 ðzp ,tÞ @t

@h0 ðtÞ Jq h0 ðtÞef qr0 ðh0 ðtÞÞ ¼ Fq ðtÞ @t þ

@hr1 ðzp ,tÞ Jp zp Jq hr1 ðzp ,tÞ ¼ f qr1 ðzp ,h0 ðtÞÞ @zp

^ @hrr1 ðzp ,tÞ @hrr1 ðzp ,tÞ þ Jp zp Jq hrr1 ðzp ,tÞ ¼ f qrr1 ðzp ,h0 ðtÞÞ @t @zp zrp :

(8a)

@hrr ðzp Þ Jp zp Jq hrr ðzp Þ ¼ f qrr ðzp Þ @zp

(8b)

(8c)

where, again, terms up to order of e are retained in the expansion. A similar set of equations was obtained by Gabale and Sinha [17] in the analysis of periodically excited nonlinear systems via normal forms. It is found that this set of equations is easier to solve compared to set (6). In the above, Eq. (8a) involves only temporal terms ðh0 ðt ÞÞ and is decoupled from the rest of equations. Therefore, any conventional technique such as perturbation, averaging, normal form or harmonic balance can be used to obtain a solution to this non-homogenous nonlinear differential equation. It also permits an investigator to choose the form of desired solution, i.e., fundamental, sub- or super-harmonic to correctly express the interaction between non-dominant dynamics and the external excitation. In the present analysis, we use method of harmonic balance for this purpose and the discussion is limited to the periodic solutions only. Substitution of h0 ðt Þ in equation set (8b) yields a set of decoupled PDEs with periodic coefficients. To solve these equations approximately, we expand the known ðf qrs ðzp ,h0 ðtÞÞÞ and the unknown ðhrs ðzp ,tÞÞ terms in multivariable Taylor– Fourier series and compare coefficients of the expansion, term-by-term. Eq. (8c) is time independent and similar to homological equations associated with the standard Poincare´ Normal form theory. Solution of this equation can be approximated by a polynomial expansion of known and unknown terms and a comparison of coefficients of the similar terms. In either case we obtain a set of linear algebraic equations which are solvable only if the corresponding linear coefficient matrices are invertible. It is important to note that this solution technique presents only one solution of PDE (5b) for any chosen solution of (8a). It also provides a greater control over the way we want to describe the interaction between external excitation and non-dominant dynamics. If matrix A has n linearly independent eigenvectors, matrix J and the linear coefficient matrices corresponding to (8b) and (8c) can be obtained in the diagonal form. Consequently the Taylor–Fourier coefficients of hrs ðzp ,tÞ and polynomial

2600

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

coefficients of hrr ðzp Þ can be put in a general form as hrs,j,ms ,l ¼

fqrs,j,ms ,l ilo þ ms Ukdom lj

(9a)

fqrr,j,mr mr Ukdom lj

(9b)

hrr,j,mr ¼

where kdom ¼ ðl1 l2    lp ÞT is a vector containing eigenvalues of the Jordan matrix Jp , lj ðj ¼ 1,2,. . .,qÞ are eigenvalues of the  m  m   P P P Pk P fqrs,j,m ,l eilot zp  s ej and f qrr zp ¼ q fqrr,j,m zp  r ej . It can be Jordan matrix Jq , f qrs ðzp ,h0 ðtÞÞ ¼ q j¼1

ms

l ¼ k

s

j¼1

mr

r

easily noticed from Eqs. (9a) and (9b) that one can find all Taylor–Fourier and polynomial coefficients of the expansion (or the linear coefficient matrices are invertible) only if the following reducibility conditions are satisfied: ilo þ ms Ukdom lj a0

(10a)

mr Ukdom lj a0

(10b)

In the event when condition (10a) and/or (10b) are not satisfied, it is called a ‘resonant’ case and the corresponding nondominant states cannot be expressed in terms of the dominant states. Reducibility condition (10a), based on the value of l, represents either the case of ‘true internal resonance’ ðl ¼ 0Þ or ‘true combination resonance’ ðla0Þ (Deshmukh et al. [18]), while, the reducibility condition (10b) always represents the case of ‘true internal resonance’. In the event of a true internal resonance, eigenvalues corresponding to the states involved in the resonance satisfy an integral relationship governed by (10a) or (10b). However, in the event of true combination resonance, these eigenvalues satisfies an integral relationship governed by (10a) which includes one of the external excitation frequencies. In such cases the proposed dominant and non-dominant dynamics cannot be separated by the suggested constraint relationship. Near resonance cases are also not desirable as they yield excessively large coefficients in the series expansion. Nevertheless, for ‘resonance’ and ‘near resonance’ cases one can always construct a reduced order model by making all the resonating states a part of either the dominant or the non-dominant state vector, in its entirety. In doing so, the unsatisfied reducibility conditions which arise due to the particular choice of dominant and non-dominant states, can be avoided. Further, it should be noted that the resonances associated with the reducibility conditions are different than the system’s resonances with the external excitation (i.e., fundamental, sub- or super-harmonic, etc.). Order reduction via proposed methodology is still possible in the presence of these primary or secondary resonances as long as Eq. (8a) is solvable. Once the constraint relationship between dominant and non-dominant states is obtained, it can be substituted into Eq. (5a) to get the desired reduced order model. The reduced order model can then be solved by any conventional nonlinear solution technique such as normal form, perturbation, etc., or can be integrated numerically. Using the similarity transformation the state response in the reduced order domain can be transformed into the original domain. It can be observed that if the external excitation is absent and the eigenvalues of Jp are critical, then the formulation reduces to the well known center manifold reduction. 3. Order reduction of nonlinear systems with periodic coefficients Next, we consider an order reduction problem for nonlinear systems with periodic external and parametric excitations given by x_ ¼ AðtÞx þ ef r ðx,tÞ þ FðtÞ

(11)

In the above equation, e is a small positive parameter, n  n matrix AðtÞ is T periodic such that AðtÞ ¼ AðT þ tÞ. f r ðx,tÞ is an n dimensional vector of nonlinear terms containing monomials in x up to the order r with k1 T periodic coefficients and FðtÞ is a k2 T periodic force vector of appropriate dimensions. k1 and k2 are assumed to be integers. However, the technique is not limited by this assumption and can be extended to systems with non-commensurate excitation, also. Proceeding along similar lines as center manifold and Poincare´ normal form theories for periodic systems [9,10,19–21], we first transform this system into its equivalent dynamic form for which the real part of the system is time-invariant. We use results from the Lyapunov–Floquet theory for this purpose. According to the L–F theory the state transition matrix of linear time-periodic systems; Uðt Þ, can be factored as [22]

UðtÞ ¼ LðtÞeR t

(12)

where LðtÞ is T periodic and R is a constant matrix. LðtÞ and R, in general, are complex. LðtÞ is called the Lyapunov–Floquet (L–F) transformation matrix and x ¼ LðtÞz converts linear part of Eq. (11) into a completely time-invariant form y_ ¼ Ry. Also, a 2T periodic LðtÞ can be determined, in which case LðtÞ and R are real. Application of the transformation x ¼ Lðt Þy to (11) yields y_ ¼ Ry þ eL1 ðtÞf r ðLðtÞy,tÞ þ L1 ðtÞFðtÞ

(13)

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

2601

where R is a real valued n  n constant matrix. For this system a similarity transformation can be found such that the linear part of (13) takes the Jordan canonical form. The resulting equation can be partitioned similar to Eq. (3) and written as ( ) " #( ) ( ) ) ( Jp 0 zp Fp ðt Þ z_ p f p r ðzp ,zq ,tÞ e ¼ þ (14) þ 0 Jq zq Fq ðt Þ z_ q f q r ðzp ,zq ,tÞ where zp is a p dimensional vector of dominant states and zq is a q dimensional ðp þ q ¼ nÞ vector of non-dominant states. Selection of dominant states, in this case, essentially follows similar guidelines as described in the previous section. For this purpose, eigenvalues of matrix R, which are often referred as ‘Floquet exponents’, can be used. Accordingly one can select dominant states based on the relative magnitudes of the real parts of Floquet exponents, or their resonances with external excitation. However, unlike systems with constant coefficients, in this case, the nonlinearities contain periodic coefficients. Therefore, in addition to the primary and secondary resonances one has to check for the possibility of combination resonances, too [18]. Moreover, the external excitation in Eq. (14) is inherently multifrequency due to the periodic nature of Lðt Þ. Eq. (14) differs from (3) in nonlinear terms f p r and f q r , which now contains time periodic coefficients. As a result, we propose a constraint relationship between dominant and non-dominant states of the form zq ¼ Hðzp ,tÞ ¼ h0 ðtÞ þ eðhr1 ðzp ,tÞ þ hr2 ðzp ,tÞ þ    þ hrr1 ðzp ,tÞ þ hrr ðzp ,tÞÞ

(15)

in which all the terms are similar to the one represented by (4) except the rth degree terms in zp , which are now assumed to have time-periodic coefficients. This constraint relationship is a p dimensional invariant manifold and must also satisfy the PDE (5b). We follow the same techniques as described before to solve this PDE, approximately. However, in this case, segregation of Eq. (5b), in terms of power of zp , yields an equation with time-periodic coefficients for the rth degree terms.   Further, when external excitation terms projected on the dominant dynamics are small Fp ðt Þ  e Fp ðt Þ one can approximate PDE (5b) by equations similar to (8a)–(8c); except that, Eq. (8c) now contains time-periodic terms. Consequently the reducibility conditions (10a) and (10b) can be represented by a single generalized reducibility condition as ilo þ ms Ukdom lj a0

(16)

Pp

where ms ¼ ðm1 ,m2 ,. . .,mp Þ, i ¼ 1 mi ¼ s ðs ¼ 1,2,. . .,rÞ. If condition (16) is not satisfied, it corresponds to the ‘resonant’ case and states involved cannot be decoupled. Nonetheless, in certain cases, one can avoid the resonance by proper choice of the dominant and non-dominant state vectors as pointed out earlier. It should be noted that when k1 and k2 are not integers, the external excitation projected on the dominant and nondominant dynamics ðFp ðtÞ,Fq ðtÞÞ as well as coefficients of nonlinearities f p r and f q r become quasi-periodic. As a result, one has to solve Eq. (8a)) as a nonlinear ordinary differential equation which is subjected to non-commensurate multifrequency excitation. It is also seen that for the solution of Eqs. (6), (8b) and (8c) one has to expand the known and unknown terms of the equations in multivariable, multifrequency Taylor–Fourier series. For small dominant excitation Fp ðt Þ, the reducibility condition can be derived which is multifrequency generalization of the condition derived before. This is not pursued any further here. At this point authors would also like to emphasize the important differences between the current methodology and the one proposed in references [15,16]. First, the spatial and temporal structure of the invariant manifold proposed in this work is different compared to the one suggested in [15,16]. Further, in the present work the invariant manifold is approximated by solving either the set of PDEs (6) or the nonlinear ordinary differential equation (ODE) (8a) and the set of PDEs (8b) and (8c). As discussed in Section 2, this allows one to correctly express interactions between the non-dominant dynamics of the system and the external excitation. The methodology proposed in [15] and [16], however, uses a set of linear ODEs and PDEs to approximate the invariant manifold. Although this methodology is simpler compared to the one proposed here, reduced order models thus obtained yields good results only when the external excitation is in no resonance or in primary resonance with the non-dominant dynamics of the system. 4. Applications To demonstrate possible applications, we first determine a reduced order model for a 2-dof spring–mass–damper system with nonlinear springs, as shown in Fig. 1. The equations of motion (with constant coefficients) can be written in the state-space form as 9 8 9 38 9 8 8_ 9 2 0 0 0 1 0 0 x x1 > > > > > > > > > > > > > > 6 > > > > 1> > > > > > > > > > 7 = < = < x_ 2 = 6 0 0 0 0 1 0 7< x2 = < 7 þ ¼6 þ e 3 3 6 7 ð x x Þ b =m x b=m _  k ð þ k Þ=m k=m c =m 0 f cos ð O t Þ > > > > > > > x x 1 2 1 1 1 1 > > 1 1 1 1 1 1 1 > > > 5> > > > > 4 > > > > 3> > 3> ; ; > : : ; > ; : b=m ðx x Þ3 b =m x3 > : f2 cosðO2 t Þ > ðk2 þ kÞ=m2 0 c2 =m2 k=m2 x4 x_ 4 2 2 1 2 1 2 (17) The nonlinearities are assumed to be of the cubic form. A similar system was considered by Jiang et al. [14] for order reduction via construction of NNM. It involved augmentation of the system to its higher order homogenous form before the order reduction procedure can be applied. Moreover, the analysis was limited to a single frequency excitation ðf2 ¼ 0Þ

2602

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

f1cos (Ω1t)

f2cos (Ω2t) b2

b

b1

k2

k1 c1

m1

k

m2

c2

Fig. 1. A 2-dof spring–mass–damper system.

and they considered only the case of primary resonance near the linear modes. In the present example, we directly construct a two dimensional time-varying invariant manifold using the proposed technique and the 2-dof system is approximated by a single-dof system. Cases of multifrequency excitation are also considered which include secondary resonances near one of the linear modes. Application of modal transformation ðx ¼ MzÞ to Eq. (17) yields an equation similar to (3). For the present case, n ¼ 4, r ¼ 3 and Fp ðt Þ and Fq ðt Þ are given by ( 1 ) ( 1 ) 1 1 f2 cosðO2 t Þ f2 cosðO2 t Þ M13 f1 cosðO1 t Þ þ M14 M33 f1 cosðo1 t Þ þM34 Fp ðt Þ ¼ ð t Þ ¼ and F q 1 1 1 1 f1 cosðO1 t Þ þ M24 f2 cosðO2 t Þ f1 cosðo1 t Þ þM44 f2 cosðO2 t Þ M23 M43 1 represents elements M1 . Mab Since the system contains cubic nonlinearity we propose an invariant manifold of the form (cf., Eq. (4))

zq ¼ Hðzp ,tÞ ¼ h0 ðtÞ þ eðh31 ðzp ,tÞ þ h32 ðzp ,tÞ þh33 ðzp ÞÞ

(18)

As discussed before, Eq. (18) must satisfy equation set (6) or (8). If all elements of (18) can be determined, then a reduced order model of the system can be obtained by replacing the non-dominant states from the dominant dynamics (5a) as functions of dominant sates, alone. In the following, some case studies for various system parameter values are presented. Parameters are chosen such that the system satisfies all the reducibility conditions and the order reduction using proposed constraint relationship is possible. In contrast to guidelines proposed earlier for the selection of dominant/ non-dominant states, states which are near resonance (primary or secondary) with the external excitation are assumed as the non-dominant states. This is intentionally done in order to highlight effectiveness of the technique in expressing unmodeled dynamics of the system via reduced order models. First we consider a case of primary resonance. For system parameters m1 ¼ m2 ¼ 1, k ¼ 1, k1 ¼ 4, k2 ¼ 13, c1 ¼ 0:2, c2 ¼ 0:2, b ¼ 1, b1 ¼ 1, b2 ¼ 5, f1 ¼ 0:5, f2 ¼ 1:1, the undamped natural frequencies of the system are o1 ¼ 2:2114 rad=s, o2 ¼ 3:7563rad=s and matrix J in Eq. (2) can be written as diag½0:1 þ 2:2091i,0:12:2091i,0:1 þ 3:7550i,   0:1þ 3:7550i. We chose the frequency of external excitation to be 4:25 rad=s O1 ¼ O2 ¼ 4:25 rad=s . The states corresponding to eigenvalues 0:1 72:2091i are taken as the dominant states while the states corresponding to eigenvalues 0:17 3:7550i are assumed as the non-dominant states ðp ¼ q ¼ 2Þ and the original 2-dof system is approximated by a 1-dof system. Invariant manifolds are calculated by solving, both, equation set (6) as well as (8) using proposed techniques. Since frequency of excitation ð4:25 rad=sÞ is away from the 3o2 and o2 =3 terms, we seek a fundamental solution to Eq. (8a), which is found by using the method of harmonic balance. Of course, one can use any other analytical method such as perturbation, averaging or normal form to find solution of (8a). The reduced order model is integrated numerically and using the modal transformation all states in x domain are recovered. The original system is also integrated numerically for the comparison. The initial conditions are taken such that the system’s trajectories are initiated on the invariant manifold (i.e., xð0Þ ¼ Mfzp ð0Þ,Hðzp ð0Þ,0ÞgT ). Fig. 2 shows time histories for states x1 and x2 as obtained from reduced and original models which show very good agreement. In order to gain insight into the frequency content of the response, fast Fourier transforms (FFTs) of the time histories of the original and reduced order systems are also shown (Fig. 3). In this case, the total response is dominated by the fundamental solution of the system, alone. Therefore, only one dominant peak is observed (frequency 4:25 rad=s) in the FFTs. Variation in amplitude of stable, steadystate response of the system near main resonance frequency is also shown in Fig. 4, which shows close agreement between the reduced and full-order model response. Thus both reduced order models, as derived by solving equation sets (6) and (8), are found equally successful in approximating response of the original model in this case. For parameter values m1 ¼ m2 ¼ 1, k ¼ 0, k1 ¼ 4, k2 ¼ 10, c1 ¼ 0:2, c2 ¼ 0:2, b ¼ 4, b1 ¼ 2, b2 ¼ 3, f1 ¼ 0:5, f2 ¼ 3:4, the undamped natural frequencies of the system are o1 ¼ 2 rad=s, o2 ¼ 3:1623rad=s and matrix J can be written in the form diag½0:1þ 1:9975i,0:11:9975i,0:1þ 3:1607i,0:13:1607i. We choose the external excitation frequencies as O 1 ¼ 3:0 rad=s and O 2 ¼ 1:2 rad=s ð3O 2  o2 Þ such that the system response contains a dominant super-harmonic component of order 3. States corresponding to eigenvalues 0:1þ 1:9975i are taken as the dominant states and to those corresponding to 0:1 þ 3:1607i as the non-dominant states. The non-dominant dynamics of the system is near

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

Fig. 2. System response ðO1 ¼ O2  o2 Þ. —, full order model; + , reduced order model;

J,

Fig. 3. Fast Fourier transforms ðO1 ¼ O2  o2 Þ. —, full order model; + , reduced order model;

2603

reduced order model assuming Fp to be small.

J,

reduced order model assuming Fp to be small.

Fig. 4. Amplitude variation vs. excitation frequency. —, full order model; + , reduced order model;

J,

reduced order model assuming Fp .

super-harmonic resonance with the external excitation. Therefore, we seek a third order super-harmonic solution for Eq. (8a), which is obtained using the method of harmonic balance. It should be noted that the use of Eq. (6) and the following numerical solution technique, does not offer such control on the construction of the invariant manifold. Figs. 5 and 6 show

2604

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

Fig. 5. System response ðO1 aO2 ,3O2  o2 Þ. —, full order model; + , reduced order model;

J,

Fig. 6. Fast Fourier transforms ðO1 aO2 ,3O2  o2 Þ. —, full order model; + , reduced order model;

reduced order model assuming Fp to be small.

J,

reduced order model assuming Fp to be small.

time histories and FFTs, respectively, of the system response as obtained from simulation of the reduced and full order models. It is observed form Fig. 6 that the system response is dominated by three distinct harmonics at 1:2, 3 and 3:6 rad=s. The first two harmonics corresponds to the external excitation frequencies while the third harmonic is due to super-harmonic resonance of non-dominant states with the external excitation. It is evident that the both reduced order models capture these components of the system dynamics successfully. Fractional harmonics at 0:6 and 1:8 rad=s are also rendered well by these models. The reduced order model obtained by assuming small dominant forcing, however, overestimates harmonics at 2:4rad=s in x1 and at 4:8rad=s and 5:4 rad=s in x2 . Since amplitudes of these harmonics are much smaller compared to the dominant one it has negligible effect on the overall system response. Next, we choose parameter values m1 ¼ m2 ¼ 1, k ¼ 0, k1 ¼ 43, k2 ¼ 4, c1 ¼ 0:1, c2 ¼ 0:1, b ¼ 3, b1 ¼ 3, b2 ¼ 7:972, f1 ¼ 2, f2 ¼ 7:5 such that the undamped natural frequencies of the system are o1 ¼ 6:5574, o2 ¼ 2 and the matrix J can be written in the form diag½0:05 þ 6:5572i,0:056:5572i,0:05 þ1:9994i,0:051:9994i. The forcing frequencies are chosen as O1 ¼ 4:4 rad=s and O2 ¼ 6:6 rad=s ðO 2  3o2 Þ, which, for the above parameter values, introduce sub-harmonic components of order 1/3 in the system response. The states corresponding to the eigenvalues 0:05 7 6:5572i are assumed as the dominant sates and those corresponding to 0:05 7 1:9994i are assumed as the non-dominant states and once again the system is approximated by a 1-dof system. Since the non-dominant undamped natural frequency o2 is close to O2 =3, we seek a sub-harmonic solution of order 1/3 of Eq. (8a). Fig. 7 shows the system response as calculated from the original and reduced order models for the current parameter values, while Fig. 8 shows FFTs of the steady-state response. In this case, too, both reduced order models estimates response of the original system quite well. Next we consider reduced order modeling of a dynamic system with periodic coefficients. Fig. 9 shows a 2-dof coupled inverted pendulum in the horizontal plane. The pendulum is subjected to periodic loading and torques as shown in the figure. The equations of motion for this system, when expanded up to cubic terms, about the fixed point fy1 , y2 , y_ 1 , y_ 2 g ¼ f0,0,0,0g can be shown to be [12] ! i P ðt Þ h k k h y3 2 3 1 c11 ðy1 y2 Þ þ c12 ðy1 y2 Þ þ c13 ðy1 y2 Þ  y€ 1 þ 12 y_ 1 þ t12 y1 þ y1  1 ¼ M1 cosðO1 tÞ 4m ml 6 ml ml ! 3 h i €y þ h2 y_ þ kt2 y þ k c ðy y Þ þc ðy y Þ2 þ c ðy y Þ3  P2 ðt Þ y  y2 ¼ M cosðO t Þ (19) 2 2 2 11 2 1 12 2 1 13 2 1 2 2 2 4m ml 6 ml2 ml2

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

Fig. 7. System response ðO1 aO2 , O2  3o2 Þ. —, full order model; + , reduced order model;

Fig. 8. FFTs ðO1 aO2 , O2  3o2 Þ. —, full order model; +, reduced order model;

J,

J,

2605

reduced order model assuming Fp to be small.

reduced order model assuming Fp to be small.

Fig. 9. A 2-dof coupled inverted pendulum.

where Pi ðtÞ ¼ Pi1 þPi2 cosðotÞ. kti and hi denote torsional stiffness and damping of the system, respectively ði ¼ 1,2Þ. Above set of equations contains time-periodic coefficients in linear as well as nonlinear terms. Typical, nonlinear decoupling of its states using traditional NNM technique requires transformation of this system into a higher order system which is homogenous as well as autonomous. With the proposed techniques, however, one can directly construct a 2-dimensional, time-varying invariant manifold and above 2-dof time-periodic system can be approximated by a 1-dof system as illustrated below. For a typical parameter set, Eq. (19) can be written in the state-space form as 8_ 9 2 0 x1 > > > > 6 >_ > < 0 x2 = 6 ¼6 6 ð1 þ 1cosð2pt ÞÞ _ x > > 3 > 4 > > > : ; 0 x_ 4

0

1

0

0

0

1

ð3 þ 2 cosð2pt ÞÞ

0

9 8 9 38 9 8 0 0 x1 > > > > > > > > > > > > > > > > > > 7> > > = > = < 0 0 1 7< x 2 = < 7 þ þ 3 3 7 0 5> x3 > > 0:5 cosðpt Þ > > > > > > cosð2pt Þx1 =62:5ðx1 x2 Þ > > > > > > > ; > > : > ; > ; : 2 cosð2pt Þx3 =6 þ2:5ðx x Þ3 > : 0:8 cosð0:5pt Þ > 1 x4 1 2 2 0

(20)

2606

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

T where fx1 ,x2 ,x3 ,x4 gT ¼ fy1 , y2 , y_ 1 , y_ 2 g , Pi1 ¼ 0 and ci2 ¼ 0. As described in Section 3, we first apply the L–F transformation ðx ¼ Lðt ÞyÞ to the above equation such that the linear part becomes time-invariant. Application of the modal transformation ðy ¼ MzÞ to this equation transforms it into the form represented by Eq. (21) as shown below:

8_ 9 2 z1 > 0:5 > > > > = 6 0:8739 < z_ 2 > 6 ¼6 > z_ 3 > 0 > 4 > > ; : > z_ 4 0

þ

8 > > > > <

0:8739

0

0:5

0

0

0:5

0

1:4621

38 9 z1 > > > > > > = < 7 0 7 z2 7 z > 5 1:4621 > 3 > > > ; : > z4 0:5 0

1:4111z32 1:2409z2 z22 1:0020z2 z24 þ 3:1826z22 z2 cosðpt Þ þ    0:0106z22 z4 cosðpt Þ0:0101z22 z4 cosð3pt Þ þ 0:0115z22 z3 sinðpt Þ þ   

9 > > > > =

0:9912z22 z4 þ 0:1496z23 z4 þ 0:1131z34 þ0:3610z32 cosðpt Þ0:8672z32 sinðpt Þ þ    > > > > > > > ; : 1:2276z2 z3 0:1892z3 0:1495z2 z2 þ 1:0494z3 cosðpt Þ þ0:2983z3 sinðpt Þ þ    > 2 3 4 2 2 9 8 0:8296cosðpt Þ þ0:0113cosð3pt Þ þ    > > > > > > > > = < 0:0031 sinðpt Þ þ 0:0031 sinð3pt Þ þ    þ 0:1615cosð0:5pt Þ þ 0:1615cosð1:5pt Þ0:4212sinð0:5pt Þ þ    > > > > > > > ; : 0:4696 cosð0:5pt Þ þ0:4696 cosð1:5pt Þ þ 0:1449sinð0:5pt Þ þ    >

(21)

We use the approach introduced by Sinha et al. [23] to calculate the real L–F transformation matrix Lðt Þ which is 2T periodic in this case. The states corresponding to exponents (eigenvalues of matrix R)0:5 7 0:8739 are assumed as the dominant while those corresponding to 0:57 1:4621 are assumed as the non-dominant. Due to the cubic nature of the nonlinearity, we propose a constraint relationship between the dominant and non-dominant states in a form similar to (18) with exception of the last term ðh33 ðzp ÞÞ which, in this case, is assumed to have time-periodic coefficients. As described in Section 3, this constraint relationship is governed by equations similar to equation set (6) or (8). In the present example, the invariant manifolds are constructed by solving both sets of governing equations and thus two separate reduced order models are obtained. Eq. (8a) is solved using the method of harmonic balance where the solution is expanded up to 10 harmonic terms. Simulated response of the original system and its reduced order models in shown in Fig. 10. Fig. 11 shows FFTs of the steady-state trajectories.

Fig. 10. System response (example with periodic coefficients). —, full order model; + , reduced order model; be small.

Fig. 11. FFTs (Example with periodic coefficients). —, full order model; + , reduced order model;

J,

J,

reduced order model assuming Fp to

reduced order model assuming Fp to be small.

A.P. Gabale, S.C. Sinha / Journal of Sound and Vibration 330 (2011) 2596–2607

2607

It is interesting to note that, in this case, a typical multivariable Taylor–Fourier expansion of the governing PDEs similar to (6), in which 10 Fourier terms are kept in the expansion, yields a set of 420 coupled nonlinear algebraic equations. An identical expansion of the terms, when used in equation set (8), results in 42 coupled nonlinear algebraic equations and 378 linear algebraic equations, where the nonlinear equations are decoupled from the linear equations. Thus, if Fp ðt Þ is assumed to be small, the effort required to construct the invariant manifold is considerably smaller. However, such simplification may come with a tradeoff (as seen in this example), where the reduced model calculated assuming small Fp ðt Þ overestimates the 2p frequency component in x2 . 5. Conclusions An order reduction technique for nonlinear systems subjected to external periodic excitations is proposed by constructing a time-dependant invariant manifold. This direct approach is in contrast to procedures where such systems are transformed into higher order homogenous forms by defining additional coordinates. The PDE governing the geometry of the manifold is nonlinear and contains time-periodic coefficients. An approximate solution of this PDE is suggested in terms of a multivariable Taylor–Fourier series. It is shown that in the case when the external excitation projected on the dominant dynamics is small, the equation governing the manifold can be reduced to a set of equations which are simpler to solve and yield reducibility conditions in a closed form. These reducibility conditions, which represent interactions among dominant, non-dominant states and the excitation terms, are indicators of the feasibility and extent of accuracy of such order reductions. The proposed technique does not put any constraint on the selection of dominant and nondominant states, provided that the chosen partitioning satisfies all reducibility conditions. Therefore, an investigator is free to choose any set of states of practical importance as dominant states. The resulting reduced order models can be numerically integrated or can be analyzed using traditional nonlinear analysis techniques. Systems with time-periodic coefficients are handled via an application of the L–F transformation that converts the linear parts of systems into timeinvariant forms. Due to the periodic nature of the L–F transformation matrix, however, the system nonlinearities in the transformed domain contain additional time-periodic coefficients which warrant a modification in the structure of the invariant manifold. The reducibility conditions are obtained in terms of Floquet exponents in this case. The effectiveness of the method is demonstrated by constructing reduced order models for some typical engineering systems. Results obtained from reduced order models are compared with numerical simulations of full order models. On the basis of these results, it is concluded that the proposed technique provides accurate reduced order models for nonlinear systems of fairly general form provided that the reducibility conditions are satisfied. In the case of systems with periodic coefficients, the technique is not limited by the magnitude of periodic terms which is an added advantage of this technique in terms of its versatility. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

J. Carr, Applications of Center Manifold Theory, Springer-Verlag, New York, 1981. J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983. R.M. Rosenberg, On nonlinear vibrations of systems with many degrees of freedom, Advances in Applied Mechanics 9 (1966) 155–242. S.W. Shaw, C. Pierre, Normal modes for non-linear vibratory systems, Journal of Sound and Vibration 164 (1) (1993) 85–124. G.S. Agnes, D.J. Inman, Performance of nonlinear vibration absorbers for multi-degrees-of-freedom systems using nonlinear normal modes, Nonlinear Dynamics 25 (2001) 275–292. E. Pesheck, C. Pierre, S.W. Shaw, Modal reduction of a nonlinear rotating beam through nonlinear normal modes, Journal of Vibration and Acoustics 124 (2002) 229–236. P. Kokotovic´, H.K. Khalil, J. O’Reilly, Singular Perturbation Methods in Control Analysis and Design, SIAM, Philadelphia, 1999. H.K. Khalil, Nonlinear Systems, Prentice-Hall, New Jersey, 2002. I.G. Malkin, Some basic theorems of the theory of stability of motion in critical cases, Stability and Dynamic Systems, Translations, American Mathematical Society, Series 1 5 (1962) 242–290. R. Pandiyan, S.C. Sinha, Analysis of time-periodic nonlinear dynamic systems undergoing bifurcation, Nonlinear Dynamics 8 (1995) 21–43. A. Da´vid, S.C. Sinha, Versal deformation and local bifurcation analysis of time-periodic nonlinear systems, Nonlinear Dynamics 21 (2000) 317–336. S.C. Sinha, S. Redkar, E.A. Butcher, Order reduction of nonlinear systems with time periodic coefficients using invariant manifolds, Journal of Sound and Vibration 284 (2005) 985–1002. S.W. Shaw, C. Pierre, E. Pesheck, Modal analysis-based reduced-order models for nonlinear structures—an invariant manifold approach, Shock and Vibration Digest 31 (1) (1999) 3–16. D. Jiang, C. Pierre, S.W. Shaw, Nonlinear normal modes for vibratory systems under harmonic excitation, Journal of Sound Vibration 288 (2005) 791–812. S. Redkar, S.C. Sinha, A direct approach to order reduction of nonlinear systems subjected to external periodic excitations, Journal of Computational and Nonlinear Dynamics 3 (2008) 031011-1-8. S. Redkar, S.C. Sinha, Order reduction of parametrically excited nonlinear systems subjected to external periodic excitations, Proceedings of the ASME 2009 International Design Engineering Technical Conferences, San Diego, USA, 2009. A.P. Gabale, S.C. Sinha, A direct analysis of nonlinear systems with external periodic excitations via normal forms, Nonlinear Dynamics 55 (2009) 79–93. V. Deshmukh, E.A. Butcher, S.C. Sinha, Order reduction of parametrically excited linear and nonlinear structural systems, Journal of Vibration and Acoustics 128 (2006) 458–468. D.K. Arrowsmith, C.M. Place, An Introduction to Dynamical Systems, Cambridge Press, Cambridge, 1990. V.I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, Springer-Verlag, New York, 1998. S.C. Sinha, R. Pandiyan, Analysis of quasilinear dynamical systems with periodic coefficients via Liapunov–Floquet transformation, International Journal of Non-Linear Mechanics 29 (5) (1994) 687–702. V.A. Yakubovich, V.M. Starzhinskii, Linear Differential Equations with Periodic Coefficients—Parts I and Part II, John Wiley and Sons, New York, 1975. S.C. Sinha, R. Pandiyan, J.S. Bibbs, Liapunov–Floquet transformation: computation and application to periodic systems, Journal of Vibration and Acoustics 118 (1996) 209–219.