Comparison of Galerkin, POD and Nonlinear-Normal-Modes Models for Nonlinear Vibrations of Circular Cylindrical Shells

,


INTRODUCTION
Reduced-order models are an attractive research topic in nonlinear dynamics of fluid and solid systems.From far, the two most popular methods used to build reduced-order models (ROMs) are the proper orthogonal decomposition (POD) and the nonlinear normal modes (NNMs) methods.The first one (POD, also referred to as the Karhunen-Loève method) uses a cloud of points in phase space, obtained from simulations or from experiments, in order to build the reduced subspace that will contain most information (Zahorian and Rothenberg 1981;Aubry et al. 1988;Sirovich 1987;Breuer and Sirovich 1991;Azeez and Vakakis 2001;Sarkar and Païdoussis 2003;Kerschen, et al. 2003;Amabili et al. 2003;Sarkar and Païdoussis 2004;Kerschen, et al. 2005;Amabili et al. 2006).The method is, in essence, linear, as it furnishes the best orthogonal basis, which decorrelates the signal components and maximizes variance.Amabili et al. (2003Amabili et al. ( , 2006) ) compared Galerkin and POD models of a water-filled circular cylindrical shell from moderate to extremely large vibration amplitudes.Accurate POD models can be build by using only POD modes with significant energy.In particular, Amabili et al. (2006) found that more proper orthogonal modes are necessary to reach energy convergence using time series extracted from more complex responses (chaotic or quasi-periodic) than from the periodic ones.Therefore by using complex response it is possible to build models with larger dimension, suitable to 1 Copyright © 2006 by ASME describe with accuracy large variations of the system parameters.
The second method (NNMs), allows to construct and define the researched subspaces from specific properties of the dynamical systems, by adapting the reduction theorems provided by the mathematics: center manifold theorem (Carr 1981;Guckenheimer and Holmes 1983) and normal form theory (Poincaré 1892;Iooss and Adelmeyer 1998;Elphick et al. 1987).Their application to vibratory systems led to two definitions of Nonlinear Normal Modes, which are equivalent in a conservative framework: either a family of periodic orbits in the vicinity of the equilibrium point (Rosenberg 1966;Mikhlin 1995;Vakakis et al. 1996), or an invariant manifold containing these periodic orbits (Shaw and Pierre, 1991).Numerous asymptotic methods have been proposed for their computation, by application of the center manifold theorem (Shaw and Pierre 1993), the normal form theory (Jézéquel and Lamarque, 1991, Touzé et al. 2004, Touzé and Amabili 2005), the conservation of energy for conservative systems (King and Vakakis 1994), or the method of multiple scales (Lacarbonara et al., 2003).Numerical procedures have also been proposed, recently by Jiang et al (2005a, b), who extended the method proposed by Pesheck et al. (2002) for conservative cases.Bellizzi and Bouc (2005) propose a numerical resolution of an extended KBM (Krylov-Bogoliubov-Mitropolsky) method, while Slater (1996) used continuation techniques to generate the NNM.
Application of the POD and NNMs methods to reducedorder modelling enabled to show that a few degrees of freedom are generally enough to catch the nonlinear behaviour of many structures, versus the several necessary in the corresponding Galerkin models.
The aim of the present study is to provide a full comparison of the results given by the two reduction methods on a realistic example: a water-filled circular cylindrical shell.The reference solution is obtained by the Galerkin method.Its convergence has been carefully checked (Pellicano et al. 2002), and experimental comparisons have been performed (Amabili 2003).The construction of the POD model has been exhaustively explained in Amabili et al. (2003Amabili et al. ( , 2006)), whereas the asymptotic NNMs procedure used here is fully explained in Touzé and Amabili (2005).The peculiarity of the NNMs formulation is that damping is taken into account via an improvement of the real normal form calculation presented in Touzé et al. (2004).Comparisons will be drawn on two different cases.First, the ability of the methods to recover frequency-response curves will be investigated, for moderate values of the amplitude of the external force.Then, bifurcation diagrams for varying amplitude of the forcing, leading to complex dynamics, will be discussed.

BASIC EQUATIONS FOR NONLINEAR VIBRATIONS OF FLUID-FILLED SHELLS
A cylindrical coordinate system (O; x, r, θ) is chosen, with the origin O placed at the centre of one end of the shell.The displacements of the middle surface of the shell are denoted by u, v and w, in the axial, circumferential and radial directions, respectively; w is taken positive inwards.By using Donnell's nonlinear shallow-shell theory, the equation of motion for finite-amplitude transverse deflection is given by Amabili and Païdoussis (2003) where is the flexural rigidity, E Young's modulus, ν the Poisson ratio, h the shell thickness, R the mean shell radius, ρ the mass density of the shell, c the coefficient of viscous damping, p the radial pressure applied to the surface of the shell by the contained fluid, and f is an external local excitation: where δ is the Dirac delta function, f is the magnitude of the localized (point) force, and θ and give the angular and axial coordinates of the point of application of the force, respectively.The viscous damping model introduced in equation ( 1) is unrealistic; it will later be replaced by modal damping coefficients in the equations of motion. x In equation (1) the overdot denotes a time derivative, and F is the in-plane Airy stress function, which is given by the following compatibility equation (Amabili and Païdoussis 2003): In equations ( 1) and (3), the biharmonic operator is defined as Donnell's nonlinear shallow-shell equations are accurate only for modes with n ≥ 4. Attention is focused on simply supported, circumferentially closed circular cylindrical shells of length L.Moreover, u, v and w must be continuous in θ.
Excitations with frequency close to the natural frequency of the lowest modes of the shell are considered; low-frequency modes are associated with predominantly radial motion and are identified by the pair (m, n), where m is the number of axial half-waves and n is the number of circumferential waves.
The contained fluid is assumed to be incompressible, inviscid and irrotational, so that potential theory can be used to describe fluid motion.Liquid-filled shells vibrating in the lowfrequency range satisfy the incompressibility hypothesis very well.Nonlinear effects in the dynamic pressure and in the boundary conditions at the fluid-structure interface are neglected.The shell prestress due to the fluid weight is also neglected.The fluid motion is described by the velocity potential Φ, which satisfies the Laplace equation; cavitation is assumed not to occur.Both ends of the fluid volume, 2 Copyright © 2006 by ASME corresponding to the shell edges, are assumed to be open, so that a zero pressure is imposed there; this physically corresponds to a long shell periodically supported (e.g. with ring stiffeners) or it approximates a shell closed by very thin circular plates.The dynamic pressure p exerted by the contained fluid on the shell is given by Amabili (2003) ( ) where ρ F is the mass density of the internal fluid, I n is the modified Bessel function of order n, and I its derivative with respect to the argument.The Galerkin method, employing any set of basis functions ϕ i , approximates the nonlinear partial differential equation (PDE) by transforming it into a finite set of coupled ordinary differential equations (ODEs), with the solution being expressed as 1 ( , ) ( ) ( ), where t is time, is the vector of spatial coordinate ( ξ , x θ ) describing the shell middle surface Ω, q i (t) are the generalized coordinates, and K is the number of generalized coordinates (degrees of freedom), i.e. the number of basis functions assumed.The linear modal base is the best choice for discretizing the shell, as these are the eigenfunctions of the linear operator of the PDE.The orthogonality property of the eigenmodes allows decoupling the ODEs at the linear stage.Other sets of basis functions may be used, with the consequence that the ODEs are linearly coupled, and more functions are needed to attain convergence.The key question in the Galerkin method is the convergence of the solution.In order to have a reasonable number of degrees of freedom, it is important to use the most significant modes.In addition to the asymmetric mode directly driven into vibration by the excitation (driven mode), it is necessary to consider (i) the orthogonal mode having the same shape and natural frequency but rotated by /(2 ) n π (companion mode), (ii) additional asymmetric modes, and (iii) axisymmetric modes.In fact, it has clearly been established that, for large-amplitude shell vibrations, the deformation of the shell involves significant axisymmetric oscillations inwards.According to these considerations, the radial displacement w is expanded by using the eigenmodes of the empty shell, which are unchanged for the completely filled shell (Amabili 2003): where kn is the number of circumferential waves, m is the number of longitudinal half-waves (only odd values are used for symmetry), m m L λ π = , and t is the time; A m,n (t), B m,n (t) and A m,0 (t) are the generalized coordinates that are unknown functions of t; the mode driven in resonance is (1, n), i.e. the mode for m = k = 1.The number of degrees of freedom used in the present numerical calculations is 16 (Amabili 2003).
The presence of pairs of modes having the same shape but different angular orientations, the first one described by cos( ) nθ (driven mode for the excitation given by equation ( 2)) and the other by sin( ) nθ (companion mode), in the periodic response of the shell leads to the appearance of a travelling wave around the shell in the θ direction when both these modes are active and when they have a relative time shift.This phenomenon is due to the axial symmetry of the system.
When the excitation has a frequency close to the resonance of a particular mode, say (m = 1, n), results for relatively low amplitude excitation (case of periodic response) show that (i) the generalized coordinates A 1,n (t) and B 1,n (t) have the same frequency as the excitation, (ii) the coordinates A 1,2n (t), B 1,2n (t), A 3,2n (t), B 3,2n (t) and all the coordinates associated with axisymmetric modes have twice the frequency of the excitation, and (iii) the coordinates A 3,n (t), B 3,n (t), A 1,3n (t), B 1,3n (t), A 3,3n (t) and B 3,3n (t) have three times the frequency of the excitation.
Expansion (6) used for the radial displacement w satisfies identically the out-of-plane boundary conditions and the continuity of the circumferential displacement.The boundary conditions for the in-plane displacements give very complex expressions when transformed into equations involving w.Therefore, they are modified into simpler integral expressions that satisfy boundary conditions on the average (Amabili 2003).
When the expansion of w, equations ( 6), is substituted in the right-hand side of equation ( 3), a partial differential equation for the stress function F is obtained, composed of the homogeneous and the particular solution.
Equations ( 4) and equation ( 6) present the same spatial distribution on the shell surface.Therefore, fluid pressure gives only an inertial effect, which is different for each mode of the expansion.Hence, the fluid is expected to change the nonlinear behaviour of the fluid-filled shell, as a consequence of the fundamental interaction among asymmetric and the axisymmetric modes.Usually the inertial effect of the fluid is larger for axisymmetric modes, thus enhancing the nonlinear softening type behaviour of the shell.
By use of the Galerkin method, 16 second-order, ordinary, coupled nonlinear differential equations are obtained for the variables A m,kn (t), B m,kn (t) and A m,0 (t), for m = 1,… M and k = 1, …3, by successively weighting the original equation (1) with 3 Copyright © 2006 by ASME the functions that describe the shape of the modes retained in equation ( 6).These equations have very long expressions containing quadratic and cubic nonlinear terms and have been obtained by using the Mathematica 4 computer software (Wolfram 1999), in order to perform analytical integrals of trigonometric functions.The generic jth Lagrange equation is divided by the modal mass associated with , taking the following form where ˆj f is projection of the nondimensionalized force, which must be set equal to zero in all the equations where are coefficients of quadratic terms and are coefficients of cubic terms; ζ ( ) j i pk h j is the modal damping ratio, substituting here the unrealistic viscous damping introduced in equation ( 1).In equation ( 7) each generalized coordinate q j (and therefore the modal damping j ζ and circular frequency j ω ) has to be referred to mode (m, n), i.e. q j = A m,n or B m,n .For computational convenience a non-dimensionalization of variables is also performed: the time is divided by the period of the resonant mode and the vibration amplitudes are divided by the shell thickness h.It can be observed that nonlinear terms do not involve time derivative of q j .By introducing a dummy variable, the K second-order equations are transformed into 2×K first-order nonlinear differential equations that are studied by using (i) the software AUTO 97 (Doedel et al. 1998) for continuation and bifurcation analysis of nonlinear ordinary differential equations, and (ii) direct integration of the equations of motion by using the DIVPAG routine of the Fortran library IMSL.Continuation methods allow following the solution path, with the advantage that unstable solutions can also be obtained; these are not ordinarily attainable by using direct numerical integration.The software AUTO 97 is capable of continuation of the solution, bifurcation analysis and branch switching by using pseudo-arclength continuation and collocation methods.

PROPER ORTHOGONAL DECOMPOSITION (POD) METHOD
The POD method optimally extracts the spatial information necessary to characterize the spatio-temporal complexity and inherent dimension of a system, from a set of temporal snapshots of the response, gathered from either numerical simulations or experimental data.In the present context, the temporal responses are obtained via the conventional Galerkin solution.The proper orthogonal modes obtained by the POD method will be used as a basis in conjunction with the Galerkin approach.The solution can be expressed by using the basis of the proper orthogonal modes (ξ where a i are the proper orthogonal coordinates and K is the number of proper orthogonal modes (degrees of freedom) used to build the POD model (in general, significantly lower than K in equation ( 5) necessary for the conventional Galerkin method).
The displacement field w is divided into its time-mean value ( ) w ξ and the zero-mean response ( ) In the POD method, the proper orthogonal modes are obtained by minimizing the objective function Ω ∀ ∈ ξ with denoting the time-averaging operation and (ξ) ψ the generic POD mode.Minimizing of the objective function ( 9) is obtained, after some mathematics, by solving the following eigenvalue problem: where is the time-averaged spatial autocorrelation function.
A Galerkin projection scheme for determining proper orthogonal modes (POMs) semi-analytically, and in parallel to approximate the solution of the PDE has been developed in (Sarkar and Païdoussis 2003).
The optimal number of terms K to be retained can be estimated by 1 1 0.99 in equation (10); however, for each problem this cut-off value can be different.It would be useful to check the convergence of the solution by increasing the value ; over a certain value, the results can become less accurate, because the additional terms introduced in the expansion may be highly noise-polluted.

K
In some applications, it may be better to use time responses obtained for different system parameters in order to produce better proper orthogonal modes.
The POD method gives a set of equations of motion with exactly the same structure of those obtained with the conventional Galerkin method, equation (7), with the difference of a reduced order of degrees of freedom, i.e.K K < .Both travelling waves around the shell in opposite directions have been taken into account (one obtained by direct integration and the other by changing the sign to the generalized coordinates associated to sin(n θ) terms in equation ( 6)) to construct the POD reduced-order model; this is fundamental to reproduce the axisymmetry of the shell with no preferential direction.

4
Copyright © 2006 by ASME 3.3 ASYMPTOTIC NONLINEAR NORMAL MODES (NNMS) METHOD Nonlinear Normal Modes (NNMs) are here defined as invariant manifold of the state space.They are moreover chosen tangent at the origin, which corresponds to the position of the structure at rest.An asymptotic procedure, based on the normal form theory, is used to compute the NNMs of the system.The method is here briefly recalled, the interested reader is referred to Touzé and Amabili (2005) for a complete description.In particular, it allows to take modal damping into account in the derivations, hence extending previous results obtained for conservative systems (Touzé et al. 2004).A thirdorder asymptotic development is applied, in order to perform a nonlinear change of coordinates for the system of damped unforced equation of motion, corresponding to equation ( 7) with right-hand-side equal to zero.A real formulation is used, so that normal forms are expressed with oscillators.The dummy variable j j y q = for the nondimensional velocity permits to recast the system of equations into first-order.The nonlinear change of coordinates is:

∑∑∑
(11) where r j is the transformed nondimensional displacement and s j is the transformed nondimensional velocity; other symbols are the transformation coefficients, known by analytical expressions.After substitution of ( 11) into (7), the dynamics, written with the newly introduced variables (r j , s j ), is expressed in an invariant-based span of the state space by cancelling all the invariant-breaking terms (i.e. the non-resonant terms).As a result, proper truncations can now be realized, as all invariantbreaking terms between oscillators in Eqs.(7), have been cancelled.The reduction can now be applied by simply selecting the most important normal coordinates for simulation (master coordinates), and cancelling all the others.In the case considered here, the minimum model must retain the NNMs corresponding to the driven mode (r 1 , s 1 , that are the continuations of A 1,5 and 1,5 A ) and the companion mode (r 2 , s 2 , that are the continuations of B 1,5 and 1,5 B ), as these two modes have the same eigenfrequency (1:1 internal resonance).Finally, the reduced-order model (ROM) built by selecting these two pairs of coordinates writes: where are the coefficients of cubic terms in equation ( 7), ( ) where are the coefficients of quadratic terms in equation ( 7).Forcing term has been added at the end of the process and now appears in Eq. ( 12a).This is the second approximation used for building this ROM, as a time-invariant manifold is used.The most accurate solution would have consisted in constructing a periodically forced invariant manifold, see e.g.(Jiang et al. 2005b).However, this results in a very complicated formulation and time-consuming numerical calculations for constructing the ROM.The proposed method has the advantage of the simplicity, the quickness of the computation, and allows deriving a differential model that could be used easily for parametric studies.However, it is valid, strictly speaking, only for small values of ( ) j il g f .
With the NNMs method, the original 16 degree of freedom of the conventional Galerkin model have been reduced to two in equation ( 12).However, differently to POD method, the structure of the equations of motion is changed.In fact, quadratic nonlinear terms have been cancelled, but cubic terms involves both the transformed nondimensional displacement and the transformed nondimensional velocity.

5
Copyright © 2006 by ASME 3.4 DISCUSSION After presentation of the two reduction methods, a first discussion on their theoretical settings is here provided in order to underline their abilities and limitations.
The POD method, which consists in finding the best orthogonal hyper-planes that contains most information, is essentially a linear method.This can be seen as an advantage since few manipulations, involving linear algebra only, are needed to construct the ROM.The key formula of the method, Eq. ( 14), is an eigenvalue problem.On the other hand, the linear essence of the method may be a drawback, as curved subspaces are generally more suitable to capture clouds of points with complicated shapes.A NNM, being an invariant manifold in state space, is a curved subspace, so that the NNMs reduction method is essentially nonlinear.The invariance property is the key that allows finding the lowest dimensional subspaces that contains dynamical properties, since dynamical motions do not stay within any other subspace that do not share this invariance property.This is the main advantage of the NNMs method as compared to the POD.It is expected that fewer NNM are necessary than POD modes.This will be illustrated in sections.
The POD method is global in the sense that it is able to capture any motion in state space and furnishes the adapted basis for decomposing it.This is an advantage as compared to the asymptotic NNMs method used here, which relies on a local theory.The third-order development, Eq. ( 11), is valid for small values of the modal amplitudes.Moreover, the use of a time-invariant manifold prevents for reliable results for large values of the amplitude of external forcing.When increasing the non-linearities by feeding more energy into the system, the results provided by the asymptotic NNMs method are awaited to deteriorate.This is not the case for the POD, if one has taken care to construct its POD-based ROM with clouds of points that are significant for a large range of values of the nonlinearity.In this context, it has already been argued (Kerschen 2003, Amabili et al. 2006) that a chaotic response is the best candidate for building the POD.
Finally, the two methods differ radically in the way the ROM is built.For the POD, one must mandatory have a response of the system to be able to build the ROM.In the present context of a completely theoretical model, this is a drawback since one must compute time responses to be in position of reducing.Moreover, it has been underlined by Amabili et al. (2003Amabili et al. ( , 2006) that the choice of these time responses is not an easy task that could not be done blindly.By comparison, the asymptotic NNMs method do not need any response of the system, but dynamical properties only, that are here provided by Eq. ( 7), i.e. after projection of the PDE with the Galerkin method.With at hand the eigenvalues of the linear part (eigenfrequencies ω j and damping coefficients ζ j ) and the nonlinear coefficients ( ) j il g and , the nonlinear change of coordinates, Eq. ( 11), can be applied directly to obtain the ROM.As the coefficients in Eq. ( 11) are computed once and for all, application of the method is easy and do not ask for computation time.In the next section, all these conclusions will be illustrated with the numerical results.ω ω = 0.9714 and at 1.0018, where branch "2" appears.
This new branch corresponds to participation of both 1, ( ) and , giving a travelling wave response.Branch "2" undergoes two Neimark-Sacker (torus) bifurcations (TR), at 1, ( ) / n ω ω = 0.9716 and 0.9949.Amplitude-modulated (quasi- periodic) response is indicated in Figure 1 for 1, ( ) The quasi-periodic time response of the shell for excitation of 3 N at frequency 1, 0.991 n ω ω = , branch "2", is reported in Figure 2 for the most significant generalized coordinates.This time response, which is more suitable than simple periodic responses to construct accurate POD models (Amabili et al. 2003), has been used to build a POD model.Both travelling waves around the shell in opposite directions have been taken into account (one obtained by direct integration and the other by changing the sign to the generalized coordinates associated to sin(n θ) terms in equation ( 8)) to construct the POD reducedorder model.The optimal number of proper orthogonal modes K to be retained in the reduced-order model can be estimated by plotting 1 ∑ ∑ as a function of K ; three proper orthogonal modes absorb practically all of the shell energy for the response at 1, / n ω ω = 0.991; therefore three proper orthogonal modes are used in the POD model.The first proper orthogonal mode is the driven mode, the second is the companion mode and the third is the axisymmetric mode.
Responses obtained by using the conventional Galerkin method (16 dofs) and by using the POD method (3 dofs) compare very well for excitation of 3 N, as shown in Figure 3 for both driven and companion modes; the main difference is a slight shift on the right of the first bifurcation point of branch "1".It can also be observed that the natural frequency computed with the POD model is practically identical to the one computed with the Galerkin model.Figure 3 also shows the response computed with the NNMs method with only two degrees of freedom.It can be observed that also the response computed with the NNMs method compare very well with the original Galerkin model, with the curves just very slightly shifted in the left and with exact qualitative behaviour.In this case the maximum vibration amplitudes reach about 1.5 h for the driven mode and 0.9 h for the companion mode.
In order to compare results also in time domain, the quasiperiodic responses ( 1, / n ω ω = 0.991, excitation 3 N) computed with the POD and NNMs models are reported in Figures 4 and  5, respectively; this response is more critical to be reproduced by reduced-order models than simple periodic responses.Whereas the response computed with the POD is in reasonably good agreement with the one in Figure 2 obtained with the Galerkin model, the response calculated by using the NNMs model is practically coincident with this.
Figure 6 has been obtained with the same three models (Galerkin, NNMs and POD for which the same equations obtained with time response for excitation of 3 N has been used), but for excitation of 8 N.In this case the maximum vibration amplitudes reach about 3 h for the driven mode and 2.5 h for the companion mode; differences among the three models become much more significant than in the previous case.In particular, the POD model is relatively close to the original Galerkin model, the main difference being the first bifurcation point of branch "1", which is now significantly shifted on the right giving rise to a significant difference in the qualitative behaviour of the two models.The NNMs model has qualitatively the same behaviour of the original Galerkin model, but the response is significantly shifted on the left giving rise to model overestimating the softening nonlinearity of the system.It can be observed here that the POD model could be improved by using a time response computed for excitation of 8 N to find the proper orthogonal modes; however it is interesting here to investigate the robustness of a reduced order model to changes in the system parameters, therefore it is convenient to use the same model.On the other hand, the NNMs model is built once and for all, and may not been changed when varying the amplitude of the forcing.The observed differences with the reference solution are the consequences of the two approximations used to build it: asymptotic development and time-invariant manifold.

CONCLUSIONS
Results show that both the proper orthogonal decomposition (POD) and the nonlinear normal modes (NNMs) methods are suitable to build reduced order model of a water-filled shell.In particular, a larger reduction of the model is possible by using the NNMs method.However, the asymptotic formulation used here in the NNMs method does not make it suitable to study very large vibration amplitudes, where the POD model performs better.These results have to be related to the properties of the two methods.The nonlinear character of the invariant manifolds that defines NNMs allows better reduction than the POD method, which is a linear decomposition.On the other hand, the global nature of the POD provides very robust results even for complex dynamics..However, it has been found that, although being an approximation, the qualitative behaviour of the NNMs compares very well with the original solution, until the validity limits are attained.It has been found numerically that these limits are not that small: amplitude of vibration up to 3h.
Construction of the NNM-based reduced-order model with the asymptotic method is direct and does not need for intensive computations, as a single nonlinear change of coordinates, computed once and for all, is required.The method can be thus blindly applied, provided the Galerkin projection has been performed on a large number of modes.On the other hand, for the POD, particular care must be placed in the choice of the time responses used to build it, as already discussed in Amabili et al. (2003).
To conclude, the investigations conducted here show that for moderate vibration amplitude, the asymptotic NNMs method provides more reduced equations that always recover the qualitative behaviour.However, the method must be modified in order to bypass its main limitation, which is due to the asymptotic development.Unfortunately, only numerical solutions are possible, thus leading to an intense numerical effort in order to build the reduced-order model.Hence, for very large vibration amplitude and large range of parameter variations, the POD method still performs better due to its global nature.