Large-amplitude non-linear normal modes of piecewise linear systems

A numerical method for constructing non-linear normal modes (NNMs) for piecewise linear autonomous systems is presented. These NNMs are based on the concept of invariant manifolds, and are obtained using a Galerkin-based solution of the invariant manifold’s non-linear partial differential equations. The accuracy of the constructed non-linear modes is checked by the comparison of the motion on the invariant manifold to the exact solution, in both time and frequency domains. It is found that thisconstruction approach can accurately capture the NNMs over a wide range of amplitudes, including those with strong non-linear effects. Several interesting dynamic characteristics of the non-linear modal motion are found and compared to those of linear modes. A two-degree-of-freedom example is used to illustrate the technique. The existence, stability and bifurcations of the NNMs for this example are investigated.


Introduction
The concept of non-linear normal modes was originated by Rosenberg [1][2][3], based on the analysis of discrete, symmetric systems with smooth non-linearities. Since then, the existence [4,5], stability [6][7][8][9][10], and construction [11][12][13][14] of non-linear normal modes have been among the topics of investigation in this field. More recently, an alternative definition for NNMs was introduced by Shaw and Pierre [15][16][17], based on invariant manifolds. With asymptotic expansions, the NNMs can be constructed symbolically [16,17], but are accurate only in a neighborhood of the original equilibrium position. A more recent approach [18] extends the construction aspects of the invariant manifold approach to strongly non-linear regions by using a Galerkin projection method to solve the invariant manifold equations. This makes the accurate construction of NNMs possible for a wide range of non-linear dynamical systems.
Previous work on non-linear normal modes dealt primarily with systems with smooth nonlinearities. However, many engineering systems involve components with contact, clearance, or different elastic materials. Such systems are often conveniently modelled by equations of motion with piecewise linear (PWL) terms. Due to their practical importance, the dynamic behavior of piecewise linear systems has been the subject of many investigations, which mainly focus on forced response under periodic excitation [19][20][21], including dynamic behaviors such as bifurcations and chaos [22][23][24].
The non-linear modal behavior of PWL systems was first considered by Zuo and Curnier [25], based on autonomous, piecewise linear, multi-degree of freedom (d.o.f), gyroscopic and nongyroscopic systems, which are the simplest models of cracked rotating shafts and cracked beams. However, the non-linear modal motions were found by a direct approach and a more general construction method for NNMs was not pursued. Moreover, in Ref. [25] the switching hyperplane of the PWL systems considered passes through the origin, so that it is a special class of general PWL systems. Chen and Shaw [26] investigated a construction method of NNM for PWL systems based on asymptotic expansions. Since the switching hyperplane is not at the origin, the class of PWL systems in Ref. [26] is more general, but the asymptotic expansion can no longer be initiated at the static equilibrium position. The NNM is expanded in a series form in a neighborhood of an invariant disk, which makes it applicable near the switching hyperplane, but it is not valid at large amplitudes beyond the switching plane. Chati et al. [27] constructed the NNMs of a two-d.o.f. PWL system using perturbation methods. The system they considered is a simplified model of the vibrations of a cantilever beam with a transverse edge crack. The NNMs obtained are only accurate for PWL systems with a small clearance spring located at the switching hyperplane, due to the approximations of perturbation methods. They also utilized the idea of bilinear frequency to compute the natural frequencies of non-linear normal mode motions. The frequencies obtained with the bilinear formula are good approximations when the difference between the linear regions is small.
In this paper, we focus on a construction method of NNMs for the class of systems considered in Ref. [26]. The general dynamic behavior of the non-linear modal motions is also discussed using the NNM results. Since asymptotic series expansions are not naturally suited for this type of nonlinearity, the Galerkin-based approach developed in Ref. [18] is extended to PWL systems and applied to a sample problem. The NNMs constructed in this manner are accurate over a large amplitude range. Also, no specific analysis is needed to account for the switching hyperplane, along which the system changes form, thereby making the approach much less cumbersome than the expansions used in Ref. [26]. Once the NNMs are constructed, the non-linear modal dynamics for the individual NNMs can be determined.
The paper is organized as follows. The class of PWL systems studied is described in Section 2. In Section 3, the Galerkin-based approach for the construction of large-amplitude NNMs is briefly reviewed and adapted to PWL systems. A two-d.o.f. example system is illustrated in Section 4, which demonstrates the individual NNMs and their general dynamic behavior (including stability calculations). Finally, some conclusions are drawn in Section 5.

Piecewise linear systems
Following the definition of Chen and Shaw [26], the dynamic system considered here is an unforced, undamped, autonomous N-d.o.f. system of the form where d > 0 is a real scalar constant, Z; h; bAR N ; and M; K 1 ; K 2 are real symmetric positive definite N Â N matrices. Obviously, this system has two distinct linear regions separated by a hyperplane, fZAR N : h T Z ¼ dg: We denote the first region as fZAR N : h T Zodg and the second region as fZAR N : h T Z > dg: For small amplitudes, solutions remain in the first region and the response is simply linear and well understood. As the amplitude increases, the motions begin to pass through the switching hyperplane and enter the second region, in which case the response is non-linear and no longer simple [22,25]. System (1) represents a large variety of vibration systems with clearance or impact, which can be modelled by piecewise linear springs. However, there is only a single switching hyperplane in system (1). For PWL systems with more than one surface of discontinuity, the general behavior is much more complicated [28], and this is not considered here. Moreover, system (1) is conservative and non-gyroscopic.
In Eq. (1), linear modal analysis can be applied to the first region, fZAR N : h T Zodg; containing the static equilibrium point. The eigenvector matrix, Q; can be found and normalized with respect to the mass matrix M: For simplicity, it is assumed that all the eigenfrequencies of the subsystem in the first region are distinct. After a linear modal transformation, Z ¼ QZ; system (1) takes the standard form where ZAR N is the vector of modal co-ordinates, and the N Â N diagonal matrix L 1 has entries that are the squares of the small-amplitude (first region) natural frequencies. The piecewise linear force vector, f ðZÞAR N ; is given by where Hð:Þ is the heaviside function that stands for the switching hyperplane.

Non-linear normal modes and the invariant manifold approach
Following the concept of invariant manifolds [15], a non-linear normal mode is ''a family of motions which lies on a two-dimensional invariant manifold in the system phase space''. Here, the term invariant indicates that any motion initiated on the manifold will remain on it for all time. In this formulation, a single pair of displacement-velocity co-ordinates is chosen as the master coordinates, which characterize the individual non-linear modal motion that occurs on the manifold. The remaining d.o.f.s are represented by slave co-ordinates, composed of displacementvelocity pairs for those d.o.f.s, which are functionally dependent on the master co-ordinates. The dynamics of the master co-ordinates dictates the response of the slave co-ordinates through these relationships. The bulk of the work for determining NNMs lies in the determination of these constraint functions, which describe the geometry of the NNM invariant manifold in the system state space. In previous work, asymptotic series have been used to obtain approximate solutions of the manifolds, whereas here we employ a numerical solution that provides better accuracy over a much larger amplitude range.
In order to construct accurate non-linear normal modes for this class of piecewise linear systems, the Galerkin-based approach [18] is applied to the invariant manifold equations. Considering system (2), for the kth NNM the master co-ordinates are taken to be ðZ k ; ' Z k Þ: These are transformed to amplitude and phase co-ordinates, ða; fÞ; as follows: where o k is the kth natural frequency of the linear system in the first region. The master coordinates ða; fÞ have bounded domains for amplitude and phase, aA½0; a max and fA½0; 2p; respectively, which makes a Galerkin-based approach feasible [18]. On the invariant manifold, all of the slave co-ordinates are expressed as functions of the master co-ordinates as follows: These functions describe constraint relationships between the master and slave co-ordinates, thereby providing a functional form for the NNM invariant manifold. If the system is globally linear, f ðZÞ is zero in Eq. (2), and the constraint functions, P i and Q i ; are also zero. For the piecewise linear system, the constraint functions are no longer zero when the master amplitude coordinate a is sufficiently large such that the system enters the second linear region. In other words, the constraint functions, P i and Q i ; capture the geometry of the invariant manifold as the system passes between the first and second regions.
The invariant manifold equations are formulated as follows. A first order, state-space formulation of the equations of motion (2) are used, into which Eqs. (4) and (5) are substituted for the dynamic variables. The use of the chain rule on the constraint functions P i and Q i results in partial differential equations (PDE) which govern the P i 's and Q i 's. These are given by [18] These are obviously valid when the force, f AR N ; is smooth. However, one must be cautious when dealing with non-smooth non-linearities, since these may result in constraint functions that are also non-smooth, in which case some of the terms in the PDEs, Eq. (6), may not exist. It is found that for piecewise linear systems, these governing PDEs are also valid, since all the derivatives, @P i =@a; @P i =@f; @Q i =@a; and @Q i =@f; exist in each of the two regions. At the switching hyperplane, the geometry of the invariant manifold is continuous, since the force is continuous, but derivatives are not necessarily continuous. Therefore, global expansion functions for a Galerkin approach may not be well suited for the task, and a special discretization is used to construct the invariant manifold as described below.
In order to solve the PDEs, Eq. (6), a Galerkin projection is carried out over the chosen domain, aA½0; a max and fA½0; 2p: In the f direction, the constraint equations P i and Q i are periodic, and thus Fourier series are the natural choice for the expansion. Furthermore, a half-Fourier basis is sufficient for system (2), due to its conservative, non-gyroscopic nature [18]. Specifically, cosine functions are used for the position constraints P i ; and sine functions are used for the velocity constraints Q i : Because of the nature of the manifold in the a direction, the domain aA½0; a max is divided into n equal segments, defined by aA½a j ; a jþ1 ; a j ¼ ja max n ; j ¼ 0; 1; y; n: In each segment, piecewise linear functions are used as the expansion functions. Then, the unknown position and velocity constraint relations are expanded over each segment as where the C's and D's are the unknown expansion coefficients, and T l;m and U l;m are tensor products of tent functions (defined below) in the a direction and trigonometric functions in f: Hence, for a given segment, where A 1 ðaÞ ¼ a À a j a jþ1 À a j and A 2 ðaÞ ¼ a jþ1 À a a jþ1 À a j are the tent functions employed in the segment aA½a j ; a jþ1 : The expansion functions given in Eq. (7) are substituted into the PDEs given in Eq. (6), and a Galerkin projection is carried out over the chosen segment, resulting in for i ¼ 1yN; iak; p ¼ 1; 2; and q ¼ 1yN f : This results in 2ðN À 1Þ2N f non-linear algebraic equations in the C's and D's. Note that such a system of equations must be solved for each a interval, resulting in a total of n such systems of equations that must be solved to obtain the entire manifold over the desired amplitude range. In order to search numerically for these unknown coefficients, a local optimization algorithm (the Hybrid Powell method embedded in the commercial algorithm package NAG) is applied using an initial guess. Since the manifold geometry is continuous at the switching hyperplane, zero is a good initial guess for a segment crossing the switching plane. Results for subsequent amplitude intervals are obtained in a sequential manner, where the results for the C's and D's of a preceding segment are used as the initial values for the following segment. In this manner the procedure is self-starting and no complicated initial guessing algorithms are necessary. Note that each segment on which solutions are obtained is an annular strip in the state space. These strips are pieced together to form the invariant manifold. It is important to note that a system with N degrees of freedom will generally have N NNMs that are continuations of the modes of the linear system, and that each manifold is solved for individually. Bifurcations of the NNMs can lead to more NNMs than d.o.f. [10], but these cases are not considered here.
The discretization in the a direction is analogous to the finite element method. In order to ensure boundary conforming conditions at the interface between neighboring segments, one ought to run the optimization algorithm once again after the local optimization results have been obtained for each of the individual segments. In other words, the results obtained from each segment are deemed as the initial guess for the optimization over the whole region. However, since the manifold geometry does not change rapidly, simple term-by-term averaging over the interface of contiguous segments has acceptable accuracy and is applied here.
Once the constraint functions (5) are obtained over the entire domain, the non-linear modal dynamics on the invariant manifold can be reduced to a pair of first order ordinary differential equations (ODEs) expressed in terms of the master co-ordinates a and f [18], as follows: where f k depends on a and f only, since the remaining dynamic states have been replaced by the constraint relations. Solutions of this relatively simple oscillator equation capture the dynamics of the full system restricted to the NNM manifold of interest. There will be N such NNMs, each of which is obtained independently.
given by where k 1 ; k 2 ; and k 3 are the stiffnesses of the linear springs and d 0 is the distance from the static equilibrium position of m 2 to its contacting position with spring k 3 : System (11) is obviously piecewise linear and of the form under consideration.
The system is non-dimensionalized by introducing the following non-dimensional variables and parameters, where the double derivative . ðÞ now denotes d 2 ðÞ=dt 2 : Compared to the form of Eq. (1) in Section 2, the mass matrix M; and the stiffness matrices K 1 and K 2 are Linear modal analysis can be performed on the sub-system M .
, where Z ¼ 0 is the static equilibrium position. For the generalized eigenvalue problem K 1 q ¼ o 2 Mq; the two real positive eigenvalues are found to be A two-degree-of-freedom piecewise linear system. and the eigenvector matrix is given by * Q ¼ ½q 1 q 2 ; where q 1 and q 2 are the eigenvectors corresponding to the non-dimensional natural frequencies o 1 and o 2 : After applying the modal transformation Z ¼ * QZ and some manipulations, Eq. (12) is expressed in the standard form of Eq. (2), and is thus ready for the construction of its NNMs.

Case one
The non-dimensional stiffness ratios are taken to be a ¼ b ¼ 1:5 and the mass ratio is g ¼ 1:0: The non-dimensional linear natural frequencies in Eq. which has not been normalized with respect to the mass matrix.

The invariant manifold
The constraint equations that describe the two NNMs can be constructed following the procedure described in Section 3. For the first NNM, Z 1 and ' Z 1 are chosen as master co-ordinates, where Z 1 ¼ a cos f and ' Z 1 ¼ Àao 1 sin f: The constraint equations Z 2 ¼ P 2 ða; fÞ and ' Z 2 ¼ Q 2 ða; fÞ are numerically constructed over the domain aA½0; a max ; fA½0; 2p: In order to show the first NNM over a large amplitude region, the parameter a max is set equal to 40:5; and the domain ½0; a max is divided into 81 equally sized segments. In order to ensure good numerical convergence, the number of harmonic terms is taken to be N f ¼ 64: For the second NNM, the master coordinates are ðZ 2 ; ' Z 2 Þ; and the constraint equations are Z 1 ¼ P 1 ða; fÞ and ' Z 1 ¼ Q 1 ða; fÞ: The parameters set in the numerical solution algorithm are: a max ¼ 60; with 120 segments in the a direction, and N f ¼ 32: The position constraints P 2 ða; fÞ and P 1 ða; fÞ for the two NNMs are shown in Fig. 2. The geometry is flat and zero for small values of a; but is no longer planar after a crosses the switching hyperplane. For the first non-linear mode, the switching plane is at aE2:11; corresponding to a displacement of m 2 of z 2 ¼ 1:0: As the maximum amplitude a max is reached, the displacement z 2 is about 19: Hence, the Galerkin approach is applicable into amplitude regions of strong nonlinearity, where asymptotic analyses [26] are not applicable. For the second NNM, the switching amplitude is also at aE2:11; and the maximum displacement z 2 is about 28 when a max is reached.
It should be noted here that the switching position, aE2:11; is represented by a circular line inside one of the strip segments used in the solution procedure. Therefore, the precise manifold geometry near the transition plane is not caught by this coarse discretization in a direction. However, it also shows that the Galerkin approach is robust, in the sense that it is not necessary to know the switching condition in advance. For the reader with an interest in further details about the manifold characteristics near this transition, two approaches can be taken. First, one could carry out an asymptotic analysis based on the Poincar! e map at the transition plane, as done in Ref. [26]. Or, one can refine the manifold discretization near the switching hyperplane into smaller strip segments. The first approach is complicated, but theoretically interesting. On the other hand, a small mesh in the a direction is feasible and relatively simple to implement, and was carried out here. The refined region for the first NNM is shown in Fig. 3, where the region aA½2:1082; 2:5 is evenly divided into 100 segments with N f ¼ 64: As described below, the Galerkin approach correctly captures the details of the switching plane with this refined mesh.
In order to better understand the nature of the invariant manifold, the manifold geometry is also shown in the original physical co-ordinate system ðz i ; ' z i Þ: Fig. 4 displays results based on the coarse segment discretization. The mesh shows the overall geometry of the manifold and the continuous solid curve depicts a representative motion on the manifold for a given initial value. Since the example system, Eq. (12), is conservative, all motions on the individual NNM manifolds must be periodic. Thus, the search for periodic solutions is essentially that of determining a set of initial conditions which ensure periodic response. These initial conditions, which characterize individual NNMs, can be found with numerical time integration and the search can be accelerated by the evaluation of the Jacobian matrix [25,29]. This approach is an alternative method for determining NNM invariant manifolds for this class of problems.
It is interesting to note that the geometry of the first NNM shown in Fig. 4(a) has a kink near the negative z 1 axis. This phenomenon-can be explained by the time histories shown in Fig. 5. As can be seen, in each period of the first NNM response shown, there exists a time interval over which the velocity ' z 1 is nearly zero, the displacement z 1 is almost constant, the velocity ' z 2 changes its sign, and the displacement z 2 changes quite rapidly. This behavior, which is reminiscent of sticking due to dry friction, leads to the kink in the manifold geometry. Fig. 5 displays a set of time responses on the individual NNM manifolds, obtained using two different approaches. The first category of responses is obtained from the time integration of the reduced equations of motion in the master co-ordinates a and f; which describe the dynamics on the manifold, Eq. (10). Initial conditions að0Þ > 0 and fð0Þ ¼ 0 are used, and since the motion on the manifold is periodic, its amplitude is equal to að0Þ: Then, the modal responses Z i ðtÞ and ' Z i ðtÞ are obtained from aðtÞ and fðtÞ based on the master co-ordinate definition, Eq. (4), and the slave co-ordinate constraint functions, Eq. (5). Finally, the displacements z i ðtÞ and velocities ' z i ðtÞ are reassembled via the linear modal transformation. The second category of responses consists of the periodic responses simulated from the equations of motion of the original system, Eq. (12). The initial conditions for periodic motions are obtained using a numerical direct search method. One can observe in Fig. 5 that the simulations restricted to the Galerkin-based manifold match very closely those of the original system with the same energy level.
Another check of the Galerkin-based invariant manifold can be performed, and some insight into the NNM dynamics can be gained by considering the response in the frequency domain, as shown in Fig. 6. The frequency f 0 of periodic motions on the manifold is defined as the fundamental frequency associated with the basic period t 0 ; f 0 ¼ 2p=t 0 : In Fig. 6, the frequency f 0 increases rapidly after the amplitude a crosses the transition hyperplane. Then f 0 tends to a limiting value as a increases. This limit is not equal to any linear modal frequency of sub-systems in Eq. (1), but is close to the bilinear frequency defined in Ref. [27] for PWL systems with zero gap, where the transition hyperplane is at the equilibrium position. In the present system, the gap becomes irrelevant at large amplitudes, and the results approach those of Chati and Rand [27]. In  Fig. 6, note the excellent agreement between response frequencies obtained from the simulation of the system dynamics and that of the dynamics restricted to the manifold.
The comparisons carried out in both the time and frequency domains clearly demonstrate that the Galerkin-based invariant manifolds accurately represent individual NNMs over a large range of amplitudes, and that the dynamics of the individual NNMs can be accurately reduced to a single DOF, given by Eq. (10).
Some interesting features of the dynamic behavior of the NNMs of the system can be observed by examining numerical responses for various initial energy levels in the phase plane and in the configuration space. Fig. 7 depicts phase plane diagrams, closed curves correspond to periodic responses whose amplitude depends on the initial energy level. These loops are symmetric with respect to the displacement axis (due to the conservative nature of the system), but not with respect to the velocity axis (due to the asymmetry of the restoring force). In the configuration space, shown in Fig. 8, the trajectories of solutions are represented by curves, but not by straight lines as in the linear case. 1 Moreover, the displacements z 1 and z 2 do not vanish simultaneously, 2 but they reach their maximum and minimum positions at the same time. Also, note that when the amplitude is increased, a completely new modal curve is followed, which is not simply an  extension of a curve from a lower amplitude. These deviations from linear system dynamics arise from the non-symmetric, non-linear nature of the example system considered here, and were also reported in Ref. [25].

Stability of non-linear normal modes
It is well known that non-linear systems can exhibit a wide range of behaviors, including instabilities with bifurcations, chaos, etc. In order to investigate the bifurcation characteristics of the NNMs of the example system, the stability of the periodic motions on the NNM manifolds is examined as the amplitude a increases. Two methods, characteristic multipliers and Poincar! e map, are employed to explore the stability of the NNMs and some additional features of the response.
The equation of motion of the two-d.o.f. example system, Eq. (12), are expressed in standard state space form, where Y ¼ ½z 1 ' z 1 z 2 ' z 2 T and the right-hand side is given by The characteristic multipliers can be determined from the monodromy matrix M; defined by MðaÞ ¼ Fðt 0 Þ; where t 0 is the period of the motion, a is the amplitude of the periodic motion, and the matrix Fðt 0 Þ is determined from the matrix initial value problem, where Y p is the periodic solution with amplitude a; and the period t 0 is obtained while numerically searching for periodic solutions. 3 The matrix F Y ðY p Þ is the Jacobian matrix of the right-hand side of the state equation (14), which can be calculated at each time step during a numerical integration. The stability of the periodic solution can be checked by the eigenvalues, mðaÞ; of the monodromy matrix MðaÞ: For a given periodic response of this conservative system, a necessary condition for stability is that all these eigenvalues lie on or inside the unit circle in the complex plane [29]. The monodromy matrix is typically used for systems with smooth restoring forces. In this example system, the force FðY Þ is piecewise smooth on either side of transition hyperplane. Since the solution intersects the surface of discontinuity without tangency and the initial time does not correspond to a crossing of the surface of discontinuity, the monodromy matrix can still be constructed from Eq. (15) [30]. For simplicity, an explicit forward Euler method is applied here to approximate the resulting matrix, so that it is not necessary to know the switching time in advance. The accuracy of the monodromy matrix is controlled by the time step in the numerical integration.
The movement of the eigenvalues of the monodromy matrix, mðaÞ; in the complex plane as the amplitude is varied is schematized in Fig. 9 for both NNMs. The multipliers for both NNMs stay on the unit circle for small amplitudes. As a increases, a pair of complex conjugate multipliers move on the circle towards À1: At a critical value, a ¼ 6:52 for the first mode, and a ¼ 4:04 for the second mode, this pair merges at À1 and then one of those multipliers escapes the unit circle, yielding unstable behavior. For the second mode, the motion on the manifold remains unstable above the  critical amplitude. For the first mode, the pair of separated multipliers on the negative real axis merge again at À1 as the amplitude increases further, and stability is recovered for aX8:62: The results from the characteristic multipliers indicate that the system loses stability when the multipliers satisfy mða c Þ ¼ À1; where a c is the critical amplitude. This corresponds to a flip bifurcation, or a subharmonic bifurcation [29], or a period doubling bifurcation. This period doubling is illustrated in Fig. 10 by computing the FFT of the response at two amplitudes close to, and on either side of, the bifurcation point. This flip bifurcation was also found in the piecewise linear system studied by Zuo and Curnier [25], where the instability was parameterized by a stiffness parameter, similar to the stiffness ratio b in Eq. (12).
In order to explore the dynamics in the unstable region, a Poincar! e map is now employed [31]. Since there is only a single transition hyperplane in the system, it provides a natural Poincar! e section [26], defined by the two components 4 S ¼ fY AR 4 : z 2 ¼ 1; and ' z 2 > 0g; where Y is the vector of state variables in Eq. (14). The dynamics in both sub-regions are linear and analytical solutions can be locally obtained. Given the initial conditions of the exact solution for each non-linear mode, the Poincar! e mapping can thus be constructed numerically [26]. From a cross-sectional view of the Poincar! e map, the motion on the first NNM is stable until the critical amplitude, a c ¼ 6:52; is reached, and is represented by the pair of points in Fig. 11(a). As the amplitude increases beyond the critical value, the motions become quasi-periodic, represented by the loops in Fig. 11(b,c), and then chaotic, as shown in Fig. 11(d). This sequence indicates that the period doubling at the critical point is subcritical. Also, the period doubling implies that the post-bifurcation responses will not exist in a NNM manifold, since period doubling cannot occur in a planar, autonomous system. Therefore, all post-critical responses associated with flip bifurcations must include both NNMs, as they are defined here. At the second critical amplitude, the periodic motion regains stability, as shown in Fig. 11(e). For the second mode, the periodic response remains unstable after the initial bifurcation at a c ¼ 4:04; as shown in Fig. 12(b). Here the bifurcation appears to be supercritical, and another period doubling occurs, resulting in the transition from two to four points. At larger amplitudes the response becomes quasiperiodic, as shown in Fig. 12(c).

Case two
Here the non-dimensional parameters in Eq. (12) are set as follows: stiffness ratios a ¼ 1:5 and b ¼ 3:0; and mass ratio g ¼ 1:0: Compared to the first case, the stiffness of the clearance spring k 3 has been doubled. Therefore, the non-linearity is larger than in the first case, and should result in more pronounced distortions of the NNM invariant manifolds. Since the linear modal parameters in Eq. (13) are independent of the parameter b; they are the same as in the first case.
The two NNM manifolds are constructed numerically following the same procedures as the first case: (1) for the first NNM, the maximum amplitude is a max ¼ 15; which is evenly divided into 30 segments, and N f ¼ 64 harmonic terms are used; (2) for the second NNM, the maximum amplitude is a max ¼ 60 with 120 segments, and N f ¼ 32 harmonic terms. The switching hyperplane is located at the amplitude aE2:11 for both modes. The NNM manifolds are shown in Fig. 13 in the co-ordinate system ðz i ; ' z i Þ: For the second NNM, the manifold looks similar to that for the first case. However, the kink in the first NNM is so pronounced in this case that the manifold is no longer single valued beyond a certain amplitude. This results from the stronger non-linearity.
The time history of the motion on the first NNM is shown in Fig. 14. Time integrations of the system equations of motion and of the dynamics restricted to the Galerkin-based manifold are compared. In all cases the periodic motions match very precisely. In the kink region of the manifold, the motion of mass m 1 has a slight oscillation superposed to the overall motion, instead of the near-stick phenomenon-observed in the first case. This results in the multiple valued region, or loop, shown in the invariant manifold in Fig. 13.
The accuracy of the NNMs is further verified in the frequency domain, as shown in Fig. 15. The frequency-amplitude relationship is similar to that in the first case. Comparison of the  Galerkin-based manifold simulation results and solutions from simulations of the system equations of motion demonstrate the excellent accuracy of the numerically computed NNMs. The stability and large amplitude dynamics of these NNMs were also checked using characteristic multipliers and Poincar! e maps. Similar results for the NNM bifurcations are found, as follows: (1) for the first NNM, the manifold is stable until a critical amplitude a c =3.82, where a pair of complex conjugate characteristic multipliers merge at À1; solutions on the manifold remain unstable as the amplitude increases up to another critical amplitude, a c ¼ 5:62; where the multipliers merge again at À1; and solutions on the manifold regain stability. (2) For the second NNM, solutions on the manifold are unstable above the critical amplitude of a c ¼ 2:30: Sample Poincar! e sections are shown in Fig. 16 to illustrate the bifurcation of the first NNM.

Conclusions
From this study of NNMs for piecewise linear autonomous systems, and the example system studied in detail, the following conclusions are drawn. (1) The Galerkin-based method, originally developed for dynamic systems with smooth non-linearities, can be extended to piecewise linear systems and used to accurately construct NNM invariant manifolds. The transformation of the master co-ordinates to polar form and the discretization in the amplitude make the Galerkinbased approach applicable in strongly non-linear regions, as well as in the transition region. (2) The dynamic response on individual NNM manifolds can be reduced to a single-d.o.f. system described in terms of the master co-ordinates. The dynamic response of all slave co-ordinates can be recovered from the simulation results of the master co-ordinates using the constraint relations.
(3) Although numerical results were obtained for a two-d.o.f. system, the Galerkin-based approach can be easily applied to multi-d.o.f. piecewise linear systems, so long as the single switching hyperplane condition is satisfied. (4) For response amplitudes beyond the transition hyperplane, the dynamic behavior of piecewise linear systems can be quite complicated. This includes non-trivial aspects of the periodic response, such as loops in the manifolds, as well as instabilities leading to a variety of system responses. (5) The stability and post-critical dynamics of the non-linear normal modes were investigated using characteristic multipliers and Poincar! e maps. Flip bifurcations were found to occur for both modes, as well as transitions to quasiperiodic responses. For the first non-linear mode, the NNM motions regained stability beyond a second bifurcation amplitude.