# Non-Linear Real And Complex Modes Of Conewise Linear Systems

L Zuo, Alain Curnier

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

## To cite this version:

L Zuo, Alain Curnier. Non-Linear Real And Complex Modes Of Conewise Linear Systems. Journal of Sound and Vibration, 1994, 174 (3), pp.289-313. 10.1006/jsvi.1994.1277. hal-01338607

HAL Id: hal-01338607

https://hal.science/hal-01338607

Submitted on 28 Jun 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NON-LINEAR REAL AND COMPLEX MODES OF
CONEWISE LINEAR SYSTEMS

L. ZUO AND A. CURNIER

*Applied Mechanics Laboratory, Department of Mechanics, Swiss Federal Institute of
Technology at Lausanne, CH-1015 Lausanne, Switzerland*

Conewise (piecewise) linear, non-gyroscopic and gyroscopic autonomous systems are studied in detail. First, the general properties of such systems are examined. Then the non-linear real and complex modes are obtained for both non-gyroscopic and gyroscopic cases, by means of analytical or numerical tools. Several interesting characteristics of non-linear modes are found and compared to those of linear modes. Accurate analytical estimates of non-linear mode periods are formulated for both non-gyroscopic and gyroscopic systems. The validity of such estimates is discussed in relation to a non-linear mode bifurcation phenomenon, which is located by using both characteristic multipliers and Poincaré mapping.

# 1. INTRODUCTION

Actual engineering systems often involve components which are in unilateral contact, operate with clearances, are made of bimodular elastic materials, etc. Such systems are governed by equations of motion containing piecewise linear terms. The example which has motivated this study is that of rotating shafts with a transverse crack. Indeed, when a cracked shaft revolves, the crack alternatively opens and closes and the shaft stiffness can be roughly modelled by a half-spacewise constant function [1-5]. Natural frequencies of rotating shafts are important in system design in order to avoid contingent resonances [6]. Crack detection is even more crucial for obvious safety reasons. The reduction of the "natural frequencies" with the crack depth could serve as an indicator of crack initiation and propagation [1, 7, 8].

Owing to its practical and theoretical importance, the dynamic behavior of piecewise
linear (PL) systems has been the subject of many investigations. The free responses of
single-degree-of-freedom (1-DOF) PL oscillators were mentioned by Timoshenko and
Young in their book on vibrations [9]. In subsequent studies, most authors considered
periodical excitations and used classical analytical tools for their investigations. When
dealing with a beam in intermittent contact with a spring, Masri *et al.* [10] sought an
analytical solution separately in the contact and gap regions and then matched these
solutions through a condition of continuity at the switching point. A similar approach was
followed by Natsiavas [11, 12], who examined some particular oscillators with bilinear and
trilinear restoring forces, and studied the existence and stability of their motions. For other
PL models, Maezawa *et al.* [13-15] used Fourier series for determining the superharmonic
resonances up to second order; Comparin and Singh [16] employed the harmonic balance
method to study the main resonance; Choi and Noah [17] combined the fast Fourier
transform (FFT) with a Newton-Raphson scheme for approximating the arbitrarily forced periodic vibrations. Recently, Kim and Noah [18] have modified this method in order to enhance its reliability and efficiency.

During the past ten years, further studies on the dynamics of forced 1-DOF PL systems have been carried out with the help of more elaborate non-linear analysis tools: Poincaré section, averaging method, local bifurcation theory, center manifolds, normal forms, etc. [19, 20]. It has been found that, in certain circumstances, the periodic responses of such systems under excitations may go through a bifurcation and exhibit complicated and unpredictable behavior known as chaos. A detailed analysis of a periodically forced PL oscillator was presented by Shaw and Holmes [21] in 1983, while extensive experimental and numerical results on the characteristics of such systems, in particular the period-doubling bifurcation, were supplied by Thompson et al. [22-25] and Mahfouz and Badrakhan [26]. Similar observations in either mechanical or electrical systems were also reported by Schulman [27], Natsiavas [11], Chua and Hasler [28], Chua and Matsumoto [29, 30], and Kim and Noah [18]. In addition, some methods for eliminating chaos in PL systems were proposed by Elvey [31], Thompson and Elvey [32] and Shaw [33].

Up to now, few results on the dynamics of multi-degree-of-freedom (MDOF) PL systems have been published. Recently, the frequency response characteristics of a MDOF system with clearances were investigated by Comparin and Singh [34] with the help of the harmonic balance method. Although their study provides a general method for the development of simplified approximate solutions, it is limited to harmonically forced systems.

In this work, we study autonomous conewise linear multi-degree-of-freedom gyroscopic and non-gyroscopic systems which, to our knowledge, have not been treated before. Such systems are of both theoretical and practical importance because they are the simplest models of cracked beams (non-gyroscopic) and cracked rotating shafts (gyroscopic) as well as of many other mechanical systems subject to unilateral contacts [1-5].

Our approach consists of combining analytical methods with numerical ones for investigating the fundamental dynamic characteristics of such systems; namely, their non-linear real and complex modes. In the light of the analytical results obtained for a 1-DOF model and a particular 2-DOF model, the concept of non-linear normal mode introduced by Rosenberg for smooth and symmetric non-linear systems [35-39] is extended to conewise linear systems. An efficient numerical algorithm is developed to extract all the non-linear modes of gyroscopic and non-gyroscopic systems and compare them with linear modes. Then, an analytical estimate of the natural frequencies of conewise linear systems is established in terms of the system parameters. This estimate not only provides additional insight into system behavior, but also a check for the numerical results obtained. Finally, the stability of the non-linear modes is discussed by using modern analysis tools: characteristics multipliers and Poincaré mapping. The bifurcation of 2-DOF non-linear modes is revealed.

## 2. CONEWISE LINEAR SYSTEMS

### 2.1. DEFINITION

The dynamic equation of the n-DOF piecewise linear gyroscopic systems considered here can be stated, together with its initial conditions, as follows. Find $x(t): [0, T] \to \mathbb{R}^n$ such that

$$ M\ddot{x} + G\dot{x} + f(x) = 0, \quad x(0) = x_0, \quad \dot{x}(0) = v_0. \quad (1) $$
In this equation, $M = M^T$ is a symmetric positive definite matrix representing the system mass, $G = -G^T$ is a skew matrix describing the gyroscopic effect and $f(x)$ is a half spacewise linear elastic force: i.e., a force defined, by means of two constant symmetric matrices $S_+ = S_+^T$ and $S_- = S_-^T$ over two complementary half-spaces $\mathbb{R}_+^n$ and $\mathbb{R}_-^n$ respectively, as

$$f(x) = S(x)x = \begin{cases} S_+ x, & \forall x \in \mathbb{R}_+^n \equiv \{x | x^T n > 0\} \\ S_- x, & \forall x \in \mathbb{R}_-^n \equiv \{x | x^T n \le 0\} \end{cases}. \quad (2)$$

Thus, the two regions of linearity are half-spaces separated by a hyperplane $x^T n = 0$, oriented by its unit normal vector $n$ pointing inward $\mathbb{R}_+^n$. Strictly speaking, models defined by equations (1) and (2) should hence be called *half-spacewise linear* systems. In such systems, the non-linearity is concentrated at the interface which includes the origin. It occurs when two system components are in unilateral contact in the rest position and one of the components has a negligible mass in comparison to the other. Systems involving several such unilateral contacts can give rise to more complicated subdivisions of $\mathbb{R}^n$ into polyhedral convex cones centered at the origin [40], by means of several hyperplanes. These systems would be more appropriately called *conewise linear* than piecewise linear. All the results presented in this paper on half-spacewise linear systems for didactic reasons can be directly extended to conewise linear systems.

A simple 2-DOF “gyroscopic” example is shown in Figure 1. It consists of a mass point attached to a rotating table by two springs and in unilateral contact with a third (massless) one. For this example, equations (1) and (2) are more explicitly defined by

$$M = m I_2, \quad G = \begin{pmatrix} 0 & -2m\Omega \\ 2m\Omega & 0 \end{pmatrix}, \quad n = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \\ S(x) = K(x) - m\Omega^2 I_2, \quad K(x) = \begin{cases} K_+ = K_+^T, & x \in \mathbb{R}_+^2 \equiv \{x | x_1 > 0\} \\ K_- = K_-^T, & x \in \mathbb{R}_-^2 \equiv \{x | x_1 \le 0\} \end{cases}. \quad (3)$$

In equations (3), $K_+$ and $K_-$ represent the stiffnesses of the springs with and without contact, respectively, and $I_2$ is the two-dimensional identity matrix. It should be emphasized that the symmetric matrix $S(x)$ is not necessarily positive definite, due to the gyroscopic term $-m\Omega^2 I_2$.

In phase space, equations (1) and (2) respectively become

$$\dot{y} = h(y) = H(y)y = 0, \quad y(0) = y_0, \quad (4)$$

where

$$H(y) = \begin{cases} \begin{bmatrix} 0 & I_n \\ -M^{-1}S_+ & -M^{-1}G \end{bmatrix} \equiv H_+, & \forall y \in \mathbb{R}_+^{2n} \equiv \{y | y^T m > 0\} \\ \begin{bmatrix} 0 & I_n \\ -M^{-1}S_- & -M^{-1}G \end{bmatrix} \equiv H_-, & \forall y \in \mathbb{R}_-^{2n} \equiv \{y | y^T m \le 0\} \end{cases}, \quad (5)$$

<figure><img src="image_3.png" /><figcaption>Figure 1. The simplest gyroscopic model. X, fixed reference; x, rotating reference.</figcaption></figure>
and

$$y = \begin{pmatrix} x \\ \dot{x} \end{pmatrix} \in \mathbb{R}^{2n}, \quad I_n \equiv \begin{bmatrix} 1 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & 1 \end{bmatrix} \in \mathbb{R}^{n \times n}, \quad \mathbf{m} = \begin{pmatrix} \mathbf{n} \\ 0 \end{pmatrix} \in \mathbb{R}^{2n}.$$

## 2.2. PROPERTIES

The main properties of MDOF half-spacewise linear systems are now deduced from their definition, as a preliminary analysis. All of these properties remain valid for conewise linear systems with minor adjustments.

### 2.2.1. Continuity of $f(x)$: existence and uniqueness of the solution $x(t)$

For mechanical and mathematical consistency, the internal force $f(x)$ is assumed to be continuous across the interface. The following characterization is crucial.

**Proposition.** The half-spacewise linear force $f(x)$ defined in equation (2) is continuous across the interface $x^T n = 0$ if and only if the half-spacewise constant stiffness $S(x)$ satisfies the jump condition

$$S_+ = S_- + \delta \mathbf{n} \mathbf{n}^T, \quad (6)$$

where $\delta$ is a real coefficient. As a consequence, the solution $x(t)$ of equations (1) exists and is unique.

*Proof.* Continuity across the hyperplane $\mathbb{R}_+^n \cap \mathbb{R}_-^n$ requires that $\forall t | t^T n = 0, S_+ t = S_- t$. It follows that the difference matrix $S_+ - S_-$ must be singular and in fact of rank one. It can thus be expressed as the dyadic product of a vector $b$ by the normal vector $n: S_+ - S_- = \mathbf{bn}^T$. Owing to the symmetry of $S_+$ and $S_-$, it is further deduced that $b = \delta \mathbf{n}$. Conversely, post-multiplying equation (6) by $t$ yields $(S_+ - S_-)t = \mathbf{bn}^T t = 0$, which implies continuity. Now, it can be shown that the elastic force $f(x)$ is Lipschitz continuous [5], with constant max $(\|S_+\|, \|S_-\|)$; so is the function $h(y)$ with constant max $(\|H_+\|, \|H_-\|)$. It follows from the theory of differential equations that the solution $y(t)$ of equations (4) exists and is unique. Therefore, the solution $x(t)$ of equations (1) exists and is unique. $\square$

Equation (6) expresses that the jump in the stiffness matrix must be normal to the plane interface $x^T n = 0$. There is no tangential jump. The coefficient $\delta$ measures the magnitude of the jump and hence the intensity of the non-linearity. For $\delta = 0$, the linear case is recovered. The non-differentiability of $f(x)$ across an interface is characteristic of the piecewise linear systems (1) and (2). This means that all the non-linearity of such systems is concentrated at the interface, and implies that the equation cannot be linearized about the origin for the purpose of finding an approximate solution.

### 2.2.2. Positive homogeneity of $f(x)$: similarity of solutions $x(t)$

The function $f(x)$ is neither additive nor homogeneous, in contrast to a linear function. However, it remains half-spacewise (more generally conewise) additive and positively homogeneous. This has a major consequence for the solutions $x(t)$.

**Proposition.** The half-spacewise linear elastic force $f(x)$ is positively homogeneous of degree 1:

$$f(\lambda x) = \lambda f(x), \quad \forall \lambda > 0, \quad \forall x \in \mathbb{R}^n. \quad (7)$$

As a consequence, if $x(t)$ is the solution of equations (1) with $x(0) = x_0$, $\dot{x}(0) = v_0$, then $\lambda x(t)$ is the solution of equations (1) with $x(0) = \lambda x_0$, $\dot{x}(0) = \lambda v_0$,

$$M(\lambda \ddot{x}) + G(\lambda \dot{x}) + f(\lambda x) = 0, \quad \forall \lambda > 0. \quad (8)$$
*Proof.* Since $\mathbb{R}_+^n$ and $\mathbb{R}_-^n$ are half-spaces, for any coefficient $\lambda > 0$ and any $x \in \mathbb{R}_+^n$ or $x \in \mathbb{R}_-^n$, one has $\lambda x \in \mathbb{R}_+^n$ or $\lambda x \in \mathbb{R}_-^n$. Therefore, in view of its definition (2), the function $f(x)$ is positively homogeneous of degree 1, as defined in equation (7). Now, let $x(t)$ be the solution of equations (1) for $x(0) = x_0$, $\dot{x}(0) = v_0$: i.e., $M\ddot{x} + G\dot{x} + f(x) = 0$. Multiplying this equation by any $\lambda > 0$ gives, in view of equation (7), $\lambda M\ddot{x} + \lambda G\dot{x} + \lambda f(x) = M(\lambda\ddot{x}) + G(\lambda\dot{x}) + f(\lambda x) = 0$. Therefore $\lambda x(t)$ is also a solution of equations (1). Since the solution is unique, it must be the solution for the initial conditions $x(0) = \lambda x_0$, $\dot{x}(0) = \lambda v_0$. $\square$

Equation (8) shows that solutions corresponding to homothetic initial conditions are homothetic (similar). This property is the most interesting characteristic of the half-spacewise linear systems under consideration and, more generally, of conewise linear systems, which is reminiscent of the superposition principle. In phase space, positive homogeneity of the solution with respect to the initial condition is expressed by

$$y(\lambda y_0, t) = \lambda y(y_0, t), \quad \forall \lambda > 0. \quad (9)$$

### 2.2.3. Conservation of energy $q$: motion constant $c$

Pre-multiplying equation (1) by $\dot{x}^T$ yields, upon taking the anti-symmetry of $G$ into account (i.e., $\dot{x}^T G \dot{x} = 0$),

$$\dot{x}^T M \ddot{x} + \dot{x}^T S(x) x = 0.$$

It follows by integration that

$$q(x, \dot{x}) = \frac{1}{2} \dot{x}^T M \dot{x} + \frac{1}{2} x^T S(x) x = c, \quad \dot{q}(x, \dot{x}) = 0, \quad (10)$$

where $q$ is a half-spacewise quadratic form representing the total energy of the system and $c$ a motion constant depending on the initial conditions. Expression (10) implies that the system is *conservative*, and that the gyroscopic matrix $G$ does not contribute to the system energy. In phase space, equations (10) can be written as:

$$q(y) = \frac{1}{2} y^T B(y) y = c, \quad \dot{q}(y) = 0, \quad (11)$$

$$B(y) = \begin{cases} B_+, & \forall y \in \mathbb{R}_+^{2n} \\ B_-, & \forall y \in \mathbb{R}_-^{2n} \end{cases}, \quad B_+ = \begin{bmatrix} S_+ & 0 \\ 0 & M \end{bmatrix}, \quad B_- = \begin{bmatrix} S_- & 0 \\ 0 & M \end{bmatrix}.$$

### 2.2.4. Positiveness of $q(y)$: stability of the equilibrium point $y_e = 0$

Obviously, the only equilibrium point of system (1) and (2) is the origin $y_e = 0$. From equations (11), it follows that if $S_+$ and $S_-$ are positive definite, the half-spacewise quadratic form $q$ satisfies

$$q(y) > 0, \quad \dot{q}(y) = 0, \quad \forall y \neq 0. \quad (12)$$

According to stability theory, $q(y)$ is a Lyapunov function. Therefore the origin is a stable equilibrium point for the system. It should be pointed out that for gyroscopic systems ($G \neq 0$), the positiveness of $S_+$ and $S_-$ is only a sufficient condition of stability.

Having finished with the general features of conewise linear systems, we progress to a more detailed investigation of particular cases.
<figure><img src="image_4.png" /><figcaption>Figure 2. Free response of 1-DOF piecewise linear model. (a) Displacement (—), velocity (—○—) and acceleration (---); (b) phase portrait.</figcaption></figure>

### 3. PIECEWISE LINEAR 1-DOF CASE

To begin with, let us consider a piecewise linear autonomous 1-DOF model:

$$ 
\begin{gathered}
M\ddot{x} + K(x)x = 0, \quad x(0) = x_0, \quad \dot{x}(0) = v_0, \\
K(x) = K_+ > 0, \quad \forall x > 0; \quad K(x) = K_- > 0, \quad \forall x \le 0.
\end{gathered}
\qquad (13)
$$

The exact solution of equations (13) can be written as follows: for $n = 1, 2, \dots$,

$$ x(t) = A(x_0, v_0)p(t_0, t), \quad A(x_0, v_0) = \sqrt{\omega_+^2 x_0^2 + v_0^2}, \qquad (14a, b) $$

$$ 
\begin{aligned}
p(t_0, t) = & \max_{t \in (t_{2n-2}, t_{2n-1})} \left( 0, \frac{1}{\omega_+} \sin \omega_+ (t - t_{2n-2}) \right) \\
& + \min_{t \in (t_{2n-1}, t_{2n})} \left( 0, -\frac{1}{\omega_-} \sin \omega_- (t - t_{2n-1}) \right),
\end{aligned}
\qquad (14c) $$

$$ t_0 = -\frac{1}{\omega_+} \tan^{-1} \left( \frac{\omega + x_0}{v_0} \right), \quad t_{2n-1} = t_{2n-2} + \frac{\pi}{\omega_+}, \quad t_{2n} = t_{2n-1} + \frac{\pi}{\omega_-}, \qquad (14d) $$

$$ \omega_+ = \sqrt{K_+/M}, \quad \omega_- = \sqrt{K_-/M}. \qquad (14e) $$
Here $x_0$ and $v_0$ are initial conditions (with $x_0 > 0$, without loss of generality). Figure 2(a) shows the solution as a function of time, and Figure 2(b) its trajectory in phase space.

Notice that the general solution $x(t) = A(x_0, v_0)p(t_0, t)$ in equations (14) is a twice differentiable periodic piecewise sinusoidal function (the acceleration is continuous but not differentiable). Its period is given by

$$T = (T_+/2) + (T_-/2) = (\pi/\omega_+) + (\pi/\omega_-) = \pi\sqrt{M}\{(1/\sqrt{K_+}) + (1/\sqrt{K_-})\}. \quad (15)$$

An equivalent natural frequency $\omega_{eq}$ of the system can be defined as follows:

$$\omega_{eq} = 2\pi/T = 2\omega_+\omega_-/(\omega_+ + \omega_-) = 2\sqrt{K_+}\sqrt{K_-}/\sqrt{M}(\sqrt{K_+} + \sqrt{K_-}). \quad (16)$$

It is obviously a dominant component of the spectrum of $x(t)$ in the frequency domain. It was found numerically that if the system is excited by a sinusoidal force with frequency $\omega_{eq}$,

$$M\ddot{x} + K(x)x = A \sin(\omega_{eq}t + \phi), \quad (17)$$

then the maximum amplitude of the motion increases linearly to infinity (see Figure 3). Such a resonance phenomenon suggests that $\omega_{eq}$ plays a role similar to that of the natural frequency $\omega$ in a 1-DOF linear system. This is indeed the case when the piecewise linear system is excited by any function $f(t) = C\dot{p}(t_0, t)$, where $p(t_0, t)$ is defined in equations (14) and $C$ and $t_0$ are arbitrary coefficients:

$$M\ddot{x} + K(x)x = C\dot{p}(t_0, t). \quad (18)$$

The forced response of the system is:

$$x(t) = (C/2M)tp(t_0, t). \quad (19)$$

<figure><img src="image_5.png" /><figcaption>Figure 3. Resonance of 1-DOF model. (a) Displacement history; (b) phase portrait.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. 2-DOF linear model: —, first mode; —, second mode; —○—, maximum equipotential. (a) Real non-gyroscopic modes; (b) complex gyroscopic modes.</figcaption></figure>

Clearly, $x(t)$ diverges linearly as $t \to +\infty$. It is well known that under an excitation $A \sin(\omega t + \phi)$, the response of a linear system with natural frequency $\omega$ is $-At \cos(\omega t + \phi)/2$, which is comparable to equation (19) and to the graph in Figure 3(a).

### 4. COMMUTATIVE 2-DOF CASE

For further reference, let us recall some basic facts about linear systems. According to classical modal analysis theory, the general solution of a 2-DOF linear system is a linear combination of two fundamental solutions called linear modes. Each mode is composed of a harmonic eigenfunction $e^{i\omega t}$ and a corresponding eigenvector $v \neq 0$, defined by

$$(-\omega^2 M + i\omega G + S)v = 0. \quad (20)$$

Equation (20) results in real eigenvalues $\omega$ and real eigenvectors $v$ if $G=0$ and complex conjugate ones if $G \neq 0$. The system launched in the direction of an eigenvector oscillates in a sinusoidal motion on this eigenvector. In configuration space, its trajectory is a segment if $G=0$ (see Figure 4(a)) and an ellipse if $G \neq 0$ (see Figure 4(b)).

Most of the above characteristics of linear systems can be extended to a particular category of conewise linear systems; namely, non-gyroscopic ones with commutative conewise constant stiffnesses. Indeed, let $G=0$ (thus $S=K$) and $M = mI_2$ for simplicity (for $M \neq mI_2$, it suffices to replace $K_+$ and $K_-$ everywhere by $M^{-1}K_+$ and $M^{-1}K_-$ or, better, by $M^{-1/2}K_+ M^{-1/2}$ and $M^{-1/2}K_- M^{-1/2}$ to preserve their symmetry). Assume further that $K_+$ and $K_-$ are not only positive definite and continuous, but also commutative:

$$K_+ K_- = K_- K_+ . \quad (21)$$

The commutativity condition (21) and the continuity condition (6) yield $K_+ nn^T = nn^T K_+$ and $K_- nn^T = nn^T K_-$, which imply that $K_+ n = (n^T K_+ n)n$ and $K_- n = (n^T K_- n)n$, respectively. It follows that $n$ is an eigenvector of both $K_+$ and $K_-$. Due to their symmetry, the other eigenvector is $t$, orthogonal to $n$. Therefore, $V \equiv (n, t)$ is the real orthogonal matrix which diagonalizes $K_+$ and $K_-$ simultaneously,

$$K_+ = V L_+ V^T, \quad K_- = V L_- V^T, \quad (22)$$
$L_+$ and $L_-$ being diagonal matrices containing the eigenvalues $(\omega_n^{+2}, \omega_l^{+2})$ of $K_+$ and $(\omega_n^{-2}, \omega_l^{-2})$ of $K_-$, respectively. It can be shown that $L_+$ and $L_-$ are related through the continuity relation

$$L_+ = L_- + \delta e_l e_l^T, \quad e_l = (1, 0)^T, \qquad (23)$$

which implies that $\omega_n^{+2} = \omega_n^{-2} + \delta$ and $\omega_l^{+2} = \omega_l^{-2} \equiv \omega_l^2$. In this special case, the system equations (1) can be uncoupled into two independent equations by means of the linear transformation $x = Vz, z = (z_n, z_l)^T$:

$$\ddot{z}_n + \omega(z_n)z_n = 0, \quad \ddot{z}_l + \omega_l^2 z_l = 0, \qquad (24)$$

where $\omega(z_n) = \omega_n^{+2}, \forall z_n > 0, \omega(z_n) = \omega_n^{-2}, \forall z_n \le 0$. The motions in the directions $z_n$ and $z_l$ are independent (see Figure 5(a)). According to the results for the 1-DOF system, the first mode $z_n = A_n p(t_0, t)$ is piecewise sinusoidal as described in equation (14), whereas the second mode $z_l = A_l \sin(\omega_l t + \phi)$ is sinusoidal (see Figure 5(b)). It follows that any solution $x(t)$ can be represented by a linear combination of $nz_n(t)$ and $tz_l(t)$,

$$x(t) = nA_n p(t_0, t) + tA_l \sin(\omega_l t + \phi_l), \qquad (25)$$

where $A_n, A_l, t_0$ and $\phi_l$ are real numbers determined by initial conditions. Inspection of equation (25) reveals that the periodic functions $p(t_0, t)$ and $\sin(\omega_l t + \phi)$ play the role of two *eigenfunctions* associated to the two *orthogonal eigenvectors* $n$ and $t$. These eigenfunctions and eigenvectors are fundamental, since the former determine the natural frequencies of the system and the latter uncouple its equations. It is straightforward to extend these properties to MDOF commutative half-spacewise linear systems. It can be concluded that the “rectilinear mode” concept remains meaningful for these particular non-linear systems. In other terms, modal decomposition is not strictly restricted to linear systems, which are both additive and homogeneous. It also applies, in the specific form (25), to half-spacewise additive, positively homogeneous, commutative systems.

Similar motions were studied by Rosenberg in a slightly different context [35-38]. He investigated discrete, undamped, strongly non-linear, *symmetric*, *smooth* systems satisfying

$$\ddot{x} + \partial U(x)/\partial x = 0, \qquad (26)$$

<figure><img src="image_7.png" /><figcaption>Figure 5. Non-linear modes of 2-DOF commutative conewise linear model: $--$, $z_n(t)$; $—$, $z_l(t)$; $—$, maximum equipotential. (a) Configuration space; (b) time histories.</figcaption></figure>
where $x \in \mathbb{R}^n$ and $U(x): \mathbb{R}^n \to \mathbb{R}$ is a twice continuously differentiable potential energy function satisfying the central symmetry conditions $U(-x) = U(x)$ and $U(0) = 0$. He was able to characterize the periodic solutions by the following properties: all co-ordinates vary according to a periodic function (not necessarily a sinusoidal one); all displacements and velocities vanish simultaneously once every half period of the motion; all co-ordinates of the system can be expressed, at any given time, as simple-valued functions of any one co-ordinate. He called such periodic solutions *normal modes*. Notice that the periodic motions $z_n(t)$ and $z_t(t)$ of a half-spacewise linear commutative system satisfy all the above conditions, despite the fact that the system is not symmetric with respect to the origin and non-smooth. To distinguish these periodic motions from the normal modes of Rosenberg, we refer to them as *non-linear modes* for reasons that will be further explained later.

### 5. GENERAL CONEWISE LINEAR MDOF SYSTEMS

Motivated by the results of the previous sections and inspired by Rosenberg's normal mode concept, we now look for the periodic motions of more general conewise linear systems. We aim, for example, at answering the following questions: what do the periodic solutions of a general (non-commutative) half-spacewise linear 2-DOF system look like? How many are there? What are their characteristics? Since, for non-commutative matrices, $S_+$ and $S_-$, there is no regular matrix V simultaneously diagonalizing them, the system equations (1) cannot be uncoupled (even when $G=0$), and the global free response of general systems can no longer be obtained analytically. A numerical approach must be used instead.

First of all, we characterize a periodic solution of MDOF conservative systems (1) with the following five conditions of periodicity and energy conservation:

$$x(x_0, \dot{x}_0, T) = x_0, \quad \dot{x}(x_0, \dot{x}_0, T) = \dot{x}_0, \quad q(x_0, \dot{x}_0) = c. \qquad (27)$$

Here $x_0$ and $\dot{x}_0$ are the initial conditions inducing the periodic solution, $T$ is its period, $q$ is the quadratic form of the system energy and $c$ its value for $x_0$, $\dot{x}_0$. By introducing an augmented phase vector $\mathbf{a}$, and a corresponding residual vector $\mathbf{r}$ as

$$\mathbf{a} \equiv \begin{bmatrix} \mathbf{x}_0 \\ \dot{\mathbf{x}}_0 \\ \mathbf{T} \end{bmatrix}, \quad \mathbf{r} \equiv \begin{bmatrix} \mathbf{x}(\mathbf{x}_0, \dot{\mathbf{x}}_0, T) - \mathbf{x}_0 \\ \dot{\mathbf{x}}(\mathbf{x}_0, \dot{\mathbf{x}}_0, T) - \dot{\mathbf{x}}_0 \\ q(\mathbf{x}_0, \dot{\mathbf{x}}_0) - c \end{bmatrix}, \qquad (28)$$

and by using the phase space vector $\mathbf{y}$ defined in equation (5), we can put conditions (27) in the form

$$\mathbf{r}(\mathbf{a}) = \mathbf{r}(\mathbf{y}_0, T) = \begin{pmatrix} \mathbf{y}(\mathbf{y}_0, T) - \mathbf{y}_0 \\ q(\mathbf{y}_0) - c \end{pmatrix} = \mathbf{0}. \qquad (29)$$

Then, we develop an efficient numerical algorithm based on Newton's iteration method and Newmark's integration technique to extract all the periodic solutions at the predetermined energy level $c$. More specially, the solution of equations (29)—a system of five non-linear equations in five unknowns in a 2-DOF case—is obtained by applying the Newton's recurrence formula

$$\mathbf{a}^{j+1} = \mathbf{a}^j - (\partial \mathbf{r} / \partial \mathbf{a})_j^{-1} \mathbf{r}^j, \qquad (30)$$

where $j$ is the iteration counter, $\mathbf{r}^j$ and $(\partial \mathbf{r} / \partial \mathbf{a})_j$ denote the residual vector and the Jacobian matrix evaluated at $\mathbf{a}^j$. At each iteration, Newmark's implicit unconditionally stable finite difference method (trapezoidal rule),

$$\mathbf{y}_{k+1} = \mathbf{y}_k + \frac{1}{2}\Delta t(\dot{\mathbf{y}}_k + \dot{\mathbf{y}}_{k+1}), \quad k = 1, \dots, N, \quad N\Delta t = T^j, \qquad (31)$$
is used to estimate the phase vector $y' = y(a')$ after an approximate period $T'$, starting from $y_0$. The difficult part of this algorithm lies in the calculation of the Jacobian matrix (cf., Appendix A). Finally, the algorithm (30)–(31) gives us the periods $T$ and the initial conditions $y_0$ of the periodic solutions.

To demonstrate the capacity of the algorithm and study the periodic solutions, we use the 2-DOF rotating system shown in Figure 1 as a gyroscopic model. A non-gyroscopic system is obtained from the same model by setting the rotation speed to zero ($\Omega = 0$):

$$ \text{non-gyroscopic system, } \mathbf{M}\ddot{\mathbf{x}} + \mathbf{K}(\mathbf{x})\mathbf{x} = \mathbf{0}; $$

$$ \text{gyroscopic system, } \mathbf{M}\ddot{\mathbf{x}} + \mathbf{G}\dot{\mathbf{x}} + [\mathbf{K}(\mathbf{x}) - m\Omega^2\mathbf{I}_2]\mathbf{x} = \mathbf{0}. $$

$\mathbf{M}$, $\mathbf{G}$ and $\mathbf{K}(\mathbf{x})$ are given in equations (3). The next four graphs illustrate the periodic solutions of these two systems. Figure 6 presents their trajectories in configuration space, Figure 7 shows the phase portraits of $x_1$ of the solution with longest period, Figure 8 shows the time history of such a solution, and Figure 9 depicts the trajectories of two periodic solutions in the configuration space for different constants of motion.

Besides assessing the algorithm capabilities, the numerical study has revealed the following interesting characteristics of the periodic solutions of conewise linear systems.

*The trajectories of the periodic solutions of non-gyroscopic systems at a given energy level are open curves in the configuration space (instead of straight lines in the linear case), which neither pass through the origin nor are othogonal at their intersection in general (in contrast to the normal modes of Rosenberg) (see Figure 6(a)). The curvature comes from the non-linearity, the intersection offset from the lack of the central symmetry, and the loss of orthogonality is due to the lack of commutativity $S_+S_- \neq S_-S_+$. The time history of the two components of the first periodic solution (Figure 8(a)) shows that they do not vanish simultaneously, but that they do reach their maximum and minimum at the same time. In phase space, the orbits are closed curves symmetric with respect to the velocity axis, instead of ellipses (see Figure 7(a)).*

*The trajectories of the periodic solutions of gyroscopic systems are closed curves (instead of ellipses) (see Figure 6(b)). The closure is due to a phase difference between the two components of the solution, as can be seen in Figure 8(b). The deviation from ellipses, both in configuration space and in phase space, comes from the non-linearity (see Figures 6(b) and 7(b)).*

*The periods of the periodic solutions are constant and independent of the energy level (i.e., of the vibration amplitude) as for linear systems. This is because the non-linearity of the systems considered is concentrated at the interface passing through the origin.*

*The trajectories of the periodic solutions with the same period but at different energy levels form a set of homothetic curves (see Figures 9(a) and (b)). This property is due to the positive homogeneity of degree 1 of the solution with respect to the initial conditions discussed before.*

*For weak non-linearities $S_+ = S_- + \delta n n^T$, $\delta \ll \|S_-\|$, the number of structurally stable periodic solutions with different periods is equal to the number of DOF of the system, just as in the linear case. By “structurally stable” we mean that there exists no infinitesimal perturbation of the system parameters that alters the qualitative behavior of the system solution [41].*
<figure><img src="image_8.png" /><figcaption>Figure 6. Periodic trajectories in configuration space: —, period $T_1$; —, period $T_2$, —○—, maximum equipotential. (a) Non-gyroscopic system; (b) gyroscopic system.</figcaption></figure>

The characteristics of periodic solutions listed above show a close resemblance between conewise linear systems and linear systems. They suggest the following definition.

*Definition.* Any structurally stable periodic solution of a conewise linear system is called a *non-linear real mode* when $G = 0$, and a *non-linear complex mode* when $G \neq 0$. By analogy with linear modes, the periodic solutions with the *same period* and *similar trajectories* but different energy levels are said to belong to the same non-linear mode. Each non-linear mode can be considered to have a *generalized eigenvalue* that is inversely proportional to the period of the solution, and a *generalized eigenfunction* describing the non-linear modal curve.

<figure><img src="image_9.png" /><figcaption>Figure 7. Phase portraits of the periodic solution with longest period: —, $(x_1, v_1)$; —, $(x_2, v_2)$. (a) Non-gyroscopic system; (b) gyroscopic system.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8. Time histories of the periodic solution with longest period: —, $x_1$; --, $x_2$. (a) Non-gyroscopic systems; (b) gyroscopic system.</figcaption></figure>

## 6. ANALYTICAL ESTIMATE OF NON-LINEAR MODE PERIODS

Having located and characterized the non-linear real and complex modes with constant periods, we proceed to formulate an analytical estimate of these periods in terms of system parameters, using the 1-DOF system analysis as a guide. This estimate provides a simple yet precise index for crack initiation and propagation in rotating shafts.

The idea is to find an equivalent n-DOF linear system, the linear mode periods of which are equal or close to the periods of the periodic solutions of the n-DOF conewise linear system:

$$ \left\{ \begin{matrix} M\ddot{x} + G\dot{x} + S(x)x = 0 \\ x(t) = x(t + T_i) \end{matrix} \right\} \left\{ \begin{matrix} T_i = 2\pi/\omega_{eq} \\ i = 1, n \end{matrix} \right\} \left\{ \begin{matrix} M\ddot{x} + G\dot{x} + S_{eq}x = 0 \\ |-\omega_{eq}^2 M + i\omega_{eq} G + S_{eq}| = 0 \end{matrix} \right\} \quad (32) $$

<figure><img src="image_11.png" /><figcaption>Figure 9. Similarity of periodic solutions in conewise linear systems: —, period $T_1$; --, period $T_2$; —○—, maximum equipotential. (a) Non-gyroscopic system; (b) gyroscopic system.</figcaption></figure>
## 6.1. NON-GYROSCOPIC CASE

It was found in section 3 (equations (16)) that the equivalent natural frequency $\omega_{eq}$ of a 1-DOF PL system can be written as a function of the system parameters $M$, $K_+$ and $K_-$ as

$$ \omega_{eq} = 2\pi/T = 2\omega_{+}\omega_{-}/(\omega_{+} + \omega_{-}) = (2/\sqrt{M})\{(1/\sqrt{K_{+}}) + (1/\sqrt{K_{-}})\}. \quad (33) $$

This formula shows that the stiffness of the equivalent linear 1-DOF system is

$$ K_{eq} = M\omega_{eq}^2 = 4\{(1/\sqrt{K_{+}}) + (1/\sqrt{K_{-}})\}^{-2}. \quad (34) $$

This formula can be straightforwardly extended to a MDOF non-gyroscopic conewise linear system, for which $G=0$ and $S(x) = K(x)$ is the conewise constant stiffness matrix with $K_+$ and $K_-$ positive definite. Indeed, since the square root of a positive definite matrix is well defined, we can determined $S_{eq}$ by analogy with equation (34) as:

$$ S_{eq} = K_{eq} = 4(\sqrt{K_{+}}^{-1} + \sqrt{K_{-}}^{-1})^{-2}. \quad (35) $$

The spectral decompositions of $K_+$ and $K_-$ into their respective orthonormal eigenvector matrices $V_+$ and $V_-$ ($V_+^T = V_+^{-1}$, $V_-^T = V_-^{-1}$) and diagonal positive eigenvalue matrices $L_+$ and $L_-$ yields

$$ \sqrt{K_{+}} = V_{+}\sqrt{L_{+}}V_{+}^{T}, \quad \sqrt{K_{-}} = V_{-}\sqrt{L_{-}}V_{-}^{T}. $$

Consequently, $K_{eq}$ in equation (35) can be computed as

$$ K_{eq} = 4(V_{+}\sqrt{L_{+}}^{-1}V_{+}^{T} + V_{-}\sqrt{L_{-}}^{-1}V_{-}^{T})^{-2}. \quad (36) $$

Finally, the solution of the eigenvalue problem $|K_{eq} - \omega_{eq}^2 M| = 0$ gives the equivalent frequencies the inverses of which are approximations to the periods of the non-linear modes according to equations (32).

These period estimates are exact for commutative MDOF systems. To assess their accuracy in the non-commutative case, we consider at first a particular 2-DOF conewise linear system, where the natural frequencies in the two linear regions satisfy $\omega_2^+ = 3\omega_1^+$ and $\omega_2^- = 3\omega_1^-$. In this case, it can be demonstrated that one of the non-linear modes has a period equal to

$$ T = (\pi/\omega_{1}^{+}) + (\pi/\omega_{1}^{-}). \quad (37) $$

In Figure 10 it is shown that the approximate formula (36) provides a very good estimate of the non-linear mode periods for small non-linearity, while the algorithm (30) and (31) and exact formula (37) yield identical results. Incidentally, this last fact confirms that the algorithm developed in section 5 gives accurate results. The non-linear mode periods calculated from the analytical estimate (36) are also compared with the numerical results

<figure><img src="image_12.png" /><figcaption>Figure 10. Comparison of theoretical, numerical and estimation results: ——, theoretical values; ○, numerical data; ---, analytical estimates. $\omega_1^+ = \omega_1^- + \delta$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11. Analytically estimated (—) and numerically calculated (○) periods in non-gyroscopic system. (a) Period of first mode; (b) period of second mode.</figcaption></figure>

for the 2-DOF conewise linear model sketched in Figure 1 when the disc is at rest. The system parameters used are as follows:

$$ K(x) = \begin{cases} \mathbf{K}_+ = \begin{pmatrix} \lambda \cos^2 \theta + \mu \sin^2 \theta + \delta & (\lambda - \mu) \cos \theta \sin \theta \\ (\lambda - \mu) \cos \theta \sin \theta & \lambda \sin^2 \theta + \mu \cos^2 \theta \end{pmatrix} & \forall x_i \ge 0 \\ \mathbf{K}_- = \begin{pmatrix} \lambda \cos^2 \theta + \mu \sin^2 \theta & (\lambda - \nu) \cos \theta \sin \theta \\ (\lambda - \mu) \cos \theta \sin \theta & \lambda \sin^2 \theta + \mu \cos^2 \theta \end{pmatrix} & \forall x_i < 0 \end{cases}, \quad m = 1, \quad \Omega = 0, \quad \lambda = 1, \quad \mu = 6, \quad \theta = 30^\circ. \quad (38) $$

The periods calculated by the analytical estimate and the numerical iteration are plotted in Figure 11 vs. the non-linearity coefficient $\delta \in [0, 2]$. Again, for a small amount of non-linearity ($\delta < \lambda = 1$), a perfect agreement between these two methods can be observed. In this particular case, the period of the second mode varies very little with the non-linearity coefficient $\delta$.

## 6.2. GYROSCOPIC CASE

Now we consider the system with a gyroscopic anti-symmetric matrix G. In this case, the equivalent natural frequencies of the system depend not only on S(x) and M, but also on G, which is neither symmetric nor positive definite. As a consequence, S<sub>eq</sub> defined in (32) can no longer be obtained through (36). To establish a correct estimate for S<sub>eq</sub>, we put equation (1) in the form

$$ A_{\pm} y = B_{\pm} \dot{y}, \quad (39) $$

where the indexes “+” and “-” denote the two regions of linearity as usual, $B_{\pm} = B_{\pm}^{T}$ are defined in equations (11) and

$$ A_{\pm} \equiv \begin{pmatrix} \mathbf{0} & \mathbf{S}_{\pm} \\ -\mathbf{S}_{\pm} & -\mathbf{G} \end{pmatrix} = -A_{\pm}^{T}. \quad (40) $$
Without losing generality, we examine the eigenvalue problem in the positive region. It can be shown that the 2n generalized eigenvalues of the system (39) obtained from $|A_+ - \lambda B_+| = 0$ are purely imaginary and conjugate and, thus, come in pairs: $(j\omega_1^+, -j\omega_1^+), \dots, (j\omega_n^+, -j\omega_n^+)$. They correspond to n pairs of complex conjugate eigenvectors $(u_1^+, \bar{u}_1^+), \dots, (u_n^+, \bar{u}_n^+)$ which form a basis for the n-dimensional complex space. Moreover, the real and imaginary parts of these eigenvectors constitute on orthogonal basis $U_+ \equiv (u_1^{+r}, u_1^{+l}, \dots, u_n^{+r}, u_n^{+l})$, $U_+^T \neq U_+^{-1}$ for the 2n-dimensional real space satisfying

$$ U_+^T A_+ U_+ = \begin{pmatrix} \omega_1^+ E & \cdots & \omega_n^+ E \end{pmatrix} = E_\omega^+, \quad U_+^T B_+ U_+ = I_{2n}, \quad E = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}. \quad (41) $$

In equation (41), $I_{2n}$ is the 2n-dimensional identity matrix. We next perform a linear transformation on equation (39),

$$ y = U_+ z_+, \quad (42) $$

and pre-multiply it by $U_+^T$ to obtain $U_+^T A_+ U_+ z_+ = U_+^T B_+ U_+ \dot{z}_+$. According to the orthogonality relations (41), the equation is reduced to $\dot{z}_+ = E_\omega^+ z_+$, which yields a set of second order differential equations,

$$ \ddot{z}_+ + W_+ z_+ = 0. \quad (43) $$

where $W_+$ is a positive definite diagonal matrix defined by

$$ W_+ \equiv -E_\omega^{+2} = \begin{bmatrix} \omega_1^{+2} I_2 & & \\ & \ddots & \\ & & \omega_n^{+2} I_2 \end{bmatrix}. \quad (44) $$

The same procedure can be carried out in the negative region, which results in

$$ \ddot{z}_- + W_- z_- = 0. \quad (45) $$

Applying the inverse transformations $z_+ = U_+^{-1}y$ and $z_- = U_-^{-1}y$ to equations (43) and (45) respectively, yields a standard non-gyroscopic equation:

$$ \ddot{y} + R(y)y = 0, \quad R(y) = \begin{cases} \mathbf{R}_+ \equiv U_+ W_+ U_+^{-1}, & \forall y \in \mathbb{R}_+^{2n} \\ \mathbf{R}_- \equiv U_- W_- U_-^{-1}, & \forall y \in \mathbb{R}_-^{2n} \end{cases}. \quad (46) $$

The system equivalent to system (46) is then defined as

$$ \ddot{y} + R_{eq} y = 0, \quad (47) $$

where $R_{eq}$ is a function of $R_+$ and $R_-$. Since $W_+$ and $W_-$ are positive definite, the inverses of the square roots of $R_+$ and $R_-$ are equal to $U_+ \sqrt{W_+^{-1}U_+^{-1}}$ and $U_- \sqrt{W_-^{-1}U_-^{-1}}$ respectively. Therefore, $R_{eq}$ can be obtained by analogy with equation (35) as

$$ R_{eq} = 4(U_+ \sqrt{W_+^{-1}U_+^{-1}} + U_- \sqrt{W_-^{-1}U_-^{-1}})^{-2}. \quad (48) $$

Finally, it can be shown (see Appendix B) that the following two eigenvalue problems possess the same distinct eigenvalues:

$$ |-\omega^2 I_{2n} + R_{eq}| = 0 \Leftrightarrow |-\omega^2 M + i\omega G + S_{eq}| = 0. \quad (49) $$

As the dimension of the first eigenvalue problem is twice that of the second, each of its distinct eigenvalues is repeated.
<figure><img src="image_14.png" /><figcaption>Figure 12. Analytically estimated and numerically calculated periods of gyroscopic system. Analytical: —, first mode; ---, second mode. Numerical: ○, first mode; ×, second mode.</figcaption></figure>

To validate this formulation, we have re-examined numerically the previous 2-DOF gyroscopic model (Figure 1) rotating at $\Omega = 8$ rad/s, the other parameters being given by equation (38). The estimates of the non-linear modes periods have been obtained for the same interval of the non-linearity coefficient $\delta \in [0, 2]$. The analytically estimated and numerically calculated periods are presented in Figure 12. The coincidence of the two methods is even better than in the non-gyroscopic case. Note that the analytical estimate can be used as an initial guess for the numerical algorithm, which usually results in a drastic reduction in the number of iterations.

## 7. STABILITY AND BIFURCATION OF CONEWISE LINEAR SYSTEMS

The restriction on the magnitude of the non-linearity that we experienced in using the analytical formulation (36) for non-gyroscopic conewise linear systems requires further study of their characteristics. It is well known that a strongly non-linear system suffers important qualitative changes, such as instability and bifurcation of solution, chaotic motion, etc. In order to trace the entire evolution of the system characteristics as the non-linearity increases and to determine the extent to which our analytical estimate is adequate, we use again the 2-DOF system (38) to examine the stability of non-linear real modes in phase space for each value of $\delta$. For completeness, two methods, namely characteristic multipliers and the Poincaré map, are employed, although they are known to be equivalent through Floquet theory [19].

### 7.1. CHARACTERISTIC MULTIPLIERS

The stability of any non-linear mode is determined by its characteristic multipliers, also called Floquet multipliers. As a start, let us briefly recall why.

The Lyapunov definition of stability of a periodic solution $y(t, y_0^*)$ requires that the distance between the non-linear mode and the perturbed motion,

$$ \|d(t)\| = \|y(y_0^* + d_0, t) - y(y_0^*, t)\| \quad (50) $$

tends to zero when $t \to \infty$, for an arbitrary small initial perturbation $\|d(0)\| \ll \|y_0^*\|$. A Taylor expansion of $d(t)$ yields $d(t) = \{\partial y(y_0^*, t)/\partial y_0^*\}d_0$ plus terms of higher order. Because of the periodicity of $y$, $y(y_0^*, t + T) = y(y_0^*, t)$ for $t = kT$, one has

$$ d(kT) \approx \{\partial y(y_0^*, T)/\partial y_0^*\}^k d_0. \quad (51) $$
<figure><img src="image_15.png" /><figcaption>Figure 13. Bifurcation diagram of non-linear modes: —O—, first mode; ---, second mode; —, bifurcated mode. (a) Amplitudes; (b) periods.</figcaption></figure>

The 2n eigenvalues of the monodromy matrix [41, 42]

$$ \frac{\partial y(y_0^*, T)}{\partial y_0^*} \qquad (52) $$

are referred to as *characteristic multipliers*. They play an important role in deciding whether or not the initial perturbation $d_0$ decays or grows as $t \to \infty$ ($k \to \infty$). It can be observed that if $3k \in [1, 2n]$ for which $|\lambda_k| > 1$, the non-linear mode is unstable.

Notice that the monodromy matrix $\partial y/\partial y_0$ associated with a non-linear mode in expression (52) is a part of the Jacobian matrix $\partial r/\partial a$ introduced in equation (30) and evaluated at the last iteration of the numerical procedure presented there (see Appendix A for a justification). Therefore, after having located a non-linear mode by numerical iteration, it is straightforward to examine its stability.

For the previously discussed system (Figure 1, equations (38)), the variations of the maximum amplitudes and periods of the non-linear modes are plotted against the non-linearity coefficient $\delta$ in Figures 13(a) and (b) respectively. For small $\delta$, two non-linear modes determined numerically are found neutrally stable. (In this particular case, the second mode period varies very little with $\delta$.) At a critical value $\delta_c$, one of the modes becomes unstable and a third periodic solution suddenly appears.

<figure><img src="image_16.png" /><figcaption>Figure 14. Eigenvalues of monodromy matrix in complex plane. (a) Before bifurcation; (b) bifurcation; (c) after bifurcation.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

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

<figure><img src="image_19.png" /><figcaption>Figure 15. Bifurcated non-linear modes. (a) Time history: —, $x_1$; ---, $x_2$. (b) Phase portrait: —, $(x_1, v_1)$; ---, $(x_2, v_2)$. (c) Configuration space: ——, first mode; —, second mode; +—, bifurcated mode; ○, maximum equipotential.</figcaption></figure>

The instability phenomenon can be interpreted in the eigenvalue complex plane as schematized in Figure 14. For small $\delta$, all the characteristics multipliers of the mode are located on the unit circle. As $\delta$ increases, a pair of multipliers relative to one mode moves on the unit circle towards -1. At a critical value $\delta_c = 0.96$ this pair merges at -1 and $\delta > \delta_c$ one eigenvalue moves inside and the other outside the circle, on the real axis. It follows from the stability criterion that the corresponding non-linear mode becomes unstable afterwards.
The numerical results show that the additional periodic solution which appears after the critical point $\delta_c$ has a period equal to the double of that of the unstable periodic solution. This newly born periodic solution traces a different trajectory in phase space, but keeps exactly the same characteristics as the first two non-linear modes. Its time history is shown in Figure 15(a), its phase portrait in Figure 15(b), and its trajectory in configuration space together with the other two non-linear modes in Figure 15(c).

Such a bifurcation with period doubling as $\delta = \delta_c$ is a typical phenomenon, known as the flip bifurcation, which occurs in many non-linear dynamic systems [19, 24, 42, 44]. Further examination of characteristic multipliers shows that for large non-linearity ($\delta > \delta_c$), the unstable non-linear mode continues to coexist with the other two modes.

## 7.2. POINCARÉ MAP

The stability of periodic solutions can also be discussed from a more geometric point of view. Instead of examining the continuous trajectory of the system in phase space, one looks at a time-sampled sequence $z(t_1), z(t_2), \dots, z(t_n), \dots$, abbreviated by $z_1, \dots, z_n$. These points are obtained by intercepting the trajectory with a hyperplane called a Poincaré section. The stability of a periodic trajectory can then be discussed in terms of the stability of the corresponding Poincaré map.

Since the definition of this map relies on the knowledge of the solution of the differential equation (4), the sequence $z_n$ cannot be computed unless a general solution of this equation is available. In our autonomous conewise linear systems, the trajectory in any one cone can be obtained analytically, as long as the entering condition is known. Thus, for establishing the Poincaré map for a half-spacewise linear MDOF system, it is natural to use the hyperplane $\Sigma$ across which the system stiffness is discontinuous as a Poincaré section.

For illustration we use a 2-DOF model and suppose that $\mathbf{m} = (1 \ 0 \ 0 \ 0)^T$:

$$ \Sigma = \mathbb{R}_+^4 \cap \mathbb{R}_-^4 : \{y \in \mathbb{R}^4 | y_1 = 0\}. \quad (53) $$

If the solution point hits the hyperplane $2N$ times within a period (see Figure 16), then it has travelled $N$ times in each half-space. In view of the closed form solution in each halfspace, three successive points $z_{2k-2}, z_{2k-1}$ and $z_{2k}$ on the Poincaré section are related through two square $3 \times 3$ matrices [5],

$$ z_{2k-1} = P(A_+, t^{2k-2})z_{2k-2}, \quad z_{2k} = P(A_-, t^{2k-1})z_{2k-1}, \quad (54) $$

where the travelling times can be found in terms of the entering positions $z_{2k-1}$ and $z_{2k-2}$:

$$ t^{2k-2} = t^{2k-2}(z_{2k-2}), \quad t^{2k-1} = t^{2k-1}(z_{2k-1}), \quad k = 1, 2, \dots \quad (55) $$

<figure><img src="image_20.png" /><figcaption>Figure 16. Poincaré section and trajectory.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 17. Maximum eigenvalues of equations (22) and (28): —, by characteristic multipliers; ○, by Poincaré section.</figcaption></figure>

Differentiating equations (54) with respect to the entering positions yields:

$$ (\partial z_{2k-1} / \partial z_{2k-2}) = P(A_+, t^{2k-2}) + \dot{P}(A_+, t^{2k-2}) z_{2k} (\partial t^{2k-2} / \partial z_{2k-2}) \equiv Z_+^{2k-2}, \quad (56a) $$

$$ (\partial z_{2k} / \partial z_{2k-1}) = P(A_-, t^{2k-1}) + \dot{P}(A_-, t^{2k-1}) z_{2k-1} (\partial t^{2k-1} / \partial z_{2k-1}) \equiv Z_-^{2k-1}. \quad (56b) $$

Upon compiling these leaps, the Jacobian matrix DP of the Poincaré map takes the form

$$ DP = Z_-^{2N-1} Z_+^{2N-2} \dots Z_-^{2k-1} Z_+^{2k-2} \dots Z_+^0. \quad (57) $$

After localizing the travelling times $t_{2k-2}$ and $t_{2k-1}$ and intercepting points $z_{2k-2}$ and $z_{2k-1}$ ($k = 1, N$), the matrices $P(A_+, t^{2k-2})$, $P(A_-, t^{2k-1})$, $Z_+^{2k-2}$, $Z_-^{2k-1}$ and DP can be successively evaluated. They are explicitly given in reference [5].

The same criterion of stability as for the characteristic multipliers is then applied to the eigenvalues of DP. To return to the example (38), exactly the same instability phenomenon is observed as with multipliers. Another neutrally stable mode with double period appears. The maximum eigenvalue of the bifurcating mode is replotted against $\delta$ in Figure 17 to show the perfect agreement between the two methods.

To recapitulate, a 2-DOF conewise linear system has two neutrally stable non-linear modes below a non-linearity threshold ($\delta < \delta_c$) and it can have three or more beyond it.

## 8. CONCLUSIONS

Mechanical systems involving unilateral contacts (whether microscopic or macroscopic) give rise to conewise linear models. In particular, half-spacewise linear models have been widely adopted, as a first order approximation, for describing the transverse vibrations of cracked beams (non-gyroscopic case) and rotating shafts (gyroscopic case) [1-5]. The present study sheds some light on the multi-degree-of-freedom models. The main results are summarized here.

(i) The existence and uniqueness of the solution has been established on the grounds of a continuity condition at the interface and the similarity of solutions with similar initial conditions (at different energy levels) has been proved on the basis of the positive homogeneity of the governing force.

(ii) Below a critical non-linearity or bifurcation threshold ($\delta < \delta_c$), conewise linear systems bear a strong resemblance to linear systems, in terms of their fundamental periodic solutions herein called non-linear modes whether real (non-gyroscopic case) or complex (gyroscopic one): their number is equal to the number of DOF, their periods are
constant and energy independent, their trajectories are similar curves at different energy levels, etc.

(iii) Beyond the bifurcation threshold ($\delta > \delta_c$), a bifurcation of the second non-linear mode occurs and the resemblance to linear systems is reduced. The threshold is much higher in the gyroscopic case ($\delta_c > 2$, [5]) than in the non-gyroscopic one ($\delta_c = 0.96$) due to the well known stabilizing effect of rotation [42].

(iv) The analytical estimate of the non-linear mode periods is very accurate below the bifurcation threshold and, hereby, could become a valuable tool in engineering practice. Indeed, for crack detection in rotating shafts, it gives the precise relationship between the natural frequencies and the system parameters, and thus shows how the former will vary with the crack depth [5]. Note that even in the non-gyroscopic case, the bifurcation threshold is situated far enough to interfere with crack detection concerns.

## ACKNOWLEDGMENT

We wish to thank Mrs M. Schatzman for showing us how to shorten the proof of the second proposition concerning the positive homogeneity of the solutions.

## REFERENCES

1. K. B. COLLINS, R. H. PLAUT and J. WAUER 1991 *Transactions of the American Society of Mechanical Engineers, Journal of Vibration and Acoustics* **113**, 74-78. Detection of cracks in rotating Timoshenko shaft using axial impulses.

2. R. GASCH 1976 *Institution of Mechanical Engineers Conference*, 123-128. Dynamical behavior of a simple rotor with a cross-sectional crack.

3. B. GRABOWSKI 1980 *Journal of Mechanical Design* **102**, 140-146. The vibrational behavior of a turbine rotor containing a transverse valve.

4. J. WAUER 1990 *International Journal of Solid and Structures* **26**, 901-914. Modeling and formulation of equations of motion for cracked rotating shafts.

5. L. ZUO 1992 *Ph.D. Thesis No. 1006, Ecole Polytechnique Fédérale de Lausanne, Switzerland*. Etude du comportement dynamique des systèmes linéaires par morceaux—contribution à la détection des fissures dans les arbres de machines tournantes.

6. A. TONDL 1965 *Some Problems of Rotor Dynamics*. London: Chapman & Hall.

7. I. W. MAYES and W. G. R. DAVIES 1976 *Institution of Mechanical Engineering Conference*, 53-64. The vibration behavior of a rotating shaft system containing a transverse crack.

8. F. ISMAIL, A. IBRAHIM and H. R. MARTIN 1990 *Journal of Sound and Vibration* **140**, 305-317. Identification of fatigue cracks from vibration testing.

9. S. TIMOSHENKO 1974 *Vibration Problems in Engineering*. New York: John Wiley.

10. S. F. MASRI, Y. A. MARIAMY and J. C. ANDERSON, 1981 *Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics* **48**, 404-410. Dynamic response of a beam with a geometric nonlinearity.

11. S. NATSIAVAS 1989 *Journal of Sound and Vibration* **134**, 315-331. Periodic response and stability of oscillators with symmetric trilinear restoring force.

12. S. NATSIAVAS 1990 *International Journal of Non-linear Mechanics* **25**(5), 535-554. On the dynamics of oscillators with bi-linear damping and stiffness.

13. S. MAEZAWA 1961 *Bulletin of the Japan Society of Mechanical Engineers* **4**(14), 201-212. Steady forced vibration of unsymmetrical piecewise-linear system.

14. S. MAEZAWA and S. FURUKAWA 1973 *Bulletin of the Japan Society of Mechanical Engineers* **16**(96), 931-941. Superharmonic resonance in piecewise-linear system.

15. S. MAEZAWA, H. KUMANO and Y. MINAKUCHI 1980 *Bulletin of the Japan Society of Mechanical Engineers* **23**(175), 68-75. Forced vibration in an unsymmetric piecewise-linear system excited by general periodic force functions.

16. R. J. COMPARIN and R. SINGH 1989 *Journal of Sound and Vibration* **134**, 259-290. Non-linear frequency response characteristics of an impact pair.
17. Y. S. CHOI and S. T. NOAH 1988 *Journal of Sound and Vibration* **121**, 117-126. Forced periodic vibration of unsymmetric piecewise-linear systems.

18. Y. B. KIM and S. T. NOAH 1991 *Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics* **58**, 545-553. Stability and bifurcation analysis of oscillators with piecewise-linear characteristics: a general approach.

19. J. GUCKENHEIMER and P. HOLMES 1989 *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*. Applied Mathematical Sciences 42. New York: Springer-Verlag.

20. A. J. LICHETENBERG and M. A. LIEBERMAN 1983 *Regular and Stochastic Motion*. Applied Mathematical Sciences 38. New York: Springer-Verlag.

21. S. W. SHAW and P. J. HOLMES 1983 *Journal of Sound and Vibration* **90**, 129-155. A periodically forced piecewise linear oscillator.

22. J. M. T. THOMPSON, A. R. BOKAIAN and R. GHAFFARI 1983 *IMA Journal of Applied Mathematics* **31**, 207-234. Subharmonic resonances and chaotic motions of a bilinear oscillator.

23. J. M. T. THOMPSON, A. R. BOKAIAN and R. GHAFFARI 1984 *Journal of Energy Resources Technology* **106**, 191-198. Subharmonic and chaotic motions of complaint offshore structures and articulated mooring towers.

24. J. M. T. THOMPSON and H. B. STEWART 1986 *Nonlinear Dynamics and Chaos*. Chichester: John Wiley.

25. J. M. T. THOMPSON and S. BISHOP 1988 *Physics Bulletin* **39**, 232-234. From Newton to chaos.

26. I. A. MAHFOUZ and F. BADRAKHAN 1990 *Journal of Sound and Vibration* **143**, 255-288. Chaotic behavior of some piecewise linear systems, part I: systems with set-up spring or with unsymmetric elasticity.

27. J. N. SCHULMAN 1983 *Physical Review A* **28**(1), 477-479. Chaos in piecewise linear systems.

28. L. O. CHUA, M. HASLER, J. NEIRYNCK and P. VERBURGH 1982 *IEEE Transactions on Circuits and Systems*, **29**, 534-546. Dynamics of a piecewise linear resonant circuit.

29. T. MATSUMOTO 1987 *IEEE Proceedings* **75**, 1033-1057. Chaos in electronic circuits.

30. L. O. CHUA, M. KOMURO and T. MATSUMOTO 1986 *IEEE Transactions on Circuits and Systems* **33**, 1072-1118. The double scroll family. Part I: rigorous proof of chaos. Part II: rigorous analysis of bifurcation phenomena.

31. J. S. N. ELVEY 1983 *IMA Journal of Applied Mathematics* **31**, 235-252. On the elimination of de-stabilizing motion of articulated mooring tower under steady sea conditions.

32. J. M. T. THOMPSON and J. S. N. ELVEY 1984 *International Journal of Mechanical Science* **26**(6-8), 419-426. Elimination of sub-harmonic resonances of compliant marine structures.

33. S. W. SHAW 1989 *International Union of Theoretical and Applied Mechanics Conference. Stuttgart, Germany*. The Suppression of chaos in periodically forced oscillators. In *Nonlinear Dynamics in Engineering Systems*.

34. R. J. COMPARIN and R. SINGH 1990 *Journal of Sound and Vibration* **142**, 101-124. Frequency response characteristics of a multi-degree-of-freedom system with clearances.

35. R. M. ROSENBERG 1960 *Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics*, 263-268. Normal modes of nonlinear dual-mode systems.

36. R. M. ROSENBERG 1962 *Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics*, 7-14. The normal modes of nonlinear N-degree-of-freedom systems.

37. R. M. ROSENBERG and J. K. KUO 1964 *Transactions of the American Society of Mechanical Engineers, Journal of Applied Mechanics*, 283-290. Nonsimilar normal mode vibrations of nonlinear systems having two degrees of freedom.

38. R. M. ROSENBERG 1966 in *Advances in Applied Mechanics*, 155-242; New York: Academic Press. On nonlinear vibrations of systems with many degrees of freedom. Academic press.

39. J. MAWHIN 1964 *Bulletin de la Société Royale des Sciences de Liège* **33**(9-10), 540-557. Oscillations en modes normaux de systèmes dynamiques non-linéaires à plusieurs degrés.

40. R. T. ROCKAFELLAR 1972 *Convex Analysis*. Princeton, New Jersey: Princeton University Press.

41. T. S. PARKER and L. O. CHUA 1989 *Practical Numerical Algorithms for Chaotic Systems*. New York: Springer-Verlag.

42. S. RÜDIGER 1988 *From Equilibrium to Chaos—Practical Bifurcation and Stability Analysis*. New York: Elsevier.

43. K. HUSEYIN 1978 in *Mechanics of Elastic Stability* (H. LEIPHOLZ, editor). Alphen aan den Rijn: Sitthoff & Noordhoff. Vibrations and stability of multiple parameter systems.

44. F. C. MOON 1987 *Chaotic Vibration—An Introduction for Applied Scientists and Engineers*. New York: John Wiley.
## APPENDIX A: CALCULATION OF JACOBIAN MATRIX ($\partial r/\partial a$) AND MONODROMY MATRIX $\partial y(y_0^*, T)/\partial y_0^*$

The definitions of $r$, $a$ and $q$ in section 5 imply the following equalities: $\partial y_0/\partial y_0 = I_{2n} \in \mathbb{R}^{2n \times 2n}$; $\partial y_0/\partial T = 0 \in \mathbb{R}^{2n \times 1}$. The Jacobian matrix in equation (30) takes therefore the form

$$ \frac{\partial r}{\partial a} = \frac{\partial}{\partial a} \begin{pmatrix} y(y_0, T) - y_0 \\ q(y_0) - c \end{pmatrix} = \begin{pmatrix} \frac{\partial y(y_0, T)}{\partial y_0} - I_{2n} & \frac{\partial y(y_0, T)}{\partial T} \\ \frac{\partial q(y_0)}{\partial y_0} & 0 \end{pmatrix}. \quad (A1) $$

Moreover, in equation (A1)

$$ \frac{\partial y(y_0, T)}{\partial T} = \frac{\partial y(y_0, T)}{\partial T} = \dot{y}(y_0, T) \quad (A2) $$

and from equation (11)

$$ \frac{\partial q(y_0)}{\partial y_0} = \begin{cases} y_0^T B_+, & \forall y_0 \in \mathbb{R}_+^{2n} \\ y_0^T B_-, & \forall y_0 \in \mathbb{R}_-^{2n} \end{cases} \quad (A3) $$

If we denote the monodromy matrix as

$$ \frac{\partial y(y_0, T)}{\partial y_0} \equiv \mathbf{J} \in \mathbb{R}^{2n \times 2n}, \quad (A4) $$

the Jacobian matrix (A1) becomes:

$$ \frac{\partial \mathbf{r}}{\partial a} = \begin{pmatrix} \mathbf{J}(y_0, T) - I_{2n} & \dot{y}(y_0, T) \\ \frac{\partial q(y_0)}{\partial y_0} & 0 \end{pmatrix}. \quad (A5) $$

In order to obtain $\mathbf{J}$, we first differentiate equation (4) with respect to $y_0$ to obtain

$$ \dot{\mathbf{J}} = \mathbf{H}(y)\mathbf{J} \quad (A6) $$

and then solve simultaneously this matrix equation and the system equation (4):

$$ \dot{\mathbf{y}} = \mathbf{H}(y)\mathbf{y}, \quad \dot{\mathbf{J}} = \mathbf{H}(y)\mathbf{J}, \quad \mathbf{y}(0) = \mathbf{y}_0, \quad \mathbf{J}(0) = I_{2n}. \quad (A7) $$

Here we use the fact that $\mathbf{J}(0) = \partial y(y_0, 0)/\partial y_0 = I_{2n} \in \mathbb{R}^{2n \times 2n}$. $\mathbf{H}(y)$ is defined in equation (5).

Integration of equation (A7) allows us to evaluate the square matrix (A4), the column vector (A2) and the line vector (A3) at time $T$ and the initial condition $y_0$. The Jacobian matrix (A1) follows. For this purpose, we employ the unconditionally stable trapezoidal rule (31). We do not use the closed form solution in each linear region because the number of times a periodic solution crosses the border line is unknown; therefore, the switching points are not defined *a priori*.

Finally, with the predicted values $y_0$, and $T$, we obtain the Jacobian matrix (A1), as well as the residual $\mathbf{r} = \mathbf{y}(y_0, T) - y_0$ for each iteration of equation (30). Once the initial condition leading to a periodic solution $y_0^*$ and its period $T$ are determined, the last evaluation of equation (A4) at $y_0^*$ and $T$ gives the monodromy matrix:

$$ \frac{\partial y(y_0^*, T)}{\partial y_0^*} = \frac{\partial y}{\partial y_0} \Big|_{y_0^*, T}. \quad (A8) $$

## APPENDIX B: ON THE EIGENVALUES OF THE EQUIVALENT NON-GYROSCOPIC SYSTEM $\ddot{y} + R_{eq}y = 0$

The $n$-DOF gyroscopic system equivalent to the piecewise linear system (1) is defined as

$$ M\ddot{\mathbf{x}} + G\dot{\mathbf{x}} + S_{eq}\mathbf{x} = 0. \quad (B1) $$

Another equivalent $2n$-DOF system represented by a non-gyroscopic equation is

$$ \ddot{\mathbf{y}} + \mathbf{R}_{eq} \mathbf{y} = 0. \quad (B2) $$
We now demonstrate that the eigenvalues of system (B2) are double, and that its distinct eigenvalues are identical to that of the system (B1).

Transforming equation (B1) from configuration space to phase space yields

$$A_{eq} y = B_{eq} \dot{y}, \qquad (B3)$$

$$y = \begin{pmatrix} \mathbf{x} \\ \dot{\mathbf{x}} \end{pmatrix}, \quad A_{eq} = \begin{pmatrix} \mathbf{0} & \mathbf{S}_{eq} \\ -\mathbf{S}_{eq} & -\mathbf{G} \end{pmatrix}, \quad \mathbf{B}_{eq} = \begin{pmatrix} \mathbf{S}_{eq} & \mathbf{0} \\ \mathbf{0} & \mathbf{M} \end{pmatrix}.$$

Taking the time derivative of (B3) further gives

$$\mathbf{B}_{eq} \ddot{\mathbf{y}} = \mathbf{A}_{eq} \dot{\mathbf{y}} = \mathbf{A}_{eq} \mathbf{B}_{eq}^{-1} \mathbf{B}_{eq} \dot{\mathbf{y}} = \mathbf{A}_{eq} \mathbf{B}_{eq}^{-1} \mathbf{A}_{eq} \mathbf{y}. \qquad (B4)$$

Comparing equations (B2) and (B4) gives the equality

$$\mathbf{R}_{eq} = -\mathbf{R}_{eq}^{-1} \mathbf{A}_{eq} \mathbf{B}_{eq}^{-1} \mathbf{A}_{eq}. \qquad (B5)$$

According to linear complex modal analysis theory, the eigenvalues of the problem

$$A_{eq} \mathbf{u} = \lambda \mathbf{B}_{eq} \mathbf{u} \qquad (B6)$$

are purely imaginary and conjugate: $\lambda = i\omega_k$, $-i\omega_k$ ($k = 1, n$) and they are associated with $n$ pairs of complex conjugate eigenvectors $\mathbf{u} = \mathbf{u}_k, \bar{\mathbf{u}}_k$. Separation of the real and imaginary parts of $\lambda$ and $\mathbf{u}$ yields

$$\begin{cases} \mathbf{B}_{eq}^{-1} \mathbf{A}_{eq} \mathbf{u}'_k = -\omega_k \mathbf{u}'_k \\ \mathbf{B}_{eq}^{-1} \mathbf{A}_{eq} \mathbf{u}'_k = \omega_k \mathbf{u}'_k \end{cases}, \quad k = 1, \dots, n. \qquad (B7)$$

In view of the definitions of the matrices $\mathbf{U}, \mathbf{E}_\omega$ and $\mathbf{W}$ in equations (41) and (44), equations (B7) take the matrix form

$$\mathbf{B}_{eq}^{-1} \mathbf{A}_{eq} \mathbf{U} = \mathbf{U} \mathbf{E}_{\omega}. \qquad (B8)$$

Finally, $\mathbf{R}_{eq}$ can be expressed as

$$\mathbf{R}_{eq} = -(\mathbf{U}\mathbf{E}_{\omega}\mathbf{U}^{-1})(\mathbf{U}\mathbf{E}_{\omega}\mathbf{U}^{-1}) = \mathbf{U}(-\mathbf{E}_{\omega}^2)\mathbf{U}^{-1} = \mathbf{U}\mathbf{W}\mathbf{U}^{-1}, \qquad (B9)$$

where the matrix $\mathbf{W}$ is diagonal and positive definite:

$$\mathbf{W} = -\mathbf{E}_{\omega}^2 = \begin{bmatrix} \omega_1^2 \mathbf{I}_2 & 0 \\ 0 & \ddots \\ \omega_n^2 \mathbf{I}_2 \end{bmatrix}. \qquad (B10)$$

It follows that the eigenvalue problem for $\mathbf{R}_{eq}$ has double roots:

$$\begin{aligned} |\mathbf{R}_{eq} - \omega^2 \mathbf{I}_{2n}| &= |\mathbf{U}\mathbf{W}\mathbf{U}^{-1} - \omega^2 \mathbf{I}_{2n}| = |\mathbf{W} - \omega^2 \mathbf{I}_{2n}| \\ &= \begin{vmatrix} (\omega_1^2 - \omega^2)\mathbf{I}_2 & \cdots & (\omega_n^2 - \omega^2)\mathbf{I}_2 \end{vmatrix} = (\omega_1^2 - \omega^2)^2 \cdots (\omega_n^2 - \omega^2)^2 = 0. \quad (B11) \end{aligned}$$