On the convergence of nonlinear modes of a finite element model

Convergence of finite element models is generally realized via observation of mesh independence. In linear systems invariance of linear modes to further mesh refinement is often used to assess mesh independence. These linear models are, however, often coupled with nonlinear elements such as CFD models, nonlinear control systems, or joint dynamics. The introduction of a single nonlinear element can significantly alter the degree of mesh refinement necessary for sufficient model accuracy. Application of nonlinear modal analysis [1,2] illustrates that using linear modal convergence as a measure of mesh quality in the presence of nonlinearities is inadequate. The convergence of the nonlinear normal modes of a simply supported beam modeled using finite elements is examined. A comparison is made to the solution of Boivin, Pierre, and Shaw [3]. Both methods suffer from the need for convergence in power series approximations. However, the finite element modeling method introduces the additional concern of mesh independence, even when the meshing the linear part of the model unless p-type elements are used [4]. The importance of moving to a finite element approach for nonlinear modal analysis is the ability to solve problems of a more complex geometry for which no closed form solution exists. This case study demonstrates that a finite element model solution converges nearly as well as a continuous solution, and presents rough guidelines for the number of expansion terms and elements needed for various levels of solution accuracy. It also demonstrates that modal convergence occurs significantly more slowly in the nonlinear model than in the corresponding linear model. This illustrates that convergence of linear modes may be an inadequate measure of mesh independence when even a small part of a model is nonlinear.

ABSTRACT.The convergence of the nonlinear normal modes of a simply supported beam modeled using finite elements is studied in this paper.A comparison is made to a solution via a Galerkin approach.Both methods suffer from the need for convergence in power series approximations, however, the finite element modeling method introduces the additional concern of mesh independence.The importance of moving to a finite element approach is the ability to solve problems of a more complex geometry for which no closed form solution exists.This case study demonstrates that a finite element model solution does converge nearly as well as a continuous solution, and presents rough guidelines for the numbers of expansion terms and number of elements needed for various levels of solution accuracy.Vibrations of linear systems are governed by a set of simultaneous linear ordinary differential equations.A differential equation is linear if the state variables and their derivatives appear only to the first or zeroth power and are not multiplied by any other state variable.The transient and steady state response of a linear vibrating system can be ex-pressed as the sum of the system's response to the initial conditions with that of the response to the input force taken separately.The equations of motion are such that a linear combination of inputs results in the same linear combination of outputs.If this is not so then the system is nonlinear.The solutions of linear differential equations are obtained most conveniently by means of linear transformations rendering the set of equations independent.The concept of normal modes and of eigenvalues is well defined for such systems.The term normal solution is associated with a fundamental set of solutions for linear systems and a linear combination of the normal solutions yields all solutions of the system.For linear, conservative, non-gyroscopic, spatially distributed systems each normal mode has associated with it a mode shape, mathematically represented by an eigenfunction, and a natural frequency represented by the corresponding eigenvalue.These normal modes represent a special set of motions of the system in which it behaves like a system of lower order.For virtually all linear systems, the dynamics of each individual mode is governed by a second order linear modal oscillator which is uncoupled from all other modal oscillators.
The study of the modal vibrations of nonlinear systems with many degrees of freedom is concerned with the search for some or all periodic solutions for systems of nonlinear differential equations.Various disciplines within mathematics are commonly used to deduce partial results that generally can be categorized as systems that are weakly nonlinear and those that are strongly nonlinear.Meirovitch in [9] defines normal modes in terms of eigenfunctions and the expression of an arbitrary system response as a superposition of modal responses.The concept of normal modes of motion is well developed for a wide class of dynamic systems that involve linearizing the nonlinear equations of motion about an equilibrium configuration to obtain linear differential equations [2].The concept of normal modes for n-degree of freedom nonlinear systems was first defined in a general manner by Rosenberg who in [1] deals with the existence and stability of normal modes for a system of n masses interconnected by nonlinear symmetric springs, and having n degrees of freedom.Shaw and Pierre in [3] generalize the definition of normal modes for nonlinear systems to include a very wide class of systems that include nonconservative, gyroscopic and infinite dimensional systems.In [1 0, 11] Slater classifies the types of nonlinear normal modes by definitions.Methodologies developed by Shaw and Pierre [3][4][5][6] and Shaw, Pierre and Boivin [7] illustrate the use of the method of invariant manifolds to obtain normal modes for weakly nonlinear systems, and demonstrate means of generating differential equations of motion that govern the dynamics of a system undergoing a nonlinear mode motion.These methods use procedures from the center manifold theory by Carr [8].To date, two different methodologies based on invariant manifold methodologies have been developed to determine the nonlinear modes of systems.The first one treats (continuous) systems as such and obtains the nonlinear normal modes directly from the differential equations of motion [4) and is referred to as the continuous nonlinear normal mode method, while the second one first discretizes the system's equations of motion using its linear normal modes, and then applies the theory developed for nonlinear discrete systems and is referred to as the discretized nonlinear normal mode method.In [7), Shaw, Pierre and Boivin show the continuous approach to be less reliable than the discretized approach.The task of presenting a comprehensive survey of the various methodologies and results is beyond the scope of this relatively brief contribution.
Here the method presented by Shaw and Pierre [3][4][5][6] is used to determine the mode shapes of nonlinear vibratory systems.This method is fundamentally different from others developed to deal with nonlinear systems in the past.It exploits the underlying geometrical structure which exists in the phase space for typical nonlinear distributed parameter systems near an equilibrium that comprises a set of invariant manifolds.These manifolds provide a means of constructing mode shapes and modal oscillators for such nonlinear systems about equilibrium.In other words, rather than view nonlinear normal modes as synchronous motions of conservative systems, they are viewed as motions on invariant manifolds that are tangent to, and of the same dimension as, the linear eigenspaces in the system phase space [3].
The usefulness of the concept of normal mode motion is greatly impaired in nonlinear systems because the principle of superposition of solutions no longer holds true for such systems.It is very difficult, if not impossible, to find closed form solutions to most nonlinear vibration problems.While a given nonlinear system may be approximated successfully by a linear system, more often than not solutions so obtained are very restricted.When a linear modal analysis of a nonlinear system is performed, typically, the resulting set of nonlinear equations is truncated to retain a small number of linear modes.The resulting reduced order model may be inaccurate or as in the case of internal resonances, even incorrect, due to the loss of nonlinear interactions between the modeled and unmodeled linear modes (a result of the linear assumption).The nonlinear normal modes defined here are invariant subspaces for the nonlinear equations of motion.This means that if the system is given an initial condition in S, the solution of the governing equations of motion remains in S for all time.The definition of nonlinear normal modes of vibration as motions occurring on invariant manifolds allows the effects of several linear modes to be incorporated into one nonlinear normal mode.When generalized this definition would allow the determination of multimode invariant manifolds which would incorporate the effects of several nonlinear modes, allowing interactions between them [3].
Difficulties in obtaining closed-form solutions to problems of continuous (distributed parameter) systems can be attributed to the inherent difficulty of solving partial differential equations, with or without space dependent coefficients, and in satisfying the boundary conditions.These difficulties can be circumvented by eliminating the spatial dependence from the problems by discretizing them.The finite element method has been viewed as a general discretization procedure of continuum problems posed by mathematically defined statements.All such analyses follow a standard pattern which is universally adaptable to discrete systems.The finite element method may be defined as a method of approximation to continuum problems such that the continuum is divided into a finite number of parts called elements, the behavior of which is specified by a finite number of parameters and the solution of the complete system as an assembly of its elements follows precisely the same rules as those applicable to the standard discrete problems.
Another way of doing this is by expanding the solution as a finite series of given basis functions.Discretization in such a manner essentially transforms vibration problems described by partial differential equations into problems described by sets of simultaneous ordinary differential equations.Consequently it reduces the associated eigenvalue problems from differential ones to algebraic ones.One of the methods used to discretize a system is the Galerkin method which is basically a weighted residuals method, and thus applicable to both self-adjoint and non self-adjoint systems.This paper applies the finite element method to the formulation developed by Shaw and Pierre [3][4][5][6] and Boivin [7], in an attempt to demonstrate the convergence of finite element analysis for continuous nonlinear problems when used in conjunction with the invariant manifold approach and to develop an alternative approach to nonlinear finite element methodology.

FINITE ELEMENT APPROACH
The solution approach is detailed using a simply-supported, linear Euler-Bernoulli beam of a constant cross section with a nonlinear cubic spring attached to its middle as shown in Figure 1 as an example.If the spring is chosen to be purely cubic and the beam is assumed to deform in the linear range, then the normal modes of the linearized system will be the same as those of a simply supported beam (i.e.pure sine waves).The example here is selected to illustrate the method and the types of results that can be obtained.Most of the symbolic manipulation and almost all of the numerical work for the examples were done using Mathematica 2.2.The plots were generated using Matlab 4.2c.If the beam is of unit length, the continuous equation of transverse motion of the system is, where E is the Young's modulus of the beam, I is its second moment of area, m is its mass per unit length, y is the nonlinear stiffness of the spring, s is the abscissa along the beam, u(s,t) is the transverse deflection, • s denotes a derivative with respect to s (spatial derivative), an overdot represents a temporal derivative (with respect to time), and 8(.) is the Dirac function.Dividing the equations of motion (6.1) by m reduces this to a non-dimensional form, ii+au..,_, + {3u 3 where a= Ellm and f3 = ylm.The associated boundary conditions are: We now apply the finite element method to the above problem.The beam is divided into a number of discrete elements (finite elements) and the finite element method is utilized to obtain the equations of motion of the now discretized system.The aim is to show that, with a reasonable number of elements, one can approximate the continuous system with a fair degree of accuracy.The method is first formulated using two finite elements, increasing the number of elements in subsequent formulations.In general, the procedure will be to successively discretize the beam into a larger number of elements with each successive iteration, to establish a general pattern of convergence.The first iteration using two elements is developed and explained here.
Four, eight and sixteen element formulations are developed in a similar manner.

Two-element formulation
After boundary conditions have been applied the beam model has four degrees of freedom.Nodes 1 and 3 have one rotational degree of freedom each and no translational degrees of freedom because of the boundary constraints.
l/2 l/2 The mid-beam node, 2, has both rotational and translational degrees of freedom and is also shared by the nonlinear spring.Element (3) is the nonlinear spring element that possesses no rotational degrees of freedom.Consequently node 4 has no degrees of freedom because of the boundary constraint on the nonlinear spring.

Finite element solution
The matrix form of the equation of vibration, in physical coordinates, { u}, is written as, where [M] and [K] are the mass and stiffness matrices of the beam.
The eigenvectors and eigenvalues of the normalized linear system, assuming a = 1, are the eigenvalues and eigenvectors of the matrix, In this case the eigenvectors are I -1 I l 0 -0.318385 0.0546485 and the eigenvalues of the system are, [A]= diag{1920., 40320.,98.1795,12130.7}(7) Incorporating the nonlinearity of the beam into the system changes the global stiffness matrix of the beam to u, u~ u, u.

4(
The mass matrix remains unchanged.
Introducing the coordinate transformation, where { r(t)}, is the column vector of modal coordinates, and [S] is the transformation matrix that consists of the eigenvectors of the linear problem.
This yields the equation of motion written in terms of the modal coordinates, which is premultiplied by [Sf to get, For a purely linear problem, this last step serves to decouple the equations of motion.Introduction of the nonlinear stiffness term into the system however means that the resulting equations remain coupled via the nonlinear terms.We therefore, get the following four equations: This set of coupled, nonlinear ODE's describes the dynamics of the linear modal components of the motion.The first two equations that are now decoupled represent the modal equations for the even modes two and four.The solution of the decoupled equations is straightforward and is not dealt with here.The last two equations remain coupled in the nonlinear terms and represent the modal equations for modes one and three.The nonlinear terms represent nonlinear modal coupling that exists when a traditional linear analysis is performed.These equations can be numerically integrated using a fourth order accurate Runge-Kutta time marching scheme, to obtain what can be considered the 'exact' dynamics of the system.This may be computationally very demanding since the more accurate the desired solution, the larger the number of time steps.Notice that as stated earlier only the symmetric (odd numbered) modes are affected by the nonlinearity and that only the symmetric linear modes contribute to the symmetric nonlinear modes.This is because the spring is located at a nodal point of the anti-symmetric (even) modes and thus does not affect them.Therefore the anti-symmetric modes of the system will be the same as those of the linearized system.The symmetric modes will be influenced by the nonlinear spring and, furthermore their mode shapes will be comprised only of nonlinear combinations of the symmetric linear mode shapes.
One can obtain approximations of the nonlinear modes of the system by applying to this set of nonlinear ODE's the general method developed by Shaw and Pierre in [3][4][5][6].To obtain the first nonlinear mode shape we assume r:.(the first linear modal displacement) to be the modal coordinate.Note that s 3 is the first linear modal velocity and s 4 is the third linear modal velocity.A normal mode motion is assumed by requ1ring that r:,, r 4 , y 3 and y 4 are related as follows: r, =u where the coefficients a 1 , h 1o a 2 , b 2 , .•. are to be determined.
To do so, we first take the time derivative of R 4 , and 5 4 , using the chain rule, so that, The equations of motion are also expressed in u and v as, s. = ;:.

(18)
The R 4 term in the left side of equation ( 18) is equated to the equations of motion for S 4 from equation (13) and the s. term in the left side of equation ( 18) is equated to ;: 4 from equation ( 17) respectively, and in the right sides of equation ( 18) u and ~• are replaced by v and ;: 3 from equation ( 16) respectively.This yields the required conditions for a modal motion.Equation ( 18) is then expanded out to cubic order in u and v, and the coefficients of u v u 2 uv v 2 u 3 tlv, uv 2 , and v 3 are gathered together to pr~vide 'the' e~ua~ !ions for the a 1 and b; coefficients.
(The right side is subtracted from the equation in order to give equations that are equal to zero).There are two equations for each term, one from each of the two equations in (18).The solution of these equations is fairly straightforward, with most of them being linear in the unknowns.
We can now express r 4 in terms of the modal coordinate r 3 .Using (13) and the results obtained for the coefficients a~> bJ.In order to convert the modal coordinates back to physical coordinates we use equation (9).Therefore, in physical coordinates we have, {-r 3 -0.82915ri 1 +0.000434734r 3 r?L {-0.318385r 3 +0.0453118rf-0.0000237575rJ"3 2 ), {0), {r 3 + 0.82915r:J'-0.000434734r 1 r 3 2 } (20) This gives us u,,u 2 ,u,,and u 4 in physical coordinates.Knowing these four nodal displacements we can find the displacement of the whole beam, and thus the mode shape, by using the shape functions of the beam.These shape functions are used to define the displacement at any point x on the beam element, given the displacements and rotations at the nodes i and}.
The results in the time domain can also be obtained fairly easily.One of the methods of doing so is to simply use a numerical integration routine (such as a Runge-Kutta fourth order method) over the time domain of interest.The temporal solution for the first mode is obtained by substituting (19) into (12) to yield a second order differential equation in r 3 .
It may be noted here that the application of a time marching scheme such as the Runge-Kutta method to this ODE is now trivial when compared to the initial problem.

RESULTS
Several simulations were carried out using the finite element method and the invariant manifold approach.The first and third mode shapes were obtained with the methods described in the previous section using two, four, eight and sixteen elements.The results clearly illustrate the convergence of the method with the continuous solution as calculated by Boivin, Shaw, and Pierre [7] as the number of elements is increased.It may be noted here that the convergence of the continuous solution to the exact solution, computed by direct integration of the set of nonlinear, coupled, first-order ODE's has already been demonstrated in [7].The solution of [7] will be referred to as the continuous solution while the method developed herein using the finite element approach will be referred to as the finite or discrete solution.
The results (see section 5) show that a fair approximation of the beam dynamics of the first mode at low amplitudes is obtained by using just two elements.On increasing the amplitude the nonlinearity becomes more evident and the 2 element solution deviates from the continuous solution.By increasing the number of elements the beam dynamics can be obtained with better accuracy.However, at least four elements are required to somewhat capture the third mode and eight elements to show acceptable agreement.
The error between the continuous and finite element methods is depicted as a function of the amplitude along the length of the beam.It is clear from these plots that the error increases exponentially as the amplitude approaches 1.On the other hand the error can be expected to remain constant for the linear system.This emphasizes the dependence of the nonlinear solution on the amplitude.Comparison of the linear finite element solution to the linear continuous solution also establishes the fact that when nonlinearities are neglected the results for the nonlinear systems reduce to the modal analysis for the linearized system, even though the approach is entirely different from the traditional linear normal mode analysis technique.A quick comparison of the error plots also shows how the magnitude of error decreases with an increase in the number of elements.This is true for the first as well as the third mode shapes.Also note that the first mode converges earlier than the third mode.
An important characteristic of nonlinear modes, namely the effect of higher linear modes on the nonlinear modes can be observed here.This is true in the case of the first as well as the third nonlinear modes.For the third mode the two element solution appears to be more accurate than the four element solution.However closer examination of the results proves this to be a coincidence.At low amplitudes the four element formulation is more accurate than the two element solution.However the effect of the fifth linear mode on the third nonlinear mode becomes more pronounced as the amplitude is increased and at higher amplitudes the third nonlinear mode shows significant contributions from the fifth linear mode.However the eight (and sixteen) element for-mulation does not show this contribution, and neither does the continuous solution indicating that the four element mesh is not fine enough.The two element system, on the other hand does not include the effects of the fifth linear mode because it only has a total of four degrees of freedom after the application of boundary conditions.It may be noted that even for linear systems one does not use a four degree of freedom system to find the third mode, the general for rule of thumb for linear modal analysis using numerical methods being

q= (2p)
where q is the number of degrees of freedom and p is the number of eigenvalues and eigenvectors to be calculated.
From the example considered here, it is clear that for nonlinear systems more degrees of freedom need to be used (than a linear system) to successfully capture the dynamics of the system.At this stage though it is not appropriate to form a rule of thumb (on the basis of one example alone).Directions for future work include modeling systems spanning a variety of geometries and nonlinearities in order to formulate such a rule of thumb for nonlinear systems.
Table 8.1 shows the results of a convergence study conducted on the example beam considered here.
Convergence is studied for 5% accuracy in the mode shape.Note that for a linear system the convergence is not affected by a change in the amplitude.

CONCLUSION
By considering nonlinear normal modes of vibration of nonlinear systems as motions occurring on invariant manifolds the effects of several linear normal modes can incorporated into one nonlinear normal mode.For a pure nonlinear modal motion, the invariant manifold approach achieves the same accuracy as that obtained using several linear normal modes but with significantly reduced computational cost.
The finite element method serves to break a problem down to solving a set of simple linear equations.The importance of moving to a finite element approach is the ability to solve problems of a more complex geometry for which no closed form solution exists.The accuracy of the finite element method in solving the system equations was investigated by applying the methodology to determine the nonlinear mode shapes of an Euler Bernoulli beam constrained by a nonlinear spring.The convergence of the finite element method when used in combination with the invariant manifold technique to determine the amplitude-dependent mode shapes of nonlinear vibratory systems was demonstrated by obtaining the mode shapes of the beam for various meshes (number of elements).

FIGURES
Cross sectional area (of a beam) Dirac function Youngs modulus Nonlinear stiffness coefficient Identity matrix Area moment of inertia Matrix of eigenvalues Length of the beam Mass, Damping, Stiffness matrices Mass of body i Material density Matrix of modeshapes (eigenvectors) Physical coordinates (state-space notation) Modal coordinates (state-space notation)

Figure 2 :
Figure 2: Finite element model of the beam