# Normal modes for piecewise linear vibratory systems

Shyh-Leh Chen, Steven Shaw

► To cite this version:

Shyh-Leh Chen, Steven Shaw. Normal modes for piecewise linear vibratory systems. Nonlinear Dynamics, 1996, 10 (2), pp.135-164. 10.1007/BF00045454. hal-01354038

HAL Id: hal-01354038

https://hal.science/hal-01354038

Submitted on 17 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
# Normal Modes for Piecewise Linear Vibratory Systems

SHYH-LEH CHEN and STEVEN W. SHAW

Department of Mechanical Engineering, Michigan State University, East Lansing, MI 48824, U.S.A.

**Abstract.** A method to construct the normal modes for a class of piecewise linear vibratory systems is developed in this study. The approach utilizes the concepts of Poincaré maps and invariant manifolds from the theory of dynamical systems. In contrast to conventional methods for smooth systems, which expand normal modes in a series form around an equilibrium point of interest, the present method expands the normal modes in a series form of polar coordinates in a neighborhood of an invariant disk of the system. It is found that the normal modes, modal dynamics and frequency-amplitude dependence relationship are all of piecewise type. A two degree of freedom example is used to demonstrate the method.

**Key words:** Nonlinear normal modes, piecewise linear systems, invariant manifold, asymptotic expansion.

## 1. Introduction

In linear vibration theory, modal analysis is a powerful technique that provides a means of reducing the system dynamics to uncoupled individual oscillators. Through the mode shapes and corresponding modal dynamics, the behavior of a linear vibratory system can be completely understood. For this reason, researchers have tried to extend the concept of normal modes to nonlinear systems. Rosenberg took the first steps in this direction [22–24]. Since then, many studies have appeared in the literature discussing the existence [5, 15, 19, 40], stability [1, 7, 10, 17, 21], and construction [8, 18–20, 39] of nonlinear normal modes, and the possibility of some type of nonlinear modal analysis [1, 32, 35–38]. Other studies have included the investigation of nonlinear normal modes for continuous [2, 3, 12, 13, 16], nonconservative [6, 28], and internally resonant systems [14]. Shaw and Pierre introduced a constructive method for nonlinear normal modes using invariant manifolds [26, 28–30]. This approach is applied here to the case of piecewise linear systems.

While the idea of nonlinear normal modes has drawn much attention, virtually all of the previous works deal with smooth systems. In contrast, nonsmooth systems, in particular piecewise linear (PWL) systems (see, e.g., [4, 9, 11, 25, 31]), received little attention in regards to normal modes. However, PWL vibratory systems are useful models for many practical vibratory systems. For example, conservative systems with clearance, backlash, or piecewise linear springs and dissipative systems with simple coulomb friction are all of PWL type.

Recently, Zuo and Curnier reported a study on the modal motions of a class of PWL systems, where the switching hyperplane passes through the origin [41]. As one will see in the sequel, their class is a special case of the present class of PWL systems to be studied, where the switching hyperplane does not necessarily pass through the origin. In particular, the present systems may not satisfy the property of *positive homogeneity* defined in [41]. As a result, the modal frequency of the present systems depends on the amplitude, in contrast with the systems in [41] whose modal frequency is independent of the amplitude. Also, the
construction of the normal mode manifolds and their attendant dynamics were not pursued in [41]. However, they did include systems with gyroscopic terms, which are not considered here.

The objective of this paper is to demonstrate a method for the construction of nonlinear normal modes and their attendant dynamics for a class of general conservative PWL systems with N (finite) degrees of freedom. To this end, we utilize the approach of invariant manifolds and asymptotic expansions. Here the invariant manifold, and hence the normal mode, is a one-parameter family of periodic orbits in the state space. The general procedure involves building a Poincaré map with the switching hyperplane as the Poincaré section in state space, determining the fixed points of the map, and finally constructing the normal mode by asymptotic expansions. Unlike the case for a smooth system, where the normal mode is expanded in a series form near an equilibrium point of the system, the normal mode for a PWL system is expanded in a series form using polar coordinates in a neighborhood of an invariant disk of the system. This invariant disk is precisely the corresponding normal mode manifold of the linear system for energies below those at which the switching takes place.

The paper is organized as follows. We begin in Section 2 with a description of the class of PWL systems in which we are interested. The general system will then be transformed into a canonical form for easy formulation. Next, in Section 3, the Poincaré map is obtained implicitly. The procedure of constructing the normal mode manifolds for PWL systems, including the determination of the fixed points of the Poincaré map, using asymptotic expansion is described in Section 4. Their modal dynamics are established in Section 5. In Section 6, a two degrees of freedom example problem with clearance is given to show some of the calculations involved with using this method and a comparison is made with simulation. Some conclusions are drawn in Section 7.

## 2. The Piecewise Linear System and Its Canonical Form

Consider below an unforced, undamped N d.o.f. system

$$ \begin{cases} M\ddot{z} + K_1z = 0 & \text{for } h^Tz \le d' \\ M\ddot{z} + K_2z = b & \text{for } h^Tz > d' \end{cases} \quad (1) $$

where $d' > 0$ is a scalar constant, $z, h, b \in \mathbb{R}^N$, and $M, K_1, K_2$ are real symmetric positive definite $N \times N$ matrices. This is a PWL system with two linear domains. $\{z \in \mathbb{R}^N : h^Tz = d'\}$ is a hyperplane dividing $\mathbb{R}^N$ into two regions: the first region $\{z \in \mathbb{R}^N : h^Tz \le d'\}$ and the second one $\{z \in \mathbb{R}^N : h^Tz > d'\}$. Each region is governed by an N d.o.f. linear vibratory subsystem. Since $M, K_1$, and $K_2$ are real symmetric and positive definite, both linear subsystems have N natural modes. For small amplitude motions, when the motions never leave the first region, the system is nothing but the first linear subsystem and is well understood. As the motion amplitudes become large, solutions pass through the switching hyperplane and go into the second region, in which case both linear subsystems will be involved and the system behavior is no longer simple. Many complicated behaviors involved in PWL systems have been reported [4, 31, 34]. In this paper, we are not going to discuss the general behavior of such PWL systems, but only their normal mode motions.

System (1) represents a large class of physical vibratory systems such as systems with clearance, impact, or piecewise linear springs. However, we must point out that it can only represent systems with one clearance, but not those with two or more clearances, which should
be modeled by PWL systems with more domains. Normal modes for PWL systems with more
than two domains are not considered in the present paper. Also, systems with dry friction do
not belong to the class of systems in equation (1) because they are nonconservative oscillators.

Before proceeding to the next step, we make some assumptions about system (1). By premultiplying equation (1) with $M^{-1}$, the system can be rewritten as

$$
\ddot{z} = \begin{cases} A_1 z & \text{for } h^T z \le d' \\ A_2 z + b' & \text{for } h^T z > d' \end{cases}, \quad (2)
$$

where $A_1 = -M^{-1}K_1$, $A_2 = -M^{-1}K_2$, and $b' = M^{-1}b$. The assumptions are

(i) $A_1 z = A_2 z + b'$ when $h^T z = d'$;

(ii) N modes in both linear subsystems are distinct.

Assumption (i) says that the right-hand side of (2) is continuous, so that when written in state equation form, i.e., as a system of first order ODEs, the system will possess a continuous vector field. Although the vector field is not $C^1$ due to its PWL nature, it will satisfy a local Lipschitz condition. Hence, the existence and uniqueness of solutions of (2) are guaranteed. Assumption (ii) is for simplicity and can be removed easily. It is needed to insure that each mode of the linear subsystems is two dimensional in the state space. Without this assumption, we might have a four or more dimensional surface for a certain mode. Also, it can be shown that assumption (i) together with the negative definiteness of $K$ and $\hat{K}$ imply that the origin is the only equilibrium point for system (6) (this can be proved but is not obvious).

With these assumptions in hand, we transform system (2) into a simpler form. Let $x = Q^T z$,
where $Q = [q_1 \ q_2 \ \dots \ q_N]$ is an $N \times N$ matrix with $q_1 = \mu h$ and $\mu \in \mathbb{R}$ is a suitable constant
such that $QQ^T = M$. In other words, the switching hyperplane becomes $\{x \in \mathbb{R}^N : x_1 = d\}$,
where $d = \mu d'$. Such a transformation matrix $Q$ does exist and can be obtained in the following
way. Recalling that $M$ is real, symmetric and positive definite, we write $M = PDP^T$, where
$P$ is orthogonal and $D$ is diagonal. Take another orthogonal matrix $U$ with the property that
$u_1 = \mu D^{-1/2} P^T h$, where $u_1$ is the first column of $U$. Hence $\mu = 1/\|D^{-1/2} P^T h\|$ where $\|\cdot\|$ denotes the Euclidean 2-norm. Then it is easy to verify that $Q = PD^{1/2}U$ is as required.

After the transformation $x = Q^T z$, system (2) takes the form

$$
\ddot{x} = \begin{cases} Kx & \text{for } x_1 \le d \\ \hat{K}x + Q^T b' & \text{for } x_1 > d \end{cases}, \quad (3)
$$

where $K = Q^T A_1 (Q^T)^{-1}$ and $\hat{K} = Q^T A_2 (Q^T)^{-1}$. By our choice of $Q$, both $\hat{K}$ and $K$ are
real symmetric negative definite matrices. From the continuity of the vector field, we have

$$
(\hat{K} - K)x + Q^T b' = 0 \quad \text{when} \quad x_1 = d. \tag{4}
$$

Let $\hat{K} = [\hat{k}_1 \ \hat{k}_2 \ \dots \ \hat{k}_N]$ and $K = [k_1 \ k_2 \ \dots \ k_N]$. Then (4) reads

$$
d(\hat{k}_1 - k_1) + x_2(\hat{k}_2 - k_2) + \dots + x_N(\hat{k}_N - k_N) + Q^T b' = 0,
$$

which most hold $\forall x_i \in \mathbb{R}, i = 2, \dots, N$. Thus, we arrive at

$$
\hat{k}_1 = k_1 + \hat{k}, \quad \hat{k}_i = k_i, \quad i = 2, \dots, N,
$$

where $\hat{k} = -(1/d)Q^T b'$. Therefore, (3) can be expressed as

$$
\ddot{x} = \begin{cases} Kx & \text{for } x_1 \le d \\ Kx + \hat{k}(x_1 - d) & \text{for } x_1 > d \end{cases}, \quad (5)
$$
or in a more compact form

$$ \ddot{x} = \begin{cases} Kx & \text{for } x_1 \le d \\ \hat{K}(x - x_e) & \text{for } x_1 > d \end{cases}, \qquad (6) $$

where $x_e = d\hat{K}^{-1}\hat{k}$. Moreover, by the symmetry of $\hat{K}$, $\hat{k}$ must be of the form

$$ \hat{k} = [\kappa \ 0 \ \dots \ 0]^T. $$

Either equation (5) or equation (6) can be considered as the canonical form of the PWL system (3) with $\kappa$ as a parameter, and when $\kappa = 0$, they reduce to a purely linear system. Also, for $d = 0$, the system reduces to that discussed in [41].

It will be clear in Section 4 that as far as $\hat{K}$ is negative definite, $\kappa$ is not restricted in size in our analysis. However, since the method of asymptotic expansion is used, the magnitude of $\kappa$ might affect the region of validity of the expansions. This will be discussed in Section 6.

### 3. The Poincaré Map

In this section, we derive the Poincaré map for the PWL system in canonical form (6), which will be used to construct the nonlinear normal mode in the next section. The switching hyperplane provides a natural Poincaré section. To this aim, we make use of the symmetry of $K$ and $\hat{K}$, and the linearity of the subsystems. Since both subsystems are linear, we can obtain their solutions in closed form which enables one to obtain an analytical, although implicit, expression for the Poincaré map.

Let us start with the derivation of the solutions of the linear subsystems. Let $\{-\omega_i^2, v_i\}_{i=1}^N$ and $\{-\hat{\omega}_i^2, \hat{v}_i\}_{i=1}^N$ be the eigenpairs of $K$ and $\hat{K}$, respectively. Suppose that the eigenvectors are normalized in such a way that they have magnitude one and positive first entries, i.e., let $v_i = [v_{1i} \ v_{2i} \ \dots \ v_{Ni}]^T$ and $\hat{v}_i = [\hat{v}_{1i} \ \hat{v}_{2i} \ \dots \ \hat{v}_{Ni}]^T$, then $v_{1i} \ge 0$ and $\hat{v}_{1i} \ge 0$. Then the solution of the first linear subsystem is

$$ x(t) = \sum_{i=1}^{N} v_i (a_i \cos \omega_i t + b_i \sin \omega_i t) = V(C(t)a + S(t)b); \qquad (7) $$

$$ y(t) := \dot{x}(t) = \sum_{i=1}^{N} \omega_i v_i (-a_i \sin \omega_i t + b_i \cos \omega_i t) = V\Omega(-S.ta + C(t)b), \qquad (8) $$

where $V = [v_1 \ v_2 \ \dots \ v_N]$ is orthogonal since $K$ is real symmetric, $C(t) = \text{diag}(\cos \omega_1 t, \cos \omega_2 t, \dots, \cos \omega_N t)$, $S(t) = \text{diag}(\sin \omega_1 t, \sin \omega_2 t, \dots, \sin \omega_N t)$, $\Omega = \text{diag}(\omega_1, \omega_2, \dots, \omega_N)$, and $a = [a_1 \ a_2 \ \dots \ a_N]^T$, $b = [b_1 \ b_2 \ \dots \ b_N]^T \in \mathbb{R}^N$ are coefficients determined by the initial conditions. Setting $t = 0$ in equations (7) and (8), we have

$$ x(0) = Va, \quad y(0) = V\Omega b, $$

which leads to

$$ a = V^{-1}x(0) = V^T x(0), \quad b = \Omega^{-1}V^T y(0). \qquad (9) $$

Inserting (9) back into (7) and (8) one obtains

$$ x(t) = VC(t)V^T x(0) + VS(t)\Omega^{-1}V^T y(0); $$

$$ y(t) = -V\Omega S(t)V^T x(0) + V\Omega C(t)\Omega^{-1}V^T y(0); $$
or in an integrated form

$$w(t) = T(t)w(0), \qquad (10)$$

where $w(t) = [x(t) \ y(t)]^T$ and

$$T(t) = \begin{bmatrix} V(C)(t)V^T & VS(t)\Omega^{-1}V^T \\ -V\Omega S(t)V^T & VC(t)V^T \end{bmatrix},$$

is a $2N \times 2N$ matrix composed of four $N \times N$ matrices. Note that we have used the commutative property of two diagonal matrices to simplify the expression. Furthermore, we can put $T(t)$ in the clearer form

$$T(t) = \begin{bmatrix} V & 0 \\ 0 & V \end{bmatrix} \begin{bmatrix} C(t) & S(t)\Omega^{-1} \\ -\Omega S(t) & C(t) \end{bmatrix} \begin{bmatrix} V^T & 0 \\ 0 & V^T \end{bmatrix}. \qquad (11)$$

It is not difficult to see that $T(t)$ is invertible and that its inverse is given by

$$T(t)^{-1} = \begin{bmatrix} V & 0 \\ 0 & V \end{bmatrix} \begin{bmatrix} C(t) & -S(t)\Omega^{-1} \\ \Omega S(t) & C(t) \end{bmatrix} \begin{bmatrix} V^T & 0 \\ 0 & V^T \end{bmatrix}. \qquad (12)$$

For the solution of the second linear subsystem, a similar procedure gives

$$w(t) = w_e + \hat{T}(t)(w(0) - w_e), \qquad (13)$$

where $w_e = [x_e \ 0]^T$ and $\hat{T}(t)$ has a definition similar to $T(t)$ and can be obtained from $T(t)$ by replacing $\omega_i, v_i$ with $\hat{\omega}_i, \hat{v}_i$. Note that for $\kappa = 0$, $\hat{T}(t) = T(t)$ and $x_e = 0$, and equation (13) reduces to equation (10).

We are now in a position to construct the Poincaré map. Define the surfaces of section (here, hyperplanes)

$$\Sigma = \{w \in \mathbb{R}^{2N} : x_1 = d \text{ and } y_1 > 0\}$$

and

$$\Sigma^* = \{w \in \mathbb{R}^{2N} : x_1 = d \text{ and } y_1 < 0\}.$$

The Poincaré map $F: \Sigma \to \Sigma$ is thus defined as

$$F(w) = T(t_{h2})(w_e + \hat{T}(t_{h1})(w - w_e)). \qquad (14)$$

$\forall w \in \Sigma$, where $t_{h1}$ is the time required to take $w$ to a point $w^*$ on $\Sigma^*$ and $t_{h2}$ is that required to return from $w^*$ back to $\Sigma$. In general, both $t_{h1}$ and $t_{h2}$ satisfy transcendental algebraic equations and only approximate solutions are possible (see [11, 27]). Note that although the system is piecewise linear, the corresponding Poincaré map is smooth.

# 4. Normal Modes for the PWL System

## 4.1. CONSTRUCTION OF THE NORMAL MODES

With the preliminary results developed in previous sections, the procedure of constructing normal modes for the PWL system (6) by the methods of invariant manifolds and asymptotic
<figure><img src="image_2.png" /><figcaption>Figure 1. The linear normal mode.</figcaption></figure>

expansions is provided in this section. Since the systems considered here are conservative, it is well known that each normal mode is a collection of periodic motions [28]. In other words, we are seeking invariant manifolds in state space that are composed of one-parameter families of periodic orbits. Therefore, when viewed in the Poincaré map, the invariant manifold is a one-parameter family of fixed points.

In the linear case, that is, when $\kappa = 0$ in system (6), a normal mode is a two-dimensional plane in state space, which intersects the switching hyperplane along a straight line. This line constitutes the one-parameter family of fixed points for the Poincaré map (14) for the corresponding normal mode. To fix the ideas, let us take the first mode as an example. Let the corresponding eigenvector be written as $v_1 = (v_{11}/d)[d \ u_2 \ ... \ u_N]^T$, where $u_i = dv_{i1}/v_{11}$, $i = 2, ..., N$. Then the 2-D invariant plane for this mode is:

$$ \{w = [x \ y]^T \in \mathbb{R}^{2N} : x = \alpha v_1, y = \beta v_1, \alpha, \beta \in \mathbb{R}\}, $$

which is composed of a family of ellipses that are periodic solutions with different amplitudes but the same frequency $\omega_1$, as shown in Figure 1. Figure 2 shows the corresponding fixed point family for the Poincaré map, given by

$$ \{w = [x \ y]^T \in \mathbb{R}^{2N} : x = \frac{d}{v_{11}} v_1, y = \beta v_1, \beta \geq 0\}. $$

For $\kappa \neq 0$, the 2-D invariant plane in state space as well as the one-parameter family of fixed points on the Poincaré section will be distorted. The fixed point family is no longer a straight line as in the linear case, but a curve emanating from a point $w_t = [(d/v_{11})v_1 \ 0]^T$ (see Figure 2). Correspondingly, the invariant plane will become an invariant manifold consisting of two pieces. One piece is an elliptic disk near the origin which represents the small amplitude vibration and is tangent to the switching hyperplane at $w_t$. This piece is a portion of the linear normal mode of the first linear subsystem since near the origin the system is linear. The other piece represents large amplitude vibration and involves both linear subsystems. This piece again contains two parts, one in the region $x_1 \leq d$ and the other in $x_1 > d$. The two parts must match along the Poincaré section in the sense that they meet at the curve of fixed points. Figure 3 is a qualitative sketch of such an invariant manifold.

In order to obtain the normal mode, we only need to construct the large amplitude piece since the elliptic disk is already known. Because it is generally impossible to obtain the entire invariant manifold in closed form, we shall employ the method of asymptotic expansions which allows us to obtain an approximate invariant manifold in the neighborhood of the
<figure><img src="image_3.png" /><figcaption>Figure 2. The fixed point family of the Poincaré map for the linear (dashed line) and PWL (solid line) systems.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Sketch of an invariant manifold for the PWL system.</figcaption></figure>

elliptic disk. For this purpose, it is natural to use polar coordinates, in contrast with the Cartesian coordinates used by Shaw and Pierre [28].

There are several steps involved in the procedure for constructing a normal mode for a PWL systems, and the procedure is exactly the same for each mode. First, an approximation to the curve of fixed points is determined. Next, using the curve obtained as a matching condition, the large amplitude piece of the invariant manifold, which is in the neighborhood of the elliptic disk and belongs to the region $x_1 \le d$, is obtained in a series form of polar coordinates. Finally, by matching the curve from the first step again, and by making use of another set of polar coordinates, the other portion of the invariant manifold belonging to the region $x_1 > d$ is expanded in the neighborhood of $w_t$. Then, by gluing together these two pieces and the elliptic disk from the linear normal mode, we are able to construct an approximation of the entire nonlinear normal mode for the PWL system. In what follows, we will explain these steps in some detail. Again, we take the first mode as the illustrative example, where the extension to the other modes is straightforward.

The fixed points of the Poincaré map are those on $\Sigma$ satisfying

$$F(w) = w, \tag{15}$$

subject to

$$\{w_e + \hat{T}(t_{h1})(w - w_e)\}_1 = d, \tag{16}$$

where $\{\cdot\}_i$ represents the $i$-th element of $\{\cdot\}$. In general, equation (15) together with equation (16) comprises a set of $2N + 1$ transcendental algebraic equations in $2N + 1$ variables, i.e., $x_2, \dots, x_N, y_1, \dots, y_N, t_{h1}$ and $t_{h2}$.
In order to determine the one-parameter family of fixed points corresponding to the first mode, we can assume a solution in a series in terms of the velocity of the first degree of freedom at the switching point, $y_1$, up to m-th order as the following

$$y_1 = \delta; \tag{17}$$

$$x_i = u_i + \sum_{j=1}^{M} \mu_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N; \tag{18}$$

$$y_i = \sum_{j=1}^{m} \nu_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N; \tag{19}$$

$$t_{h1} = \sum_{j=1}^{m} \tau_{1j} \delta^j + O(m+1); \tag{20}$$

$$t_{h2} = \frac{2\pi}{\omega_1} + \sum_{j=1}^{m} \tau_{2j} \delta^j + O(m+1), \tag{21}$$

where $O(m+1)$ denotes the higher order terms and the $\mu_{ij}$'s, $\nu_{ij}$'s, and $\tau_{ij}$'s are unknown coefficients, depending on the parameter $\kappa$, which are to be determined. Recall that $u_i = dv_{i1}/v_{11}$ and $[d u_2 \dots u_N]^T$ is an eigenvector associated with the first mode of the first linear subsystem. Note here that we have assumed $v_{11} \neq 0$. In general, if $v_{1i} = 0$ for some $1 \le i \le N$, then the i-th normal mode of the PWL system is the same as that of the first linear subsystem, i.e., a 2-D plane. This is because on this mode, $x_1 \equiv 0$ and hence it will never exceed $d$. Note also that as a $\delta \to 0$ in equations (17) through (21), the assumed solution reduces to $w_t$, which is a fixed point for all values of $\kappa$.

Substituting equations (17)–(21) into equations (15) and (16), and matching the coefficients of the same order of $\delta$, we obtain $m(2N+1)$ equations. In doing this, we need power series expansions of $\hat{T}(t_{h1})$ and $T(t_{h2})$ about $\delta = 0$ up to order $m$. Then the solution of the $\mu_{ij}$'s, $\nu_{ij}$'s, and $\tau_{ij}$'s are obtained by solving these equations sequentially. This completes the first step in the process. Notice that there are $m(2N+1)$ equations and only $2mN$ unknowns, although some of the equations are redundant.

Before moving to the next step, we transform the coordinates $(x_1, y_1)$ to a polar form $(r, \theta)$ defined by

$$x_1 = r \cos \theta, \quad y_1 = \omega_1 r \sin \theta. \tag{22}$$

As a result, a periodic solution of the first mode of the first linear subsystem can be represented simply as $r = \text{constant}$ when projected onto $(x_1, y_1)$ plane. Then when $r \le d$, the normal mode can be expressed as

$$x_i = \frac{u_i}{d} r \cos \theta, \quad 2 \le i \le N;$$

$$y_i = \frac{u_i}{d} \omega_1 r \sin \theta, \quad 2 \le i \le N,$$

which is the elliptic disk mentioned above. Generally, in terms of $(r, \theta)$, the normal mode for the PWL system takes the form
$$x_i = \begin{cases} \frac{u_i}{d} r \cos \theta & \text{for } r \le d \\ \frac{u_i}{d} r \cos \theta + f_{i1}(r, \theta) & \text{for } r > d \text{ and } r \cos \theta \le d; \\ \frac{u_i}{d} r \cos \theta + f_{i2}(r, \theta) & \text{for } r \cos \theta > d \end{cases} \quad (23)$$

$$y_i = \begin{cases} \frac{u_i}{d} \omega_1 r \sin \theta & \text{for } r \le d \\ \frac{u_i}{d} \omega_1 r \sin \theta + g_{i1}(r, \theta) & \text{for } r > d \text{ and } r \cos \theta \le d, \\ \frac{u_i}{d} \omega_1 r \sin \theta + g_{i2}(r, \theta) & \text{for } r \cos \theta > d \end{cases} \quad (24)$$

for $2 \le i \le N$. Basically, our aim in this section is to find approximate solutions for $f_{ij}(r, \theta)$ and $g_{ij}(r, \theta)$. Clearly, these functions must satisfy the continuity condition on the boundaries.

Next, consider the region $r > d$ and $r \cos \theta \le d$, which is outside the elliptic disk but still governed by the first linear subsystem. Let $P(\delta)$ be a fixed point on $\Sigma$ with $y_1$ coordinate being $\delta$ and others being given from the first step in equations (18) and (19). Then there must be a periodic trajectory passing through $P$.<fnref n="1" /> Let $P_{\theta}(\delta)$ be a point on the trajectory in this region whose coordinates are

$$x_1 = r \cos \theta, \quad x_i = \frac{u_i}{d} r \cos \theta + \sum_{j=1}^{m} \alpha_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N; \quad (25)$$

$$y_1 = \omega_1 r \sin \theta, \quad y_i = \frac{u_i}{d} \omega_1 r \sin \theta + \sum_{j=1}^{m} \beta_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N, \quad (26)$$

where the $\alpha_{ij}$'s and $\beta_{ij}$'s are functions of $\theta$ to be determined, $r$ is a function of $\delta$ and $\theta$ and is assumed to be

$$r = d + \sum_{j=1}^{m} r_j \delta^j + O(m+1), \quad (27)$$

where the $r_j$'s are also functions of $\theta$ to be determined. In other words, $P_{\theta}$ is a point on the normal mode. Reference to Figure 4, which is the state space projected onto the $(x_1, y_1)$ plane, may clarify matters for the reader at this point. If $\delta = 0$ in equations (25)–(27), $P_{\theta}$ becomes a point on the boundary of the elliptic disk and will reach the Poincaré section $\Sigma$ at time $t = \theta/\omega_1$. Therefore, it is natural to express the hitting time $t_{h\theta}$ of $P_{\theta}$ to $P$ as

$$t_{h\theta} = \frac{\theta}{\omega_1} + \sum_{j=1}^{m} t_j \delta^j + O(m+1), \quad i = 2, \dots, N, \quad (28)$$

where the $t_j$'s are again functions of $\theta$ to be determined.

With $P_{\theta}$ as the initial point $w(0)$ and substituting expressions (25)–(28) into equation (10), we obtain $w(t_{h\theta}) = P$ or

$$P = T(t_{h\theta}) P_{\theta}. \quad (29)$$

By expanding $T(t_{h\theta})$ with respect to $\delta = 0$ up to order $m$ and matching coefficients of like order in $\delta$ in equation (29), one can solve for the $\alpha_{ij}$'s, $\beta_{ij}$'s, $r_j$'s and $t_j$'s. Then, by inverting

<footnote n="1">¹ From now on, $P$ instead of $P(\delta)$ will be used for brevity, similarly for $P_{\theta}(\delta)$ and $P_{\phi}(\delta)$ in the sequel.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 4. Polar coordinates for the region $x_1 \le d$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Polar coordinates for the region $x_1 \ge d$.</figcaption></figure>

the power series of equation (27), $\delta$ can be obtained as an approximate function of $r$ and $\theta$. Entering the resulting $\delta(r, \theta)$ into equations (25) and (26), and comparing the result with equations (23) and (24), we have

$$f_{i1}(r, \theta) \approx \sum_{j=1}^{m} \alpha_{ij}(\theta) \delta(r, \theta)^j; \qquad (30)$$

$$g_{i1}(r, \theta) \approx \sum_{j=1}^{m} \beta_{ij}(\theta) \delta(r, \theta)^j, \qquad (31)$$

which completes the second step.

Finally, to build up the piece of invariant manifold in the region $r \cos\theta > d$, we need another polar coordinate set $(\hat{r}, \hat{\phi})$ for $(x_1, y_1)$ as shown in Figure 5. Let $P_{\phi}(\delta)$ be a point on the periodic orbit passing through $P$. In view of the fact that $P_{\phi}$ is in the neighborhood of $w_t$, the coordinates of $P_{\phi}$ are assumed to be

$$x_1 = d + \hat{r} \cos \phi, \quad x_i = \frac{u_i}{d} x_1 + \sum_{j=1}^{m} \hat{\alpha}_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N; \qquad (32)$$
$$y_1 = \omega_1 \hat{r} \sin \phi, \quad y_i = \frac{u_i}{d} y_1 + \sum_{j=1}^{m} \hat{\beta}_{ij} \delta^j + O(m+1), \quad i = 2, \dots, N, \quad (33)$$

and expand $\hat{r}$ about $w_t$ as follows

$$\hat{r} = \sum_{j=1}^{m} \hat{r}_j \delta^j + O(m+1). \quad (34)$$

If $\delta = 0$ in equations (32) through (34), then $\hat{r} = 0$, and both $P_\phi$, and $P$ shrink to the same point, $w_t$. Again, the backward hitting time $t_{h\phi}$ for $P_\phi$ to $P$ is expressed as a series in $\delta$ as follows

$$t_{h\phi} = \sum_{j=1}^{m} \hat{t}_j \delta^j + O(m+1). \quad (35)$$

As in the procedure described above, the expansion coefficients $\hat{\alpha}_{ij}$'s, $\hat{\beta}_{ij}$'s, $\hat{r}_j$'s and $\hat{t}_j$'s are unknown functions of $\phi$ determined by the matching condition

$$P = w_e + \hat{T}(-t_{h\phi})(P_\phi - w_e). \quad (36)$$

Transforming $(\hat{r}, \phi)$ back to $(r, \theta)$, and from equation (34), we solve for $\delta$ as a function of $(r, \theta)$. Then, comparing equations (32) and (33) with equations (23) and (24), we obtain

$$f_{i2}(r, \theta) \approx \sum_{j=1}^{m} \hat{\alpha}_{ij}(r, \theta) \delta(r, \theta)^j; \quad (37)$$

$$g_{i2}(r, \theta) \approx \sum_{j=1}^{m} \hat{\beta}_{ij}(r, \theta) \delta(r, \theta)^j, \quad (38)$$

where one should note that $\delta(r, \theta)$ here is a different function from that given in equations (30) and (31) since they are in different regions of $(r, \theta)$. This is the final step and it completes the procedure for generating the normal modes of the PWL system (6).

This section is ended with some remarks on the order of approximation. Since we are expanding the invariant manifold around an invariant disk, solutions will be accurate close to the disk. Here the measurement of closeness is $\delta$, the $y_1$ coordinate of a fixed point on the Poincaré section $\Sigma$. Throughout the paper, the order of approximation is in terms of $\delta$. One may relate $\delta$ to the energy level of the modal motion above and near the switching hyperplane. A larger $\delta$ represents a higher energy level on the invariant manifold. However, their relationship cannot be obtained in closed form since a closed form expression for the invariant manifold is not available.

## 4.2. THE EXISTENCE AND UNIQUENESS OF THE NORMAL MODES

The problem of determining the existence and uniqueness of nonlinear normal modes has been the subject of many efforts since the concept of nonlinear normal modes was introduced [5, 15, 19, 40]. For smooth nonlinear systems, the existence problem is well understood. As for uniqueness, the Lyapunov center theorem states that the number of nonlinear normal modes is the same as that of linear modes under nonresonant conditions [30]. In the case of internal
resonance, Nayfeh and Chin reported in [14] that the number of nonlinear normal modes may exceed the number of linear ones. For PWL systems, both the existence and uniqueness problems remain open. In [41], Zuo and Curnier numerically obtained the normal modes for a special class of PWL systems and discussed the bifurcation of these modes. However, they did not address the existence problem, nor did they investigate internally resonant cases.

Although we do not have a definite answer to the problem of existence and uniqueness of normal modes for PWL systems either, some evidence (a necessary condition) to the existence of the normal modes has been determined. Specifically, the problem is related to the solution of the fixed point equations (15) and (16). If, for a given linear mode, there is a unique one parameter family of fixed points satisfying equations (15) and (16), then the nonlinear normal mode exists and is unique for the given mode. Suppose this is indeed the case. Then there must be exactly $2N$ independent equations in the set of $2N + 1$ simultaneous equations (equations (15) and (16)), i.e., essentially, only one of the equations in (15) and (16) depends on the others. This is verified by the illustrative example following in Section 6. In other words, the Jacobian matrix of equations (15) and (16) at the given linear mode with $\kappa = 0$ and $\delta = 0$ should have rank equal to $2N$. It is shown in Appendix A that this is generally true for nonresonant cases. In lieu of equations (15) and (16), a more convenient but equivalent form of fixed point equations is used in Appendix A. Under resonant conditions, the rank of the Jacobian matrix will be less than $2N$, indicating the possibility of nonunique normal modes for a given linear mode. Interestingly, one should note that the PWL system is linear locally and hence the normal mode is unique near the origin, whether resonant or not. Over the switching hyperplane, however, the normal mode can split into several pieces.

## 5. Modal Dynamics for the PWL System

With the normal modes in hand, the next step is to establish the dynamics of the normal modes. As expected, the modal dynamics are also of piecewise type.

Again, the first mode is considered here. Its local coordinate system is taken to be the polar form $(r, \theta)$ defined by equation (22). Hence

$$ \dot{r} = \dot{x}_1 \cos \theta + \frac{1}{\omega_1} \dot{y}_1 \sin \theta; \qquad (39) $$

$$ \dot{\theta} = -\frac{1}{r} \dot{x}_1 \sin \theta + \frac{1}{\omega_1 r} \dot{y}_1 \cos \theta. \qquad (40) $$

From the state equations, $\dot{x}_1$ and $\dot{y}_1$ are functions of $x$ and $y$, which are related to $r$ and $\theta$ on the first mode in piecewise fashion by equations (22)–(24). From equation (5), $\dot{x}_1$ and $\dot{y}_1$ are given in terms of $r$ and $\theta$ by

$$ \dot{x}_1 = y_1 = \omega_1 r \sin \theta; \qquad (41) $$

$$ \dot{y}_1 = \begin{cases} \{Kx\}_1 & \text{for } x_1 \le d \\ \{Kx\}_1 + \kappa(x_1 - d) & \text{for } x_1 > d \end{cases} \\ = \begin{cases} \{K\tilde{v}_1 r \cos \theta\}_1 & \text{for } r \le d \\ \{K(\tilde{v}_1 r \cos \theta + f_1(r, \theta))\}_1 & \text{for } r > d \text{ and } r \cos \theta \le d, \\ \{K(\tilde{v}_1 r \cos \theta + f_2(r, \theta))\}_1 + \kappa(r \cos \theta - d) & \text{for } r \cos \theta > d \end{cases} $$
where $\tilde{v}_1 = [1 \ (u_2/d) \ \dots \ (u_N/d)]^T$ is the eigenvector of $K$ associated with eigenvalue $-\omega_1^2$ and $f_i(r, \theta) = [0 \ f_{2i}(r, \theta) \ \dots \ f_{Ni}(r, \theta)]^T$ for $i = 1$ and $2$ are as defined in equation (23) and their approximate solutions are given in equations (30) and (37). Also, the reader should recall that $\{\cdot\}_1$ stands for the first element of $\{\cdot\}$. Therefore,

$$ \{K \tilde{v}_1 r \cos \theta\}_1 = \{-\omega_1^2 \tilde{v}_1 r \cos \theta\}_1 = -\omega_1^2 r \cos \theta, $$

which yields

$$ \dot{y}_1 = \begin{cases} -\omega_1^2 r \cos \theta & \text{for } r \le d \\ -\omega_1^2 r \cos \theta + k_1^T f_1(r, \theta) & \text{for } r > d \text{ and } r \cos \theta \le d, \\ -\omega_1^2 r \cos \theta + k_1^T f_2(r, \theta) + \kappa(r \cos \theta - d) & \text{for } r \cos \theta > d \end{cases}, \quad (42) $$

where $k_1$ is the first column of $K$. Inserting expressions (41) and (42) into equations (39) and (40), one obtains the dynamics of the first mode as

$$ \dot{r} = \begin{cases} 0 & \text{for } r \le d \\ \frac{\sin \theta}{\omega_1} [k_1^T f_1(r, \theta)] & \text{for } r > d \text{ and } r \cos \theta \le d; \\ \frac{\sin \theta}{\omega_1} [k_1^T f_2(r, \theta) + \kappa(r \cos \theta - d)] & \text{for } r \cos \theta > d \end{cases} \quad (43) $$

$$ \dot{\theta} = \begin{cases} -\omega_1 & \text{for } r \le d \\ -\omega_1 + \frac{\cos \theta}{\omega_1 r} [k_1^T f_1(r, \theta)] & \text{for } r > d \text{ and } r \cos \theta \le d. \\ -\omega_1 + \frac{\cos \theta}{\omega_1 r} [k_1^T f_2(r, \theta) + \kappa(r \cos \theta - d)] & \text{for } r \cos \theta > d \end{cases} \quad (44) $$

As expected, the modal dynamics are of piecewise type and can be regarded as a one degree of freedom nonlinear oscillator. Also, it is easy to see that the right-hand sides of (42)–(44) are continuous. Note that the first piece ($r \le d$) is for motions of the first linear subsystem, which govern for energy levels below that required to exceed the switching hyperplane, while the second and third pieces describe the motion in the two linear domains for higher energy levels. Due to the piecewise nature of the equations, however, conventional methods in determining the frequency-amplitude relationship of nonlinear oscillators fail to apply here. Conventional methods make use of perturbation techniques which require the system to be smooth, which is not the case here. Despite these difficulties, we can still determine the frequency-amplitude relationship of the modal dynamics from the analysis presented in Section 4, rather than deal with equations (43) and (44) directly.

Consider the motion on the normal mode. When the amplitude is sufficiently small that the vibration stays in the disk region $r \le d$, the system is linear with constant frequency $\omega_1$. If the vibration is large enough, it will hit the switching hyperplane at a fixed point on the Poincaré section $\Sigma$. Each fixed point is parameterized by the $y_1$ coordinate and is given by equations (17)–(19) in Section 4. For each $y_1 = \delta$, we have one such oscillation whose period is equal to $t_p = t_{h1} + t_{h2}$, which is also parameterized by $\delta$, as given in equations (20) and (21). Specifically,

$$ t_p = \frac{2\pi}{\omega_1} + \sum_{j=1}^{m} t_{pj} \delta^j + O(m+1), $$
<figure><img src="image_7.png" /><figcaption>Figure 6. The mechanical model for the example.</figcaption></figure>

where $t_{pj} = \tau_{1j} + \tau_{2j}$. Thus, the frequency can be easily calculated from

$$ \tilde{\omega}_1 = \frac{2\pi}{t_p} = \omega_1 + \sum_{j=1}^{m} \tilde{\omega}_{1j}\delta^j + O(m+1). \quad (45) $$

Equation (45) is the frequency-amplitude relationship of the first modal dynamics in terms of $\delta$. We can also convert equation (45) to the usual form of relationship between $\tilde{\omega}_1$ and $A_1$, the peak displacement of $x_1$. To do this, let $\bar{w} = [\bar{x} \ \bar{y}]^T$ be a fixed point on $\Sigma$ with $\bar{y}_1 = \delta$. The peak displacement occurs when the velocity is zero, i.e., $y_1 = 0$. With $\bar{w}$ as the initial point, we thus seek the value of $x_1$ at $y_1 = 0$. First, the time $t_{h0}$ bringing $y_1 = \delta$ to $y_1 = 0$ is obtained in terms of a power series in $\delta$ by solving

$$ \{w_e + \hat{T}(t_{h0})(\bar{w} - w_e)\}_{N+1} = 0 $$

approximately up to the $m$-th order of $\delta$. Then the peak displacement of $x_1$ is simply given by

$$ A_1 = \{w_e + \hat{T}(t_{h0})(\bar{w} - w_e)\}_1, \quad (46) $$

which is also a series form of $\delta$. Using equation (46), one can obtain $\delta$ as a function of $A_1$ in a series in terms of $\delta$. The frequency-amplitude dependence in terms of $A_1$ is then expressed as

$$ \tilde{\omega}_1 \approx \begin{cases} \omega_1 & \text{for } A_1 \le d \\ \omega_1 + \sum_{j=1}^{m} \tilde{\omega}_{1j}\delta(A_1)^j & \text{for } A_1 > d \end{cases} \quad (47) $$

These ideas are now demonstrated by an example.

## 6. Example: A Two Degrees of Freedom System with a Clearance

The system to be considered is shown in Figure 6, which is an undamped, two degrees of freedom system with a clearance. Except for the clearance, no other nonlinearities appear in the system. Following the procedure provided above, the normal modes for the system, the associated modal dynamics and the frequency-amplitude relationship will be given for both normal modes. Numerical simulations are also carried out to compare with the theoretical results.

It is important to point out that most of the manipulations encountered here, such as power series expansions, inversion of power series, and solving for the unknown coefficients of the normal modes, can hardly be done without the aid of computer symbolic manipulators. Indeed, much of the work for this example was done by Mathematica™ running on a Macintosh II
cx computer. All of the numerical simulations presented were carried out using MatLab on a
Sun station.

Assume that the springs are unstressed at the equilibrium position, where the displacements $x_1$ and $x_2$ are zero. In addition, all the values for inertia, stiffnesses, and clearance magnitude are taken to be unity except the free spring whose stiffness is assumed to be $\kappa$. The equations of motion are, therefore,

$$ \ddot{x}_1 = \begin{cases} -2x_1 + x_2 & \text{for } x_1 \le 1 \\ -2x_1 + x_2 - \kappa(x_1 - 1) & \text{for } x_1 > 1 \end{cases}; \quad (48) $$

$$ \ddot{x}_2 = x_1 - 2x_2. \quad (49) $$

In terms of the canonical form given in equation (6), the equations of motion are expressed with the following definitions:

$$ d = 1, $$

$$ x_e = \left( \frac{\frac{2\kappa}{3+2\kappa}}{\frac{\kappa}{3+2\kappa}} \right), $$

$$ K = \begin{bmatrix} -2 & 1 \\ 1 & -2 \end{bmatrix}, $$

$$ \hat{K} = \begin{bmatrix} -2 - \kappa & 1 \\ 1 & -2 \end{bmatrix}. $$

The eigenpairs of $K$ are

$$ -\omega_1^2 = -1, \quad v_1 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ 1 \end{pmatrix}; $$

$$ -\omega_2^2 = -3, \quad v_2 = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ -1 \end{pmatrix}, $$

and those for $\hat{K}$ are

$$ -\hat{\omega}_1^2 = \lambda_1 - 2, \quad \hat{v}_1 = \frac{1}{\sqrt{1 + \lambda_1^2}} \begin{pmatrix} 1 \\ \lambda_1 \end{pmatrix}; $$

$$ -\hat{\omega}_2^2 = \lambda_2 - 2; \quad \hat{v}_2 = \frac{1}{\sqrt{1 + \lambda_2^2}} \begin{pmatrix} 1 \\ \lambda_2 \end{pmatrix}, $$

where

$$ \lambda_1 = \frac{-\kappa + \sqrt{4 + \kappa^2}}{2} \quad \text{and} \quad \lambda_2 = \frac{-\kappa - \sqrt{4 + \kappa^2}}{2}. $$

Thus, we have

$$ V = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}, $$
$$ \Omega = \begin{bmatrix} 1 & 0 \\ 0 & \sqrt{3} \end{bmatrix}; $$

$$ \hat{V} = \begin{bmatrix} \frac{1}{\sqrt{1+\lambda_1^2}} & \frac{1}{\sqrt{1+\lambda_2^2}} \\ \frac{\lambda_1}{\sqrt{1+\lambda_1^2}} & \frac{\lambda_2}{\sqrt{1+\lambda_2^2}} \end{bmatrix}, $$

$$ \hat{\Omega} = \begin{bmatrix} \sqrt{2-\lambda_1} & 0 \\ 0 & \sqrt{2-\lambda_2} \end{bmatrix}. $$

## 6.1. THE FIXED POINT FAMILY ON THE POINCARÉ SECTION

Consider the first mode and note that $\omega_1 = 1$, $d = 1$ and $u_2 = 1$ for our example. We express the fixed point $w = [x \ y]^T$ on $\Sigma$ as

$$ x_1 = 1, $$

$$ x_2 = 1 + \mu_{21}\delta + \mu_{22}\delta^2 + \mu_{23}\delta^3 + O(4), \qquad (50) $$

$$ y_1 = \delta, $$

$$ y_2 = \nu_{21}\delta + \nu_{22}\delta^2 + \nu_{23}\delta^3 + O(4), \qquad (51) $$

and the hitting times as

$$ t_{h1} = \tau_{11}\delta + \tau_{12}\delta^2 + \tau_{13}\delta^3 + O(4); \qquad (52) $$

$$ t_{h2} = 2\pi + \tau_{21}\delta + \tau_{22}\delta^2 + \tau_{33}\delta^3 + O(4), \qquad (53) $$

where the approximations are done up to third order. The coefficients in equations (50)–(53) have to be chosen to satisfy the fixed point equations (15) and (16). Instead of solving equations (15) and (16) simultaneously, which requires some cumbersome calculations, we divide the job into several steps. First, $t_{h1}$ is solved (in terms of other unknowns) from equation (16). Next, the hitting point $w^*$ of $w$ on $\Sigma^*$ is obtained in terms of the unknowns. Finally, $w^*$ is put into equation (15) and solved. The resulting equations for $t_{h1}$ are

$$ \delta \text{ term : } \quad 1 - \frac{1}{2} \tau_{11} = 0, $$

$$ \delta^2 \text{ term : } \quad \mu_{21}\tau_{11} - \tau_{12} = 0, $$

$$ \delta^3 \text{ term : } \quad \frac{1}{2}(\mu_{22}\tau_{11} + \mu_{21}\tau_{12} - \tau_{13}) + \frac{1}{6}\tau_{11}^2(\nu_{21} - 2 - \kappa) + \frac{1}{24}(1 + \kappa)\tau_{11}^3 = 0, $$

which yield

$$ \tau_{11} = 2, \qquad (54) $$

$$ \tau_{12} = 2\mu_{21}, \qquad (55) $$

$$ \tau_{13} = 2\mu_{22} + 2\mu_{21}^2 + \frac{2}{3}(2\nu_{21} - 3 - \kappa). \qquad (56) $$
Then, by inserting $t_{h1}$ into equation (13), $w^*$ can be expressed as

$$x_2^* = 1 + \mu_{21}\delta + (\mu_{22} - 2 + 2\nu_{21})\delta^2 + (\mu_{23} - 8\mu_{21} + 2\nu_{22} + 2\mu_{21}\nu_{21})\delta^3 + O(4), \quad (57)$$

$$y_1^* = -\delta + \frac{2}{3}(\nu_{21} - 1)\delta^3 + O(4), \quad (58)$$

$$y_2^* = (\nu_{21} - 2)\delta + (\nu_{22} - 6\mu_{21})\delta^2 + \left[ \nu_{23} - 6(\mu_{21}^2 + \mu_{22} - 1) - \frac{2}{3}(8\nu_{21} - \kappa + 1) \right] \delta^3 + O(4). \quad (59)$$

Equations (57)–(59) are significantly more complicated if $t_{h1}$ is not obtained first. From equation (58), one can see that $y_1^* < 0$ when $\delta$ is sufficiently small and that $w^*$ is indeed on $\Sigma^*$.

Substituting equations (57)–(59) into equation (15) results in the following 12 equations with 9 unknowns:

$$\delta \text{ term :} \quad \frac{1-c}{2} \mu_{21} - \frac{s}{2\sqrt{3}} (\nu_{21}-1) = 0, \quad (60)$$

$$\mu_{21} = \mu_{21}, \quad (61)$$

$$-1 - \tau_{21} = \nu_{21} \quad (62)$$

$$\nu_{21} - 2 - \tau_{21} = 1, \quad (63)$$

$$\delta^2 \text{ term :} \quad \frac{1-c}{2} \mu_{22} - \frac{s}{2\sqrt{3}} \nu_{22} = 0, \quad (64)$$

$$(1-c)\mu_{22} - \frac{s}{\sqrt{3}} \nu_{22} = 0, \quad (65)$$

$$\nu_{22} - \tau_{22} = 0, \quad (66)$$

$$\nu_{22} - \tau_{22} = 0, \quad (67)$$

$$\delta^3 \text{ term :} \quad \frac{B}{2} + \tau_{22} - \frac{s\kappa}{3\sqrt{3}} = 0, \quad (68)$$

$$\mu_{23} - \frac{B}{2} + \tau_{22} + \frac{s\kappa}{3\sqrt{3}} = \mu_{23}, \quad (69)$$

$$\frac{(1+c)\kappa}{3} - \frac{\sqrt{3} s B}{2(1-c)} + \frac{\tau_{21}^2}{2} + \frac{\tau_{21}^3}{6} - \tau_{23} = \nu_{23}, \quad (70)$$

$$\frac{(1+c)\kappa}{3} + \frac{\sqrt{3} s B}{2(1-c)} + \frac{\tau_{21}^2}{2} + \frac{\tau_{21}^3}{6} - \tau_{23} + \nu_{23} = 0, \quad (71)$$

where $c = \cos 2\sqrt{3}\pi$, $s = \sin 2\sqrt{3}\pi$, and $B = (1-c)\mu_{23} - (s/\sqrt{3})\nu_{23}$. Note that at each stage we have made use of the results from previous stages. For instance, the relations obtained from the equations of the $\delta$ and $\delta^2$ terms have been used to simplify the equations of the $\delta^3$ term. Also note that some of the above 12 equations are dependent (or redundant), for example, equation (61) and equations (64)-(67). The solution to these equations is unique
and is given by

$$ \mu_{21} = 0, \quad \mu_{22} = 0, \quad \mu_{23} = \frac{s\kappa}{3\sqrt{3}(1-c)}, $$

$$ \nu_{21} = 1, \quad \nu_{22} = 0, \quad \nu_{23} = -\frac{\kappa}{3}, $$

$$ \tau_{21} = -2, \quad \tau_{22} = 0, \quad \tau_{23} = \frac{2+\kappa}{3}. $$

Then from equations (54)–(56), we also get

$$ \tau_{11} = 2, \quad \tau_{12} = 0, \quad \tau_{13} = -\frac{2(1+\kappa)}{3}. $$

Consequently, the one parameter family of fixed points on $\Sigma$ can be expressed in terms of $y_1$ as

$$ y_1 = \delta, \tag{72} $$

$$ x_2 = 1 + \frac{s\kappa}{3\sqrt{3}(1-c)} \delta^3 + O(4), \tag{73} $$

$$ y_2 = \delta - \frac{\kappa}{3} \delta^3 + O(4). \tag{74} $$

The corresponding set of fixed points for the second mode, obtained in the same manner, is given by

$$ y_1 = \delta, $$

$$ x_2 = -1 + \frac{s'\kappa}{27\sqrt{3}(1-c')} \delta^3 + O(4), $$

$$ y_2 = -\delta - \frac{\kappa}{27} \delta^3 + O(4), $$

where $s' = \sin(2\pi/\sqrt{3})$ and $c' = \cos(2\pi/\sqrt{3})$.

## 6.2. THE NORMAL MODES

Based on equations (72)–(74), the invariant manifold in the region $x_1 \le 1$ can be constructed. First, let $P$ denote a fixed point on $\Sigma$ with the coordinates described in equations (72)–(74). Since $\omega_1 = 1$, the coordinate transformation in equation (22) becomes

$$ x_1 = r \cos \theta, \quad y_1 = r \sin \theta. \tag{75} $$

Suppose that a point $P_\theta$ on the periodic orbit passing through $P$ has coordinates

$$ x_2 = r \cos \theta + \alpha_{21} \delta + \alpha_{22} \delta^2 + \alpha_{23} \delta^3 + O(4), \tag{76} $$

$$ y_2 = r \sin \theta + \beta_{21} \delta + \beta_{22} \delta^2 + \beta_{23} \delta^3 + O(4), \tag{77} $$

in addition to $x_1$ and $y_1$ given in equation (75). Here $r$ is not arbitrary. Depending on $\delta$ and $\theta$, it must be chosen so that $P_\theta$ lies on the periodic orbit. Hence we take

$$ r = 1 + r_1 \delta + r_2 \delta^2 + r_3 \delta^3 + O(4). \tag{78} $$
The hitting time from $P_\theta$ to $P$ is assumed to be of the form

$$t_{h\theta} = \theta + t_1\delta + t_2\delta^2 + t_3\delta^3 + O(4). \quad (79)$$

Substitution of equations (75)–(79) into equation (29) yields a set of 12 algebraic equations from which one can determine the 12 unknown coefficients in equations (76)–(79). The solution procedure is quite similar to that presented in Section 6.1 and is not repeated here. The resulting solution is

$$\alpha_{21} = 0, \quad \beta_{21} = 0, \quad r_1 = 0, \quad t_1 = -1,$$

$$\alpha_{22} = 0, \quad \beta_{22} = 0, \quad r_2 = \frac{1}{2}, \quad t_2 = 0$$

$$\alpha_{23} = \frac{\kappa}{3\sqrt{3}(1-c)} [s \cos \sqrt{3} \theta + (1-c) \sin \sqrt{3} \theta],$$

$$\beta_{23} = \frac{\kappa}{3(1-c)} [s \sin \sqrt{3} \theta - (1-c) \cos \sqrt{3} \theta],$$

$$r_3 = \frac{\kappa}{12\sqrt{3}(1-c)} \{2s - (\sqrt{3} + 1)[(1-c)\sin(\sqrt{3}-1)\theta + s\cos(\sqrt{3}-1)\theta] \\ + (\sqrt{3}-1)[(1-c)\sin(\sqrt{3}+1)\theta + s\cos(\sqrt{3}+1)\theta]\},$$

$$t_3 = \frac{2+\kappa}{6} - \frac{1}{2}(\alpha_{23} \sin \theta - \beta_{23} \cos \theta).$$

Inversion of the series in equation (78) yields

$$\delta(r, \theta) \approx \sqrt{2(r-1)} - 2r_3(\theta)(r-1).$$

Therefore, we finally have

$$f_{21}(r, \theta) \approx \alpha_{23}(\theta)[\sqrt{2(r-1)} - 2r_3(\theta)(r-1)]^3, \quad (80)$$

$$g_{21}(r, \theta) \approx \beta_{23}(\theta)[\sqrt{2(r-1)} - 2r_3(\theta)(r-1)]^3. \quad (81)$$

For the region $r \cos \theta > 1$, we set

$$x_1 = 1 + \hat{r} \cos \phi, \quad (82)$$

$$y_1 = \hat{r} \sin \phi, \quad (83)$$

$$x_2 = x_1 + \hat{\alpha}_{21}\delta + \hat{\alpha}_{22}\delta^2 + \hat{\alpha}_{23}\delta^3 + O(4), \quad (84)$$

$$y_2 = y_1 + \hat{\beta}_{21}\delta + \hat{\beta}_{22}\delta^2 + \hat{\beta}_{23}\delta^3 + O(4), \quad (85)$$

as the coordinates of a point $P_\phi$ on the periodic trajectory through $P$ (see Figure 5). Again, $\hat{r}$ is not arbitrary. Let $\hat{r}$ and the backward hitting time $t_{h\phi}$ of $P_\phi$ to $P$ be written as

$$\hat{r} = \hat{r}_1\delta + \hat{r}_2\delta^2 + \hat{r}_3\delta^3 + O(4), \quad (86)$$

$$t_{h\phi} = \hat{t}_1\delta + \hat{t}_2\delta^2 + \hat{t}_3\delta^3 + O(4). \quad (87)$$
Solving for the unknown coefficients in equations (84)–(87) from equation (36) (by repeating the solution scheme used above) yields

$$ \hat{\alpha}_{21} = 0, \quad \hat{\alpha}_{22} = 0, \quad \hat{\alpha}_{23} = \frac{s\kappa}{3\sqrt{3}(1-c)}, $$

$$ \hat{\beta}_{21} = 0, \quad \hat{\beta}_{22} = 0, \quad \hat{\beta}_{23} = 0, $$

$$ \hat{r}_1 = 0, \quad \hat{r}_2 = \frac{1}{2\cos\phi}, \quad \hat{r}_3 = 0, $$

$$ \hat{t}_1 = 1, \quad \hat{t}_2 = -0.5 \tan \phi, \quad t_3 = -\frac{1+\kappa}{3}. $$

By equation (86) and transforming $(\hat{r}, \phi)$ back into $(r, \theta)$, we have

$$ \delta(r, \theta) \approx \sqrt{2(r \cos \theta - 1)}. $$

Thus, we finally have

$$ f_{22}(r, \theta) \approx \frac{2\sqrt{2} s\kappa}{3\sqrt{3}(1-c)} (r \cos \theta - 1)^{3/2}, \quad (88) $$

$$ g_{22}(r, \theta) \approx 0. \quad (89) $$

The same procedure can be applied to the second mode, resulting in

$$ f'_{21}(r, \theta) \approx \alpha'_{23}(\theta)[\sqrt{6(r-1)} - 18r'_{3}(\theta)(r-1)]^3, $$

$$ g'_{21}(r, \theta) \approx \beta'_{23}(\theta)[\sqrt{6(r-1)} - 18r'_{3}(\theta)(r-1)]^3, $$

$$ f'_{22}(r, \theta) \approx \frac{2\sqrt{6} s'\kappa}{9(1-c')} (r \cos \theta - 1)^{3/2}, $$

$$ g'_{22}(r, \theta) \approx 0, $$

where

$$ \alpha'_{23} = -\frac{\kappa}{27(1-c')} \left[ s' \cos \frac{\theta}{\sqrt{3}} + (1-c') \sin \frac{\theta}{\sqrt{3}} \right], $$

$$ \beta'_{23} = \frac{\kappa}{27(1-c')} \left[ s' \sin \frac{\theta}{\sqrt{3}} - (1-c') \cos \frac{\theta}{\sqrt{3}} \right], $$

$$ r'_{3} = \frac{\kappa}{108\sqrt{3}(1-c')} \left\{ 2\sqrt{3}s - (\sqrt{3}+1)\left[(1-c')\sin\frac{\sqrt{3}-1}{\sqrt{3}}\theta + s'\cos\frac{\sqrt{3}-1}{\sqrt{3}}\theta\right] - (\sqrt{3}-1)\left[(1-c')\sin\frac{\sqrt{3}+1}{\sqrt{3}}\theta + s'\cos\frac{\sqrt{3}+1}{\sqrt{3}}\theta\right] \right\}. $$

Recall that the parameter $\kappa$, which is the stiffness of the free spring, accounts for the only nonlinearity in our system. In fact, it can be verified that all of the functions $f_{ij}$, $g_{ij}$, $f'_{ij}$, and $g'_{ij}$ are identically zero when $\kappa = 0$. In this case, all of the modal subspaces are simply two-dimensional planes.
<figure><img src="image_8.png" /><figcaption>Figure 7. The modal subspace of the first mode with $κ = 3$.</figcaption></figure>

Using equations (80)–(81) and (88)–(89), the modal subspaces of the first mode are plotted in Figure 7. Similarly, the modal subspaces of the second mode are shown in Figure 8. From the figures, one can see that the normal modes are of piecewise type with a flat disk at the center.

### 6.3. MODAL DYNAMICS

It follows directly from equations (43) and (44) that the dynamics of the first mode are

$$ \dot{r} = \begin{cases} 0 & \text{for } r \le 1 \\ \alpha_{23}(\theta) \sin \theta [\sqrt{2(r-1)} - 2r_3(\theta)(r-1)]^3 & \text{for } r > 1 \text{ and } r \cos \theta \le 1 \\ \kappa \sin \theta \left[(r \cos \theta - 1) + \frac{2\sqrt{2}s}{3\sqrt{3}(1-c)}(r \cos \theta - 1)\right]^{3/2} & \text{for } r \cos \theta > 1 \end{cases} \quad (90) $$

$$ \dot{\theta} = \begin{cases} -1 & \text{for } r \le 1 \\ -1 + \frac{\alpha_{23}(\theta) \cos \theta}{r} [\sqrt{2(r-1)} - 2r_3(\theta)(r-1)]^3 & \text{for } r > 1 \text{ and } r \cos \theta \le 1 \\ -1 + \frac{\kappa \cos \theta}{r} [(r \cos \theta - 1) \\ + \frac{2\sqrt{2}s}{3\sqrt{3}(1-c)}(r \cos \theta - 1)^{3/2}] & \text{for } r \cos \theta > 1 \end{cases} \quad (91) $$

while those for the second mode are

$$ \dot{r} = \begin{cases} 0 & \text{for } r \le 1 \\ \frac{\alpha'_{23}(\theta) \sin \theta}{\sqrt{3}} [\sqrt{6(r-1)} - 18r'_3(\theta)(r-1)]^3 & \text{for } r > 1 \text{ and } r \cos \theta \le 1 \\ \frac{\kappa \sin \theta}{\sqrt{3}} [(r \cos \theta - 1) - \frac{2\sqrt{6}s'}{9(1-c')} (r \cos \theta - 1)^{3/2}] & \text{for } r \cos \theta > 1 \end{cases} $$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. The modal subspace of the second mode with $\kappa = 3$.</figcaption></figure>

$$ \dot{\theta} = \begin{cases} -\sqrt{3} & \text{for } r \le 1 \\ -\sqrt{3} + \frac{\alpha'_{23}(\theta) \cos \theta}{\sqrt{3} r} [((\sqrt{6(r-1)} - 18r'_{3}(\theta)(r-1)})^3 & \text{for } r > 1 \text{ and } r \cos \theta \le 1 \\ -\sqrt{3} + \frac{\kappa \cos \theta}{\sqrt{3} r} [((r \cos \theta - 1) + \frac{2\sqrt{6} s'}{9(1-c')} (r \cos \theta - 1)^{3/2})] & \text{for } r \cos \theta > 1 \end{cases} $$

The frequency-amplitude relationships are obtained from equation (45). For the first mode in the region $r > 1$, it is

$$ \tilde{\omega}_1 = 1 + \frac{\kappa}{6\pi} \delta^3 + O(4). $$

The parameter $\delta (= y_1)$ is related to the peak displacement $A_1$ of $x_1$ via equation (46) as follows

$$ A_1 = 1 + \frac{1}{\kappa} \delta^2 + O(4). $$
Hence we have

$$ \delta(A_1) \approx \sqrt{2(A_1 - 1)} $$

and

$$ \tilde{\omega}_1 \approx \begin{cases} 1 & \text{for } A_1 \le 1 \\ 1 + \frac{\sqrt{2} \kappa}{3\pi} (A_1 - 1)^{3/2} & \text{for } A_1 > 1 \end{cases} $$

The corresponding relationship for the second mode is

$$ \tilde{\omega}_2 \approx \begin{cases} \sqrt{3} & \text{for } A_1 \le 1 \\ \sqrt{3} + \frac{\sqrt{6} \kappa}{9\pi} (A_1 - 1)^{3/2} & \text{for } A_1 > 1 \end{cases} $$

Again, observe that as $\kappa = 0$, the frequencies will be constant, corresponding to the limiting linear system.

It is also interesting to note that the modal motion for this system is synchronous in the sense that both masses simultaneously reach their peaks. However, the displacements do not necessarily vanish simultaneously. Such a phenomenon has been verified in [41] through simulations for a special case of PWL systems. To see this, observe that $\beta_{23}(\pi)$ and $\beta'_{23}(\pi)$ are zeros, which yield $g_{21}(r, \theta) \approx 0$ and $g'_{21}(r, \theta) \approx 0$ when $\theta = \pi$. Moreover, we already have $g_{22}(r, \theta) \approx 0$ and $g'_{22}(r, \theta) \approx 0$. These imply that both $y_1$ and $y_2$ are necessarily zero when $\sin \theta = 0$ by equation (24), the equation for the normal mode. Hence, both modes are synchronous as described above. Thus, we can compute the ratio of peak amplitudes $A_2/A_1$ for each mode, which are:

$$ \text{Mode 1 :} \quad \frac{A_2}{A_1} \approx \begin{cases} 1 & \text{for } A_1 \le 1 \\ 1 - \frac{\sqrt{3} (1-c_1) \kappa}{3\sqrt{3} s_1} (A_1 - 1)^{3/2} & \text{for } A_1 > 1 \end{cases} $$

$$ \text{Mode 2 :} \quad \frac{A_2}{A_1} \approx \begin{cases} 1 & \text{for } A_1 \le 1 \\ 1 + \frac{2\sqrt{6} s\kappa}{9(1-c)} (A_1 - 1)^{3/2} & \text{for } A_1 > 1 \end{cases} $$

where $c_1 = \cos \sqrt{3} \pi$ and $s_1 = \sin \sqrt{3} \pi$.

## 6.4. SIMULATION RESULTS AND COMPARISON WITH THE RITZ METHOD

In order to examine the accuracy and validity of the analysis, the exact modal solutions are needed. For a given mode and a given $\delta$, a fixed point on the Poincaré section $\Sigma$ is obtained by solving the fixed point equations (15) and (16) numerically using the fixed point approximation (72)–(74) as an initial guess. With this fixed point as an initial condition, one can then obtain the modal solution by integrating the equations of motion, equations (48) and (49).

Although conventional perturbation techniques fail to apply to PWL systems as mentioned previously, some energy-based methods, like the Ritz method, are applicable here. As a comparison, the Ritz method will be employed to the example PWL system (48) and (49). Since the Ritz method can be found in many books (for example, [33]), only final results are presented here. The detailed calculations are shown in Appendix B.

In what follows, we shall compare the present method with the Ritz method and the exact solutions obtained through numerical simulations. To guarantee that they belong to the same
<figure><img src="image_11.png" /><figcaption>*Figure 9. Relationship of modal frequencies and $\delta$ for different methods.*</figcaption></figure>

energy level, the initial conditions for different approaches are chosen in the following way. For a given mode and a fixed $\delta$, the fixed point approximation (equations (72)–(74)) is taken to be the initial condition for the present approach. The exact solution uses the exact fixed point under the same mode and $\delta$. From this exact solution, one can get the maximum displacement of the first mass $A_1$. Then, the initial condition for the Ritz method is given by $[A_1 \ 0 \ A_1b_2 \ 0]$, where $b_2$ can be found in Appendix B.

In simulations it is observed that for small $\delta$, both the present method and the Ritz method have good estimates to the modal solution. However, as $\delta$ increases, the Ritz method gives a better approximation than the present method for the periodic solutions. This is because the present method is based on asymptotic expansion which holds only for small $\delta$, i.e., in the neighborhood of the invariant disk.

The frequency and peak amplitudes relationships are given in Figures 9 and 10 for different combinations of $\kappa$ and modes. It is easy to see that good approximations are achieved for both methods for small $\delta$ ($\delta \le 0.4$), while for large $\delta$, the Ritz method is more accurate than the present method. Also, the effect of $\kappa$ can be clearly seen.

# 7. Conclusions

Based on the construction of a convenient Poincaré map and making use of invariant manifold theory, a general procedure for constructing the normal modes of a class of piecewise linear (PWL) vibratory systems has been developed. The class of systems considered contains two linear subsystems separated by a switching hyperplane in state space. It is shown that such a PWL system possesses a canonical form upon which the proposed method applies.
<figure><img src="image_12.png" /><figcaption>Figure 10. Ratio of peak amplitudes for different methods.</figcaption></figure>

There are several steps in this constructive procedure. First, using the switching hyperplane as a Poincaré section $Σ$, a one parameter family of fixed points on $Σ$ is obtained for each mode. Then, using a form of polar coordinates, the state space is decomposed into three regions. The first one ($r \le d$) is the region of small amplitude vibrations in which the usual linear modes hold. While both the second ($r > d$ and $r \cos\theta \le d$) and third ($r \cos\theta > d$) regions correspond to large amplitude motions, they are governed by two different linear subsystems. In the first region, the normal modes are invariant elliptic disks and are easily obtained. The normal modes in the other two regions are constructed by matching the fixed point family on the Poincaré section. Each mode in the second region is expanded in a series form around a neighborhood of the invariant disk in the first region. In the last region, the normal mode is expanded in a series form near the point where the invariant disk in the first region and the switching hyperplane are tangent.

From the expressions for the normal modes, the modal dynamics and corresponding frequency-amplitude dependence relationship on each mode are obtained. As expected, the normal modes, the modal dynamics, and the frequency-amplitude relationships are all of piecewise type (although not piecewise linear).

An undamped, two degrees of freedom system with a clearance is demonstrated as an example to illustrate the procedure. The numerical results demonstrate the effectiveness of the procedure. It must be admitted that the Ritz method is far easier to apply and is more accurate in terms of estimating frequencies. The advantage of the present method is that it allows one to construct the differential equations that give the modal dynamics.
# Appendix A. A Necessary Condition for the Existence of Normal Modes for PWL Systems

The fixed point equations (15) and (16) are equivalent to the following equation:

$$ \begin{pmatrix} w \\ w^* \end{pmatrix} = \begin{pmatrix} 0 \\ w_e \end{pmatrix} + \begin{bmatrix} 0 & T(t_{h2}) \\ \hat{T}(t_{h1}) & 0 \end{bmatrix} \begin{pmatrix} w - w_e \\ w^* \end{pmatrix}, $$

where $w^* = [x^* \ y^*]^T \in \Sigma^*$. This can be written as

$$ G(\xi, \kappa) = 0, \qquad (92) $$

where $\xi = (t_{h1}, x_2, \dots, x_N, y_1, \dots, y_N; t_{h2}, x_2^*, \dots, x_N^*, y_1^*, \dots, y_N^*)$.

We know that when $\kappa = 0$, equation (92) possesses the following solution for the first mode:

$$ \xi_0 = \left( 0, u_2, \dots, u_N, 0, \dots, 0; \frac{2\pi}{\omega_1}, u_2, \dots, u_N, 0, \dots, 0 \right), $$

which is the linear solution. The Jacobian matrix of equation (92) at $(\xi_0, 0)$ is given by (after some manipulations)

$$ J = \left. \frac{\partial G}{\partial \xi} \right|_{\substack{\xi=\xi_0 \\ \kappa=0}} = \begin{bmatrix} J_1 & T(\frac{2\pi}{\omega_1})J_2 \\ J_2 & J_1 \end{bmatrix}, $$

where the $2N \times 2N$ matrices $J_1$ and $J_2$ are defined by

$$ J_1 = \begin{bmatrix} 0 & \cdots & 0 \\ \vdots \\ 0 & I_{2N-1} \end{bmatrix}, $$

$$ J_2 = \begin{bmatrix} 0_N & 0 \cdots 0 \\ \alpha v_1 & I_{2N-1} \end{bmatrix}, $$

where $\alpha = d\omega_1^2/v_{11}$, $0_N$ stands for the $N$-dimensional zero vector, and $I_m$ is the identity matrix of order $m$. If $J$ is nonsingular, by the implicit function theorem, there exists only one fixed point for each $\kappa$ near 0. Then it is impossible to have a one-parameter family of fixed passing through $(\xi_0, 0)$. Thus, $J$ must be singular in order for the normal mode to exist. It will be shown below that $J$ is indeed not of full rank and exactly has rank equal to $4N - 1$ if the system has no internal resonance.

Let $j_i$ be the $i$-th column of $J$. First, replace $j_1$ and $j_{2N+1}$ in $J$ by

$$ j_1' = \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \\ -1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \quad \text{and} \quad j_{2N+1}' = \begin{pmatrix} T(\frac{2\pi}{\omega_1}) & \begin{pmatrix} -1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} \\ 1 & 0 \\ \vdots & \vdots \\ 0 & 0 \end{pmatrix}, $$
to form a new $4N \times 4N$ matrix

$$J' = \begin{bmatrix} I_{2N} & -T\left(\frac{2\pi}{\omega_1}\right) \\ -I_{2N} & I_{2N} \end{bmatrix}.$$

Let $(\begin{smallmatrix} u \\ v \end{smallmatrix})$ be a null vector of $J'$. Then

$$\left. \begin{array}{l} u - T\left(\frac{2\pi}{\omega_1}\right)v = 0 \\ -u + v = 0 \end{array} \right\} \implies \left( I_{2N} - T\left(\frac{2\pi}{\omega_1}\right) \right) u = 0,$$

where we compute

$$\begin{align*} I_{2N} - T\left(\frac{2\pi}{\omega_1}\right) &= I_{2N} - \begin{bmatrix} V & 0 \\ 0 & V \end{bmatrix} \begin{bmatrix} C\left(\frac{2\pi}{\omega_1}\right) & S\left(\frac{2\pi}{\omega_1}\right)\Omega^{-1} \\ -\Omega S\left(\frac{2\pi}{\omega_1}\right) & C\left(\frac{2\pi}{\omega_1}\right) \end{bmatrix} \begin{bmatrix} V^T & 0 \\ 0 & V^T \end{bmatrix} \\ &= \begin{bmatrix} V & 0 \\ 0 & V \end{bmatrix} \begin{bmatrix} I_N - C\left(\frac{2\pi}{\omega_1}\right) & -S\left(\frac{2\pi}{\omega_1}\right)\Omega^{-1} \\ \Omega S\left(\frac{2\pi}{\omega_1}\right) & I_N - C\left(\frac{2\pi}{\omega_1}\right) \end{bmatrix} \begin{bmatrix} V^T & 0 \\ 0 & V^T \end{bmatrix}, \end{align*}$$

where

$$I_N - C\left(\frac{2\pi}{\omega_1}\right) = \operatorname{diag}\left(0, 1 - \cos\frac{\omega_2}{\omega_1} 2\pi, \dots, 1 - \cos\frac{\omega_N}{\omega_1} 2\pi\right),$$

$$S\left(\frac{2\pi}{\omega_1}\right) = \operatorname{diag}\left(0, \sin\frac{\omega_2}{\omega_1} 2\pi, \dots, \sin\frac{\omega_N}{\omega_1} 2\pi\right).$$

Therefore, without internal resonance, we will have

$$1 - \cos \frac{\omega_i}{\omega_1} 2\pi \neq 0 \quad \text{and} \quad \sin \frac{\omega_i}{\omega_1} 2\pi \neq 0$$

for each $i = 2, \dots, N$. It follows that the null space dimension of $I_{2N} - T(2\pi/\omega_1)$ and therefore of $J'$ is 2.

Next, pull out $j'_1$ and $j'_{2N+1}$ from $J'$. The remaining matrix becomes

$$J'' = \begin{bmatrix} 0 & T\left(\frac{2\pi}{\omega_1}\right) \begin{pmatrix} 0 \\ I_{2N-1} \end{pmatrix} \\ I_{2N-1} & 0 \\ 0 & I_{2N-1} \end{bmatrix},$$

which has $4N - 2$ columns. It is obvious that $j'_1$ is independent of the columns of $J''$. Also, it is not difficult to check that $j'_{2N+1}$ is dependent on $j'_1$ and the columns of $J''$. Thus, we deduce that the rank of $J''$ is $4N - 3$.

Finally, one can prove that adding $j'_1$ and $j'_{2N+1}$ back to $J''$ will increase the matrix rank by 2, which results in our claim that $J$ has a rank equal to $4N - 1$.

If $\omega_i/\omega_1$ is an integer for some $i \neq 1$, then $I_{2N} - T(2\pi/\omega_1)$ will have a null space of dimension more than 2, which implies the rank of $J$ will be less than $4N - 1$.
# Appendix B. Detailed Calculations for the Ritz Method

The Ritz method assumes a modal solution of the form:

$$x_1(t) = A_1 \cos \omega t,$$

$$x_2(t) = A_1 b_2 \cos \omega t,$$

where $\omega$ is the modal frequency and $b_2$ is the ratio of peak amplitudes to be determined. The system equation is rewritten as:

$$\epsilon_i \equiv \ddot{x}_i + f_i(x_1, x_2) = 0, \quad i = 1, 2.$$

Then by minimizing the solution error in the sense

$$\int_{0}^{2\pi} \epsilon_i(\theta) \, d\theta = 0, \quad i = 1, 2, \qquad (93)$$

where $\theta = \omega t$, $\omega$ and $b_2$ can be obtained in terms of $A_1$. For the example system given in equations (48) and (49), equation (93) reads

$$-\omega^2 + 2 - b_2 + g(A_1) = 0,$$

$$-\omega^2 b_2 - 1 + 2b_2 = 0,$$

where

$$g(A_1) = \begin{cases} 0 & \text{for } A_1 \le 1 \\ \frac{\kappa}{2\pi} \left[ \cos^{-1} \frac{1}{A_1} - \frac{\sqrt{A_1^2 - 1}}{A_1^2} \right] & \text{for } A_1 > 1 \end{cases}.$$

Thus we have

$$b_2 = g(A_1) \pm \sqrt{g(A_1)^2 + 1},$$

$$\omega = \sqrt{2 - \frac{1}{b_2}},$$

where + corresponds to the first mode and – to the second mode.

## References

1. Atkinson, C. and Taskett, B., 'A study of the nonlinearly related modal solutions of coupled nonlinear systems by superposition techniques', *Journal of Applied Mechanics* **32**, 1965, 359–364.

2. Benamar, R., Bennouna, M., and White, R., 'The effects of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures', *Journal of Sound and Vibration* **149**, 1991, 179–195.

3. Benamar, R., Bennouna, M., and White, R., 'The effects of large vibration amplitudes on the fundamental mode shape of a clamped-clamped uniform beam', *Journal of Sound and Vibration* **96**, 1984, 309–331.

4. Chua, L. O., Komuro, M., and Matsumoto, T., 'The double scroll family', *IEEE Transactions on Circuits and Systems* **33**, 1986, 1073–1118.

5. Cooke, C. and Struble, R., 'The existence of periodic solutions and normal mode vibrations in nonlinear systems', *Quarterly of Applied Mathematics* **24**, 1966, 177–193.

6. Caughey, T. K., 'Classical normal modes in damped linear dynamics systems', *Journal of Applied Mechanics* **32**, 1965, 359–364.
7. Caughey, T. K. and Vakakis, A. F., 'A method for examining steady state solutions of forced discrete systems with strong nonlinearities', *International Journal of Non-Linear Mechanics* **26**, 1991, 89-103.

8. Greenberg, H. and Yang, T.-L., 'Modal subspaces and normal mode vibrations', *International Journal of Non-Linear Mechanics* **6**, 1971, 311-326.

9. Kim, Y. B. and Noah, S. T., 'Stability and bifurcation analysis of oscillators with piecewise linear characteristics: A general approach', *Journal of Applied Mechanics* **58**, 1991, 545-553.

10. Month, L. A. and Rand, R. H., 'An application of Poincaré map to the stability of nonlinear normal modes', *Journal of Applied Mechanics* **47**, 1980, 645-651.

11. Natsiavas, S., 'On the dynamics of oscillators with bilinear damping and stiffness', *International Journal of Non-Linear Mechanics* **25**, 1990, 311-326.

12. Nayfeh, A. and Nayfeh, S., 'On nonlinear modes of continuous systems', *Journal of Vibration and Acoustics* **116**(1), 1994, 129-136.

13. Nayfeh, A. and Nayfeh, S., 'Nonlinear normal modes of a continuous system with quadratic nonlinearities', *Journal of Vibration and Acoustics* **117**, 1995, 199-205.

14. Nayfeh, A. and Chin, C., 'On nonlinear modes of systems with internal resonance', *Journal of Vibration and Acoustics* **118**, 1996 (to appear).

15. Pak, C. and Rosenberg, R., 'On the existence of normal mode vibrations in nonlinear systems', *Quarterly of Applied Mathematics* **26**, 1968, 403-416.

16. Pak, C., Rand, R., and Moon, F., 'Free vibrations of a thin elastica by normal modes', *Nonlinear Dynamics* **3**, 1992, 347-364.

17. Rand, R., Pak, C., and Vakakis, A., 'Bifurcation of nonlinear normal modes in a class of 2 degrees of freedom systems', *Acta Mechanica* **3**, 1992, 129-145.

18. Rand, R., 'Nonlinear normal mode in 2 degrees of freedom systems', *Journal of Applied Mechanics* **38**, 1971, 561.

19. Rand, R., 'A direct method for nonlinear normal modes', *International Journal of Non-Linear Mechanics* **9**, 1974, 363-368.

20. Rand, R., 'A higher order approximation for nonlinear normal modes in two degrees of freedom systems', *International Journal of Non-Linear Mechanics* **6**, 1971, 545-547.

21. Rand, R., 'The geometrical stability of nonlinear normal modes in 2 degrees of freedom systems', *International Journal of Non-Linear Mechanics* **8**, 1973, 161-168.

22. Rosenberg, R. M., 'The normal modes of nonlinear n degrees of freedom systems', *Journal of Applied Mechanics* **30**, 1962, 7-14.

23. Rosenberg, R. M., 'On nonlinear vibrations of systems with many degrees of freedom', *Advance in Applied Mechanics* **9**, 1966, 155-242.

24. Rosenberg, R. M. and Kuo, J. K., 'Nonsimilar normal mode vibrations of nonlinear systems having two degrees of freedom', *Journal of Applied Mechanics* **31**, 1964, 283-290.

25. Shaw, S. W., 'On the dynamic response of a system with dry friction', *Journal of Sound and Vibration* **108**, 1986, 305-325.

26. Shaw, S. W., 'An invariant manifold approach to nonlinear normal modes of oscillation', *Journal of Nonlinear Science* **4**, 1994, 419-448.

27. Shaw, S. W. and Holmes, P. J., 'A periodically forced piecewise linear oscillator', *Journal of Sound and Vibration* **90**, 1983, 305-325.

28. Shaw, S. W. and Pierre, C., 'Normal modes for nonlinear vibratory systems', *Journal of Sound and Vibration* **164**, 1993, 85-124.

29. Shaw, S. W. and Pierre, C., 'Nonlinear normal modes and invariant manifolds', *Journal of Sound and Vibration* **150**, 1991, 170-173.

30. Shaw, S. W. and Pierre, C., 'Normal modes of vibration for nonlinear continuous systems', *Journal of Sound and Vibration* **169**, 1994, 319-347.

31. Sparrow, C.T., 'Chaos in a three-dimensional single loop feedback system with a piecewise linear feedback function', *Journal of Mathematical Analysis and Applications* **83**, 1981, 275-291.

32. Szemplinska-Stupnicka, W., 'The resonant vibration of homogeneous nonlinear systems', *International Journal of Non-Linear Mechanics* **15**, 1980, 407-415.

33. Szemplinska, W., *The Behavior of Nonlinear Vibrating Systems*, Kluwer, Dordrecht, 1990.

34. Thompson, J. M. T., Bokaian, A. R., and Ghaffari, R., 'Subharmonic resonance and chaotic motions of a bilinear oscillator', *IMA Journal of Applied Mathematics* **37**, 1983, 207-234.

35. Vakakis, A., 'Analysis and identification of linear and nonlinear normal modes in vibrating systems', Ph.D. Dissertation, California Institute of Technology, 1990.

36. Vakakis, A. F. and Caughey, T. K., 'A theorem on the exact nonsimilar steady state motions of a nonlinear oscillator', *Journal of Applied Mechanics* **5D**, 1992, 418-424.
37. Vakakis, A. and Rand, R., 'Normal modes and global dynamics of a 2 degrees of freedom system – I. Low energies', *International Journal of Non-Linear Mechanics* **27**, 1992, 861–874.

38. Vakakis, A. and Rand, R., 'Normal modes and global dynamics of a 2 degrees of freedom system – II. High energies', *International Journal of Non-Linear Mechanics* **27**, 1992, 875–888.

39. Van der Varst, P., 'On normal mode vibrations of nonlinear conservative systems', Doctoral Thesis, Germany, 1982.

40. Yen, D., 'On the normal modes of nonlinear dual-mass systems', *International Journal of Non-Linear Mechanics* **9**, 1974, 45–53.

41. Zuo, L. and Curnier, A., 'Nonlinear real and complex modes of conewise linear systems', *Journal of Sound and Vibration* **174**, 1994, 289–313.