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

Dongying Jiang, Christophe Pierre, Steven Shaw

► To cite this version:

Dongying Jiang, Christophe Pierre, Steven Shaw. Large-amplitude non-linear normal modes of piecewise linear systems. Journal of Sound and Vibration, 2004, 272, pp.869-891. 10.1016/S0022-460X(03)00497-8. hal-01350797

HAL Id: hal-01350797

https://hal.science/hal-01350797

Submitted on 1 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Large-amplitude non-linear normal modes of piecewise linear systems

D. Jiangª, C. Pierreª, S.W. Shawᵇ

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 this construction 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 NNM for this example are investigated.

## 1. Introduction

The concept of non-linear normal modes was originated by Rosenberg [1–3], based on the analysis of discrete, symmetric systems with smooth non-linearities. Since then, the existence [4,5], stability [6–10], and construction [11–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–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

<footnote>ªDepartment of Mechanical Engineering, University of Michigan, Ann Arbor, MI 48109-2125, USA</footnote>

<footnote>ᵇDepartment of Mechanical Engineering, Michigan State University, East Lansing, MI 48824-1226, USA</footnote>
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 non-linearities. 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–21], including dynamic behaviors such as bifurcations and chaos [22–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 non-gyroscopic 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 non-linearity, 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.
## 2. 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

$$ \begin{cases} M\ddot{Z} + K_1Z = 0 & \text{for } h^T Z \le d, \\ M\ddot{Z} + K_2Z = b & \text{for } h^T Z > d, \end{cases} \qquad (1) $$

where $d > 0$ is a real scalar constant, $Z, h, b \in \mathbb{R}^N$, and $M, K_1, K_2$ are real symmetric positive definite $N \times N$ matrices. Obviously, this system has two distinct linear regions separated by a hyperplane, $\{Z \in \mathbb{R}^N: h^T Z = d\}$. We denote the first region as $\{Z \in \mathbb{R}^N: h^T Z < d\}$ and the second region as $\{Z \in \mathbb{R}^N: h^T Z > d\}$. 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, $\{Z \in \mathbb{R}^N: h^T Z < d\}$, 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 = Q\eta$, system (1) takes the standard form

$$ \ddot{\eta} + \Lambda_1 \eta = f(\eta), \qquad (2) $$

where $\eta \in \mathbb{R}^N$ is the vector of modal co-ordinates, and the $N \times N$ diagonal matrix $\Lambda_1$ has entries that are the squares of the small-amplitude (first region) natural frequencies. The piecewise linear force vector, $f(\eta) \in \mathbb{R}^N$, is given by

$$ f(\eta) = H(h^T Q \eta - d)[(\Lambda_1 - Q^T K_2 Q)\eta + Q^T b], \qquad (3) $$

where $H(.)$ is the heaviside function that stands for the switching hyperplane.

## 3. 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 co-ordinates, 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 displacement-velocity 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 $(\eta_k, \dot{\eta}_k)$. These are transformed to amplitude and phase co-ordinates, $(a, \phi)$, as follows:

$$ 
\begin{aligned}
\eta_k(t) &= a(t) \cos \phi(t), \\
\dot{\eta}_k(t) &= -\omega_k a(t) \sin \phi(t),
\end{aligned}
\qquad (4) $$

where $\omega_k$ is the kth natural frequency of the linear system in the first region. The master co-ordinates $(a, \phi)$ have bounded domains for amplitude and phase, $a \in [0, a_{max}]$ and $\phi \in [0, 2\pi]$, 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:

$$ 
\begin{aligned}
\eta_i &= P_i(a, \phi), \\
\dot{\eta}_i &= Q_i(a, \phi),
\end{aligned}
\qquad (5) $$

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(\eta)$ 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 co-ordinate $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]

$$ 
\begin{aligned}
Q_i &= \frac{\partial P_i}{\partial a} \left( \frac{-f_k}{\omega_k} \right) \sin \phi + \frac{\partial P_i}{\partial \phi} \left( \omega_k - \frac{f_k \cos \phi}{a \omega_k} \right) \\
&\quad - \omega_i^2 P_i + f_i = \frac{\partial Q_i}{\partial a} \left( \frac{-f_k}{\omega_k} \right) \sin \phi + \frac{\partial Q_i}{\partial \phi} \left( \omega_k - \frac{f_k \cos \phi}{a \omega_k} \right)
\end{aligned}
\qquad (6) $$

These are obviously valid when the force, $f \in \mathbb{R}^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, $\partial P_i/\partial a$, $\partial P_i/\partial \phi$, $\partial Q_i/\partial a$, and $\partial Q_i/\partial \phi$, 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, $a \in [0, a_{max}]$ and $\phi \in [0, 2\pi]$. In the $\phi$ 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 $a \in [0, a_{max}]$ is divided into $n$ equal segments, defined by

$$a \in [a_j, a_{j+1}], \quad a_j = \frac{ja_{max}}{n}, \quad j = 0, 1, \dots, 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

$$P_i(a, \phi) = \sum_{l=1}^{2} \sum_{m=1}^{N_{\phi}} C_i^{l,m} T_{l,m}(a, \phi) \quad \text{for } i = 1, 2, \dots, N, \quad i \neq k, \quad (7)$$

$$Q_i(a, \phi) = \sum_{l=1}^{2} \sum_{m=1}^{N_{\phi}} D_i^{l,m} U_{l,m}(a, \phi)$$

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 $\phi$. Hence, for a given segment,

$$T_{l,m} = A_l(a) \cos[(m-1)\phi], \quad U_{l,m} = A_l(a) \sin(m\phi),$$

where

$$A_1(a) = \frac{a - a_j}{a_{j+1} - a_j} \quad \text{and} \quad A_2(a) = \frac{a_{j+1} - a}{a_{j+1} - a_j}$$

are the tent functions employed in the segment $a \in [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

$$0 = \int_{a, \phi} T_{p,q} \left[ -a \sum_{l,m} D_i^{l,m} U_{l,m} + \sum_{l,m} C_i^{l,m} \frac{\partial T_{l,m}}{\partial a} \frac{-f_k}{\omega_k} a \sin \phi \right. \\ \left. + \sum_{l,m} C_i^{l,m} \frac{\partial T_{l,m}}{\partial \phi} \left( a\omega_k - \frac{f_k \cos \phi}{\omega_k} \right) \right] da \, d\phi, \quad (8)$$

$$0 = \int_{a, \phi} U_{p,q} \left[ \omega_i^2 a \sum_{l,m} C_i^{l,m} T_{l,m} - af_i + \sum_{l,m} D_i^{l,m} \frac{\partial U_{l,m}}{\partial a} \frac{-f_k}{\omega_k} a \sin \phi \right. \\ \left. + \sum_{l,m} D_i^{l,m} \frac{\partial U_{l,m}}{\partial \phi} \left( a\omega_k - \frac{f_k \cos \phi}{\omega_k} \right) \right] da \, d\phi \quad (9)$$
for $i = 1...N$, $i \neq k$, $p = 1,2$, and $q = 1...N_{\phi}$. This results in $2(N-1)2N_{\phi}$ 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 $\phi$ [18], as follows:

$$ \dot{a} = -\frac{f_k}{\omega_k} \sin \phi, \quad \dot{\phi} = \omega_k - \frac{f_k}{a\omega_k} \cos \phi, \qquad (10) $$

where $f_k$ depends on $a$ and $\phi$ 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.

## 4. A two-d.o.f. case study

An autonomous two-d.o.f. piecewise linear system is studied here to illustrate the construction procedures described in Sections 2 and 3. The system shown in Fig. 1 is composed of two masses linked with linear springs, which is similar to the example system considered in Refs. [26,27].

Let $x_1$ and $x_2$ represent the displacements of the first and second masses, $m_1$ and $m_2$, respectively, from their static equilibrium positions. The equations of motion are
<figure><img src="image_2.png" /><figcaption>Fig. 1. A two-degree-of-freedom piecewise linear system.</figcaption></figure>

given by

$$
\begin{align}
m_1 \ddot{x}_1 + k_1 x_1 - k_2 (x_2 - x_1) &= 0 \nonumber \\
m_2 \ddot{x}_2 + k_2 (x_2 - x_1) &= 0 \nonumber \\
m_1 \ddot{x}_1 + k_1 x_1 - k_2 (x_2 - x_1) &= 0 \nonumber \\
m_2 \ddot{x}_2 + k_2 (x_2 - x_1) + k_3 (x_2 - d_0) &= 0 \tag{11}
\end{align}
$$

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, $z_1 = x_1/d_0$, $z_2 = x_2/d_0$, $\tau = t/\sqrt{m_1/k_1}$, $\alpha = k_2/k_1$, $\beta = k_3/k_1$, and $\gamma = m_2/m_1$,
yielding

$$
\begin{bmatrix} 1 & 0 \\ 0 & \gamma \end{bmatrix} \begin{Bmatrix} \ddot{z}_1 \\ \ddot{z}_2 \end{Bmatrix} + \begin{bmatrix} 1+\alpha & -\alpha \\ -\alpha & \alpha \end{bmatrix} \begin{Bmatrix} z_1 \\ z_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix} \quad \text{if } z_2 < 1,
$$
$$
\begin{bmatrix} 1 & 0 \\ 0 & \gamma \end{bmatrix} \begin{Bmatrix} \ddot{z}_1 \\ \ddot{z}_2 \end{Bmatrix} + \begin{bmatrix} 1+\alpha & -\alpha \\ -\alpha & \alpha+\beta \end{bmatrix} \begin{Bmatrix} z_1 \\ z_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ \beta \end{Bmatrix} \quad \text{if } z_2 > 1, \tag{12}
$$

where the double derivative (\ddot{)} now denotes $d^2() / d\tau^2$.

Compared to the form of Eq. (1) in Section 2, the mass matrix *M*, and the stiffness matrices *K*₁
and *K*₂ are

$$
M = \begin{bmatrix} 1 & 0 \\ 0 & \gamma \end{bmatrix}, \quad K_1 = \begin{bmatrix} 1+\alpha & -\alpha \\ -\alpha & \alpha \end{bmatrix}, \quad K_2 = \begin{bmatrix} 1+\alpha & -\alpha \\ -\alpha & \alpha+\beta \end{bmatrix}.
$$

The displacement vector is $Z = [z_1 \ z_2]^T$, and the constant offset vector is $b = [0 \ \beta]^T$. The switching plane is defined by $\{Z: h^T Z = d\}$, where $h = [0 \ 1]^T$ and $d = 1$.

Linear modal analysis can be performed on the sub-system $M\ddot{Z} + K_1Z = 0$ in Eq. (12), where $Z = 0$ is the static equilibrium position. For the generalized eigenvalue problem $K_1q = \omega^2 Mq$, the two real positive eigenvalues are found to be

$$
\omega_{1,2}^2 = \frac{1}{2\gamma}[(\alpha + \alpha\gamma + \gamma) \mp \sqrt{(\alpha + \alpha\gamma + \gamma)^2 - 4\alpha\gamma}] \quad (13)
$$
and the eigenvector matrix is given by $\tilde{Q} = [q_1 \ q_2]$, where $q_1$ and $q_2$ are the eigenvectors corresponding to the non-dimensional natural frequencies $\omega_1$ and $\omega_2$. After applying the modal transformation $Z = \tilde{Q}\eta$ and some manipulations, Eq. (12) is expressed in the standard form of Eq. (2), and is thus ready for the construction of its NNMs.

## 4.1. Case one

The non-dimensional stiffness ratios are taken to be $\alpha = \beta = 1.5$ and the mass ratio is $\gamma = 1.0$. The non-dimensional linear natural frequencies in Eq. (13) are $\omega_1 = 0.6472$ and $\omega_2 = 1.8924$, and the corresponding eigenvector matrix is given by

$$ \tilde{Q} = \begin{bmatrix} 0.3419 & 0.6581 \\ 0.4743 & -0.4743 \end{bmatrix}, $$

which has not been normalized with respect to the mass matrix.

### 4.1.1. 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, $\eta_1$ and $\dot{\eta}_1$ are chosen as master co-ordinates, where $\eta_1 = a \cos \phi$ and $\dot{\eta}_1 = -a\omega_1 \sin \phi$. The constraint equations $\eta_2 = P_2(a, \phi)$ and $\dot{\eta}_2 = Q_2(a, \phi)$ are numerically constructed over the domain $a \in [0, a_{max}]$, $\phi \in [0, 2\pi]$. 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_{\phi} = 64$. For the second NNM, the master co-ordinates are $(\eta_2, \dot{\eta}_2)$, and the constraint equations are $\eta_1 = P_1(a, \phi)$ and $\dot{\eta}_1 = Q_1(a, \phi)$. The parameters set in the numerical solution algorithm are: $a_{max} = 60$, with 120 segments in the $a$ direction, and $N_{\phi} = 32$.

The position constraints $P_2(a, \phi)$ and $P_1(a, \phi)$ 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 $a \approx 2.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 non-linearity, where asymptotic analyses [26] are not applicable. For the second NNM, the switching amplitude is also at $a \approx 2.11$, and the maximum displacement $z_2$ is about 28 when $a_{max}$ is reached.

It should be noted here that the switching position, $a \approx 2.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é 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 $a \in [2.1082, 2.5]$ is
<figure><img src="image_3.png" /><figcaption>Fig. 2. The ith non-linear mode invariant manifold for the position constraint, $\eta_i = P_i(a, \phi)$: (a) first mode, (b) second mode, as a function of the amplitude a and phase $\phi$.</figcaption></figure>

evenly divided into 100 segments with $N_{\phi} = 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, \dot{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
<figure><img src="image_4.png" /><figcaption>Fig. 3. The refined first non-linear mode invariant manifold for the position constraint, $\eta_2 = P_2(a, \phi)$, as a function of the amplitude $a$ and phase $\phi$.</figcaption></figure>

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 $\dot{z}_1$ is nearly zero, the displacement $z_1$ is almost constant, the velocity $\dot{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 $\phi$, which describe the dynamics on the manifold, Eq. (10). Initial conditions $a(0) > 0$ and $\phi(0) = 0$ are used, and since the motion on the manifold is periodic, its amplitude is equal to $a(0)$. Then, the modal responses $\eta_i(t)$ and $\dot{\eta}_i(t)$ are obtained from $a(t)$ and $\phi(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 $\dot{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
<figure><img src="image_5.png" /><figcaption>Fig. 4. The non-linear mode invariant manifold for the position constraint, $z_2(z_1, \dot{z}_1)$: (a) first mode, (b) second mode, as a function of the non-dimensional displacement and velocity ($z_1, \dot{z}_1$). A representative motion on the manifold is shown (—).</figcaption></figure>

fundamental frequency associated with the basic period $\tau_0$, $f_0 = 2\pi/\tau_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
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Time histories of the motion on the invariant manifold: —, solution from direct integration of the system equations of motion; ---, time simulation of the dynamics restricted to the Galerkin-based manifold. (a) First NNM, amplitude $a = 24$; (b) second NNM, amplitude $a = 35$.</figcaption></figure>

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
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Response frequency for the non-linear normal mode: —, solution from direct integration of the system equations of motion; ○, time simulation of the dynamics restricted to the Galerkin-based manifold. (a) first mode; (b) second mode.</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Phase portraits of the periodic motion on the invariant manifold: —, mass one; ---, mass two. (a) First mode, amplitude $a = 12, 24$; (b) second mode, amplitude $a = 17.5, 35$.</figcaption></figure>

space, shown in Fig. 8, the trajectories of solutions are represented by curves, but not by straight lines as in the linear case.<fnref n="1" /> Moreover, the displacements $z_1$ and $z_2$ do not vanish simultaneously,<fnref n="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

<footnote n="1">¹These are not tight loops, as can be shown from the symmetry of the manifold geometry, Fig. 4, or the Poincaré section in Figs. 11 and 12.</footnote>

<footnote n="2">²This is apparent in Fig. 8 for the second NNM, and is also true for the first NNM, upon close examination.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 8. Periodic trajectories in the $(z_1, z_2)$ configuration space: —, first NNM, amplitude $a = 12, 24$; ---, second NNM, amplitude $a = 17.5, 35$.</figcaption></figure>

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].

### 4.1.2. 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é 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,

$$ \dot{Y} = F(Y), \quad (14) $$

where $Y = [z_1 \dot{z}_1 z_2 \dot{z}_2]^T$ and the right-hand side is given by

$$ F(Y) = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -(1+\alpha) & \alpha & 0 & 0 \\ \alpha/\gamma & -\alpha/\gamma - \beta H(z_2-1)/\gamma & 0 & 0 \end{bmatrix} Y + \begin{Bmatrix} 0 \\ 0 \\ 0 \\ \beta/\gamma \end{Bmatrix} H(z_2-1). $$

The characteristic multipliers can be determined from the monodromy matrix $M$, defined by $M(a) = \Phi(\tau_0)$, where $\tau_0$ is the period of the motion, $a$ is the amplitude of the periodic motion, and the matrix $\Phi(\tau_0)$ is determined from the matrix initial value problem,

$$ \dot{\Phi}(\tau) = F_Y(Y_p)\Phi, \quad \Phi(0) = I, \quad (15) $$
<figure><img src="image_13.png" /><figcaption>Fig. 9. Eigenvalues of the monodromy matrix in the complex plane as the NNM amplitude *a* is increased: (a) first mode; (b) second mode.</figcaption></figure>

where $Y_p$ is the periodic solution with amplitude *a*, and the period $\tau_0$ is obtained while numerically searching for periodic solutions.<fnref n="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, $\mu(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, $\mu(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

<footnote n="3">³Note that the period $\tau_0$ depends on the amplitude *a*.</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 10. FFT of the time response $z_1$ on the first NNM: (a) amplitude $a = 6.52$; (b) $a = 6.53$.</figcaption></figure>

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 $a \ge 8.62$.

The results from the characteristic multipliers indicate that the system loses stability when the multipliers satisfy $\mu(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 $\beta$ in Eq. (12).

In order to explore the dynamics in the unstable region, a Poincaré map is now employed [31]. Since there is only a single transition hyperplane in the system, it provides a natural Poincaré section [26], defined by the two components<fnref n="4" />

$$ 
\begin{aligned}
\Sigma &= \{ Y \in \mathbb{R}^4: z_2 = 1, \text{ and } \dot{z}_2 > 0 \}, \\
\Sigma^* &= \{ Y \in \mathbb{R}^4: z_2 = 1, \text{ and } \dot{z}_2 < 0 \},
\end{aligned} 
\quad (16) 
$$

<footnote n="4">⁴The Poincaré section should be defined such that periodic motions pierce it only once during one period. In the present case, either section $\Sigma$ or section $\Sigma^*$ is sufficient for this purpose, but solutions for both sections are shown to evidence the symmetry of the periodic responses.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Poincaré section of the motion initiated on the first NNM versus amplitude *a*: (a) *a* = 6.52; (b) *a* = 6.53; (c) *a* = 7.5; (d) *a* = 8.61; (e) *a* = 8.62.</figcaption></figure>

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é mapping can thus be constructed numerically [26].

From a cross-sectional view of the Poincaré map, the motion on the first NNM is stable until the critical amplitude, *a*<sub>c</sub> = 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*<sub>c</sub> = 4.04, as shown in Fig. 12(b). Here the bifurcation appears to be supercritical, and another period doubling occurs,
<figure><img src="image_16.png" /><figcaption>Fig. 12. Poincaré section of the motion initiated on the second NNM versus amplitude *a*: (a) *a* = 4.04; (b) *a* = 4.05; (c) *a* = 6.0.</figcaption></figure>

resulting in the transition from two to four points. At larger amplitudes the response becomes
quasiperiodic, as shown in Fig. 12(c).

4.2. Case two

Here the non-dimensional parameters in Eq. (12) are set as follows: stiffness ratios $\alpha = 1.5$ and $\beta = 3.0$, and mass ratio $\gamma = 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 $\beta$, 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_{\phi} = 64$ harmonic terms are used; (2) for the second NNM, the maximum amplitude is $a_{max} = 60$ with 120 segments, and $N_{\phi} = 32$ harmonic terms. The switching hyperplane is located at the amplitude $a \approx 2.11$ for both modes.
<figure><img src="image_17.png" /><figcaption>Fig. 13. The NNM invariant manifold for the position constraint, $z_2(z_1, \dot{z}_1)$: (a) first mode, (b) second mode, as a function of the non-dimensional displacement and velocity $(z_1, \dot{z}_1)$. A representative motion on the manifold is shown (—).</figcaption></figure>

The NNM manifolds are shown in Fig. 13 in the co-ordinate system $(z_i, \dot{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
<figure><img src="image_18.png" /><figcaption>Fig. 14. Time histories of the motion on the first NNM invariant manifold at amplitude $a = 8.5 : -, $ solution from direct integration of the system equations of motion; $--$, time simulation of the dynamics restricted to the Galerkin-based manifold.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15. Response frequency for the non-linear normal mode: $-$, solution from direct integration of the system equations of motion; $O$, time simulation of the dynamics restricted to the Galerkin-based manifold. (a) First mode; (b) second mode.</figcaption></figure>

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
<figure><img src="image_20.png" /><figcaption>Fig. 16. Poincaré sections of the motion initiated on the first NNM versus amplitude *a*: (a) *a* = 3.82; (b) *a* = 3.83.</figcaption></figure>

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é 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*<sub>c</sub> = 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*<sub>c</sub> = 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*<sub>c</sub> = 2.30. Sample Poincaré sections are shown in Fig. 16 to illustrate the bifurcation of the first NNM.

**5. 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 Galerkin-based 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é 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.

## Acknowledgements

This work was supported by a grant from the Army Research Office. Dr. Gary Anderson is the grant monitor. The authors are grateful to Vincent Soumier for his help with the computations.

## References

[1] R.M. Rosenberg, The normal modes of non-linear n-degree-of-freedom systems, *American Society of Mechanical Engineers, Journal of Applied Mechanics* 30 (1962) 7–14.

[2] R.M. Rosenberg, J.K. Kuo, Nonsimilar normal mode vibrations of non-linear systems having two degrees of freedom, *American Society of Mechanical Engineers, Journal of Applied Mechanics* 31 (1964) 283–290.

[3] R.M. Rosenberg, On non-linear vibrations of systems with many degrees of freedom, *Advances in Applied Mechanics* 9 (1966) 155–242.

[4] C. Cooke, R. Struble, The existence of periodic solutions and normal mode vibrations in non-linear systems, *Quarterly of Applied Mathematics* 24 (1966) 177–193.

[5] C. Pak, R. Rosenberg, On the existence of normal mode vibrations in non-linear systems, *Quarterly of Applied Mathematics* 26 (1968) 403–416.

[6] C. Atkinson, B. Taskett, A study of the non-linearly related modal solutions of coupled non-linear systems by superposition techniques, *American Society of Mechanical Engineers, Journal of Applied Mechanics* 32 (1965) 359–364.

[7] R. Rand, The geometrical stability of non-linear normal modes in 2 degrees of freedom systems, *International Journal of Non-Linear Mechanics* 8 (1973) 161–168.

[8] L.A. Month, R.H. Rand, An application of poincaré map to the stability of non-linear normal modes, *American Society of Mechanical Engineers, Journal of Applied Mechanics* 47 (1980) 645–651.

[9] T.K. Caughey, A.F. Vakakis, A method of examining steady state solutions of forced discrete systems with strong non-linearities, *International Journal of Non-Linear Mechanics* 26 (1991) 89–103.

[10] C. Pak, R. Rand, A. Vakakis, Bifurcation of non-linear normal modes in a class of 2 degrees of freedom systems, *Acta Mechanica* 3 (1992) 129–145.

[11] H. Greenberg, T.L. Yang, Modal subspaces and normal mode vibrations, *International Journal of Non-Linear Mechanics* 6 (1971) 311–326.

[12] R. Rand, Nonlinear normal mode in 2 degrees of freedom systems, *American Society of Mechanical Engineers, Journal of Applied Mechanics* 38 (1971) 561.

[13] R. Rand, A direct method for non-linear normal modes, *International Journal of Non-Linear Mechanics* 9 (1974) 363–368.

[14] J.C. Slater, A numerical method for determining non-linear normal modes, *Nonlinear Dynamics* 10 (1996) 19–30.

[15] S.W. Shaw, C. Pierre, Non-linear normal modes and invariant manifolds, *Journal of Sound and Vibration* 150 (1991) 170–173.

[16] S.W. Shaw, C. Pierre, Normal modes for non-linear vibratory systems, *Journal of Sound and Vibration* 164 (1993) 85–124.
[17] S.W. Shaw, C. Pierre, Normal modes of vibration for non-linear continuous systems, *Journal of Sound and Vibration* 169 (1994) 319-347.

[18] E. Pescheck, C. Pierre, S.W. Shaw, A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds, *Journal of Sound and Vibration* 249 (2002) 971-993.

[19] Y.S. Choi, S.T. Noah, Forced periodic vibration of unsymmetric piecewise-linear systems, *Journal of Sound and Vibration* 121 (1988) 117-126.

[20] R.J. Comparin, R. Singh, Non-linear frequency response characteristics of an impact pair, *Journal of Sound and Vibration* 134 (1989) 259-290.

[21] R.J. Comparin, R. Singh, Frequency response characteristics of a multi-degree-of-freedom system with clearances, *Journal of Sound and Vibration* 142 (1990) 101-124.

[22] S.W. Shaw, P.J. Holmes, A periodically forced piecewise linear oscillator, *Journal of Sound and Vibration* 90 (1983) 129-155.

[23] J.M.T. Thompson, H.B. Stewart, *Nonlinear Dynamics and Chaos*, Wiley, Chichester, 1986.

[24] I.A. Mahfouz, F. Badrakhan, Chaotic behavior of some piecewise linear systems, part i: systems with set-up spring or with unsymmetric elasticity, *Journal of Sound and Vibration* 143 (1990) 255-288.

[25] L. Zuo, A. Curnier, Non-linear real and complex modes of conewise linear systems, *Journal of Sound and Vibration* 174 (1994) 289-313.

[26] S. Chen, S.W. Shaw, Normal modes for piecewise linear vibratory systems, *Nonlinear Dynamics* 10 (1996) 135-163.

[27] M. Chati, R. Rand, S. Mukherjee, Modal analysis of a cracked beam, *Journal of Sound and Vibration* 207 (1997) 249-270.

[28] A.P. Ivanov, The stability of periodic solutions of discontinuity systems that intersect several surfaces of discontinuity, *Journal of Applied Mathematics and Mechanics* 62 (1998) 677-685.

[29] R. Seydel, *From Equilibrium to Chaos—Practical Bifurcation and Stability Analysis*, Elsevier, New York, 1988.

[30] M.A. Aizerman, F.R. Gantmakher, Stability in the linear approximation of a periodic solution of a system of differential equations with discontinuous right-hand sides, *Prikladnaya Matematika I Mekhanika* 21 (1957) 658-669.

[31] J. Guckenheimer, P. Holmes, *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*, Applied Mathematical Sciences, Vol. 42, Springer, New York, 1989.