# Impact modes in discrete vibrating systems with rigid barriers

Valery Pilipchuk

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

## To cite this version:

Valery Pilipchuk. Impact modes in discrete vibrating systems with rigid barriers. International Journal of Non-Linear Mechanics, 2001, 36 (6), pp.999-1012. 10.1016/S0020-7462(00)00066-4. hal-01295245

HAL Id: hal-01295245

https://hal.science/hal-01295245

Submitted on 30 Mar 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_2.png" /></figure>
# Impact modes in discrete vibrating systems with rigid barriers

V.N. Pilipchuk

*Mechanical Engineering, Wayne State University, Detroit, MI 48202, USA*

*E-mail address: pylyp@me1.eng.wayne.edu (V.N. Pilipchuk)*

## Abstract

A class of strongly non-linear vibrating systems composed of linear elastic structures under absolutely rigid constraints condition is considered. Impact mode exact solutions are expressed through a saw-tooth time argument and, as a result, represented in a closed unit form. Based on this special representation, sufficient conditions of existence and also non-existence for the impact modes are formulated and interpreted on the spectral axes. In particular, it was shown that impact modes exist when their basic frequencies are shifted into the right small neighborhood of any natural frequency of the linearized (no barriers) system. The frequencies of the localized impact mode solutions are located at the right of the linear spectrum and have no upper boundary.

*Keywords: Normal modes; Impact systems; Non-linear localization; Non-smooth transformations*

## 1. Introduction

In many practically important cases, the displacements of a vibrating mechanical system can be spatially bounded by rigid barriers. This can be caused by the system design and/or some gained clearances due to deterioration of joints. Both situations were considered by many investigators. The related engineering applications were summarized in research monographs quite completely, see for example [1–3]. If the total energy of the system is sufficient to reach the barriers, the system immediately becomes strongly non-linear even when the differential equations of motion between the barriers are linear. As a result, even structurally simple

models can perform a very complex dynamical behavior [4,5]. This is why general analytical solutions for many degrees of freedom vibroimpact systems are rather impossible. However, an analytical insight into the multidimensional impacting systems can be reached by imposing certain restrictions on the class of motions. In this paper, we consider special vibrating regimes when all the system’s particles synchronously vibrate as a single oscillator. In linear cases, such assumption leads to the normal mode solutions, and hence, gives a complete information about vibrating systems due to the linear superposition principle. In non-linear cases, this principle does not work, however, special non-linear normal mode solutions appear to be very important. An overview of the related analytical methods as well as some original material can be found in [6]. An important physical feature is that non-linear modes frequently perform a “trend” to spatial localization [7]. From the point of view of applied mechanics, a literature survey and an
overview of different analytical methods can be found in [6,8] Spatial symmetries, stability, and localization of non-linear modes in the two- and three degrees of freedom systems with bilateral perfectly rigid barriers were considered in [9] combining the special transformation of coordinates [10] and the idea of averaging. Recently, the non-linear normal modes concept was also adapted for two degrees of freedom non-linear systems with non-linear regular (polynomial) and singular (potential barriers) components in [11].

If clearances between the barriers are sufficiently small, the differential equations of motion between the barriers can be assumed linear. Although the system still remains strongly non-linear due to the constraints condition, one can employ the temporal symmetry of the periodic motion in order to match different pieces of the linear solution between the barriers, and thus, obtain some kind of exact normal mode solutions [12]. Exact unit-form periodic solutions were obtained in recent work [13] for a two degrees of freedom “vibroimpact” system. In principal, the methodology of work [13] is based on the saw-tooth transformation of time adapted for a class of impulsively excited systems [14] (see Conclusion for related remarks).

In certain sense, the present work combines the ideas of works [12,13]. As a result, exact unit-form analytical solutions for a multidegrees of freedom system with rigid barriers are obtained. In particular, conditions of existence and localization of impact modes, and also their physical meaning are considered. The unit form solutions could be suitable for analytical manipulations with them in different perturbations schemes, including averaging techniques.

In this work, the *impact modes* of a discrete vibrating system with bilateral rigid barriers are periodic motions when some or may be all particles of the system interact with the barriers twice per one period.

It will be shown that the case of unilateral constraints can be considered as well.

## 2. An introductory example

Let us consider a one degree of freedom free harmonic oscillator between two absolutely rigid

<figure><img src="image_3.png" /><figcaption>Fig. 1. The oscillator with bilateral rigid barriers (a) is replaced by the oscillator under the periodic series of external impulses (b).</figcaption></figure>

barriers. A mechanical model of the oscillator can be represented as a massive particle attached to the free end of a massless cantilever beam as it is shown in Fig. 1a.

The interaction with the barriers at $x = \pm \Delta$ is assumed to be perfectly elastic, and the system is represented in the form

$$ \ddot{x} + \omega_0^2 x = 0, \quad |x| \le \Delta \quad (1) $$

Since the normal mode regimes are periodic by their definition, the reaction of constraints can be treated as a periodic series of external impulses acting on the masses of the system.

Applying this remark to the one degree of freedom system as it is shown in Fig. 1, the related
differential equation of motion is written in the linear form

$$ \ddot{x} + \omega_0^2 x = 2p \sum_{k=-\infty}^{\infty} [\delta(\omega t + 1 - 4k + \varphi) - \delta(\omega t - 1 - 4k + \varphi)] = p\tau''(\omega t + \varphi) \quad (2) $$

where $\delta(\xi)$ is Dirac function; $2p$, $\omega$, and $\varphi$ will be kept as arbitrary constants; $\tau(\xi)$ is a standard saw-tooth sine

$$ \tau(\xi) = \begin{cases} \xi, & -1 \le \xi \le 1 \\ -\xi, & 2, 1 \le \xi \le 3 \end{cases} \quad \forall \xi \quad \tau(\xi) = \tau(\xi + 4) \quad (3) $$

This piece-wise linear periodic function can be also represented in the unit form $\tau(\xi) = (2/\pi)\arcsin\sin(\pi\xi/2)$. The related amplitude and the period are normalized in such a manner that the expression $[\tau'(\xi)]^2 = 1$ holds at least for almost all$^2 \xi \in (-\infty, \infty)$. For further convenience, the right-hand side of Eq. (2) is expressed through second-order generalized derivative of the saw-tooth sine, where the derivative is taken with respect to the whole argument, $\omega t + \varphi$. The parameter $\omega$ will be called a frequency parameter; though the standard circular frequency is expressed as $\Omega = (\pi/2)\omega$. In this paper, both parameters, $\omega$ and $\Omega$, will be used.

In contrast to system (1), the auxiliary system (2) is linear. Representing unknown steady-state periodic solution in the form,

$$ x = X(\tau), \quad \tau = \tau(\omega t + \varphi) \quad (4) $$

one obtains the boundary value problem with no singular terms,

$$ \omega^2 X''(\tau) + \omega_0^2 X(\tau) = 0, \quad X'(\tau)|_{\tau=\pm 1} = p\omega^{-2} \quad (5) $$

and the related solution is represented in the saw-tooth time form [14]

$$ x = \frac{p \sin[(\omega_0/\omega)\tau]}{\omega\omega_0 \cos(\omega_0/\omega)}. \quad (6) $$

This solution can be verified by direct substitution of expression (6) into the equation of motion (2).

A connection between solution (6) and vibration of the original system with the rigid constraints is established by imposing the conditions (Fig 1):

* The impulses at the right-hand side of Eq. (2) must act at those time instances at which the mass is in contact with the constraints, i.e.,

$$ x = \pm \Delta \quad \text{when } \tau = \pm 1 \leftrightarrow \tau'' \neq 0 \quad (7) $$

* The system cannot penetrate through the barriers, i.e.,

$$ |x| \le \Delta \quad \text{for all } \tau \in [-1, 1]. \quad (8) $$

Substituting solution (6) into condition (7), one determines the impulses parameter $p$,

$$ p = \Delta\omega\omega_0 \cot(\omega_0/\omega) \quad (9) $$

for which solution (6) takes the final form

$$ x = \Delta \frac{\sin[(\omega_0/\omega)\tau(\omega t + \varphi)]}{\sin(\omega_0/\omega)}. \quad (10) $$

Solution (10) satisfies condition (7) automatically. The related parameter $p$ (9) will further be treated as an “eigen-value” of the non-linear (impact) problem.

The parameters, $\omega$ and $\varphi$, can be expressed through the initial conditions. Let us assume that $x(0) = 0$, i.e., $\varphi = 0$, and hence, the total energy of the oscillator per unit mass is expressed through the initial velocity as $E = [\dot{x}(0)]^2/2$. Substitution (10) into the last equality gives expression for the frequencies ratio,

$$ \frac{\omega_0}{\omega} = \pm \frac{1}{2} \arccos \left( 1 - \frac{\omega_0^2 \Delta^2}{E} \right) + k\pi, \quad k = 0, 1, \dots \quad (11) $$

The right-hand side gives a sequence of real numbers when the total energy, $E$, is sufficiently large, $E \ge E_* = \omega_0^2 \Delta^2/2$, so that the oscillator can reach the constraints.

Generally speaking, not all the magnitudes of the parameter $\omega$ given by sequence (11) lead to real motions of the original system with rigid barriers. Indeed, since the auxiliary system of type (2) does

<footnote n="2">²A periodic series of points, {$\xi$: $\tau(\xi) = \pm 1$}, at which first derivative $\tau'(\xi)$ has a step wise discontinuity, may give an exemption. However, it is not important because the differential equations of motion must be treated in terms of distributions whenever the Dirac $\delta$-function is involved in the equations (see [20]).</footnote>
not have any constraints for the coordinate, $x(t)$, hence condition (7) does not guarantee that this coordinate will remain inside the admissible region, $|x| \le \Delta$, during the whole period of vibration. That is why another condition (8) must be verified as well. Such a verification implemented for solution (10) shows that condition (8) is satisfied only for the smallest root in set (11). Fig. 2 illustrates the temporal mode shapes corresponding to the first two roots $\omega_0/\omega$.

It is seen that the second solution, which is shown in fragment (b), violates condition (8), although it satisfies condition (7). Such “phantom” solutions, however, can be provided with a certain physical meaning if we imagine removable barriers appearing periodically twice per one period and reflecting the particle into the exterior of the region $|x| \le \Delta$.

In this paper, the case of unilateral barriers will not be considered specially. The related investigation could be implemented, however, in the same way. Let us remove, for example, the left barrier and consider the oscillator (1) under the unilateral constraint condition $x \le \Delta$. In this case, the boundary conditions in (5) should be modified as $X'(\tau)|_{\tau=\pm 1} = \pm p\omega^{-2}$. Such periodic change of sign effectively switches the directions of positive $\delta$-pulses to opposite at the right-hand side of Eq. (2).

<figure><img src="image_4.png" /><figcaption>Fig. 2. Real (a) and “phantom” (b) solutions corresponding to the first (smallest) and second roots respectively, $(\omega_0/\omega)_1 = 1.1502$ and $(\omega_0/\omega)_2 = -1.1502 + \pi$. The total energy level is $E = 1.2E_*$.</figcaption></figure>

As a result, the solution takes the form,

$$x = \Delta \frac{\cos[(\omega_0/\omega)\tau(\omega t + \varphi)]}{\cos(\omega_0/\omega)} \quad (12)$$

where the period of the solution, $x(t)$, and the related basic frequency will be $T = 2/\omega$ and $\Omega = \pi\omega$, respectively, due to evenness of cosine. Analysis of the temporal mode shapes for different frequency ratio $\omega_0/\omega$ is not a difficult task.

### 3. Impact modes

Let us consider $N$-degrees of freedom conservative system described by the coordinates vector $\mathbf{x} = (x_1, \dots, x_N)^T \in \mathbb{R}^N$. We suppose that one of the system coordinates, say the $a$th one, is restricted as $|x_a| \le \Delta_a$. In the vector notation, this can be rewritten as

$$|\mathbf{I}_a^T \mathbf{x}| \le \Delta_a \quad (13)$$

where $\mathbf{I}_a$ is the $a$th vector of the basis of the system coordinates, $\mathbf{I}_a = (0, \dots, 1, \dots, 0)^T$.

Inside the domain (13), the differential equations of motion are assumed to be linear,

$$M\ddot{\mathbf{x}} + K\mathbf{x} = 0 \quad (14)$$

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Mass-spring model of a discrete elastic system with bilateral barriers.</figcaption></figure>

where M and K are respectively constant mass and stiffness N × N-matrices.

A mass-spring model of the system (13) and (14) is shown in Fig. 3.

Note that the form of matrix equation (14) does not suppose the system to be necessarily a mass-spring chain, see Example 1 later.

To obtain an explicit analytical solution for the impact modes, the essentially non-linear system (13) and (14) is replaced by an impulsively forced linear system under no constraints condition

$$M\ddot{x} + Kx = pI_a\tau''(\omega t + \varphi) \quad (15)$$

where p is an unknown “eigenvalue”, the constants $\omega$ and $\varphi$ are going to be kept as arbitrary parameters. Thus one seeks a two parameter family of periodic solutions.

A family of periodic solutions of the period $4/\omega$ can be found as a linear superposition of solutions (6) for each of the N modes of system (14) with a relevant replacement of the parameters [14],

$$x = p \sum_{j=1}^{N} \frac{(e_j^T I_a) e_j \sin[(\omega_j/\omega)\tau]}{\omega \omega_j \cos(\omega_j/\omega)} \quad (16)$$

where $\tau = \tau(\omega t + \varphi)$ is the sawtooth temporal argument, $e_j$ and $\omega_j$ are the jth normal mode and the natural frequency of linear system (3.2). The linear normal modes are normalized as $e_j^T M e_j = \delta_{ji}$, where $\delta_{ji}$ is the Kronecker symbol. We suppose that $\omega_i < \omega_j$ when $i < j$.

The impulses must act at those time instances when the a-th mass interacts with the constraints. Such spatially-temporal “location of the impacts” can be assigned by condition

$$I_a^T x = \pm \Delta_a \quad \text{when } \tau = \pm 1. \quad (17)$$

Substituting (16) into (17), one obtains the related “eigenvalue,”

$$p = \Delta_a \omega \left[ \sum_{j=1}^{N} (e_j^T I_a)^2 \frac{\tan(\omega_j/\omega)}{\omega_j} \right]^{-1} \quad (18)$$

where $e_j^T I_a$ is simply the a-th component of the j-th linear mode vector. Substituting (18) into (16) gives a two parameter family of the periodic solutions for the impact modes. The parameter $\varphi$ is an arbitrary phase shift, whereas the frequency parameter $\omega$ implies some restrictions due to condition (13).

**Proposition 1.** Sufficient conditions for the impact modes existence can be summarized as follows:

* The basic frequency of the vibration, $\Omega = (\pi/2)\omega$, pertains to a sufficiently small right neighborhood of any frequency $\omega_j$ of the linear spectrum, i.e.,

$$\Omega = \omega_j/(1 - \sigma), \quad 0 < \sigma \le 1 \quad (19)$$

* For all $\omega_i$ such that $i \ne j$ the linear frequencies ratios $\omega_i/\omega_j$ are not odd numbers, $\omega_i/\omega_j \ne 2k - 1$ ($k = 1, 2, ...$).

Under these conditions, solution (16), (18), satisfies inequality (13), and hence, describes a two parameter family of the real impact mode regimes, $x = x(t; \omega, \varphi)$.

An idea of the proof is to find out such cases when $x$ is a monotonic function of $\tau$ on the interval $-1 \le \tau \le 1$, and hence condition (17) at the boundaries guarantees that inequality (13) holds inside the whole interval as well, see Appendix for further details of the proof. The above conditions are sufficient but not necessary. Namely, if the conditions are true then the impact mode exists. Existence of an impact mode does not guarantee, however, that the conditions hold. For example, some impact modes can exist in the resonance case $\omega_i/\omega_j = 2k - 1$, when the parameter,

$$\lambda = (e_i^T I_a / e_j^T I_a)^2 / (2k - 1)^2 \quad (20)$$

is sufficiently small. This parameter shows a relative weight of the non-monotonic term which appears in the expression due to the resonance (see appendix).

Example 2 below shows that the impact modes may possess a quite complicated spectral structure.
Thus, the global investigation is needed to formulate necessary and sufficient conditions. On the other hand, a sufficient condition of non-existence can be formulated by making use the physical meaning of the impact interaction parameter $p$.

**Proposition 2.** If an impact mode exists then the inequality $p(\omega) > 0$ holds.

**Proof.** The parameter $p$ (18) cannot be negative for any real impact vibrating regime as a reaction of constraint, because it cannot be directed to a barrier. Thus impact modes cannot exist when $p(\omega) < 0$.

The above Propositions 1 and 2 are illustrated below on special examples.

**Example 1.** The top mass $m_1$ of a free doublependulum (see Fig. 4) oscillates between the two absolutely rigid constraints providing small angle amplitude of the top pendulum to the vertical,

$$|\varphi_1| \le \Delta_1 \ll 1.$$

Assuming that another angle, $\varphi_2$, is sufficiently small as well, one can take the differential equations of motion between the barriers in the linear form

$$\mu^2 \ddot{\varphi}_1 + \ddot{\varphi}_2 + \mu^2 \varphi_1 = 0,$$

$$\ddot{\varphi}_1 + \ddot{\varphi}_2 + \varphi_2 = 0,$$

where $\mu^2 = 1 + m_1/m_2$, overdot denotes differentiation with respect to a new time parameter, $\bar{t} = (g/l)^{1/2}t$, and $l$ is the length of the rods.

In this case, the linear modes and natural frequencies are, respectively,

$$\mathbf{e}_1 = \frac{1}{\sqrt{2\mu(\mu+1)}} \begin{pmatrix} 1 \\ \mu \end{pmatrix}, \quad \mathbf{e}_2 = \frac{1}{\sqrt{2\mu(\mu-1)}} \begin{pmatrix} 1 \\ -\mu \end{pmatrix}$$

and

$$\omega_1 = \sqrt{\frac{\mu}{\mu+1}}, \quad \omega_2 = \sqrt{\frac{\mu}{\mu-1}}$$

where the modal vectors satisfy the orthogonality condition, $\mathbf{e}_j^T \mathbf{M} \mathbf{e}_i = \delta_{ji}$, with respect to the mass matrix

$$M = \begin{pmatrix} \mu^2 & 1 \\ 1 & 1 \end{pmatrix}.$$

<figure><img src="image_7.png" /><figcaption>Fig. 4. Double pendulum with bilateral constraints.</figcaption></figure>

Taking into account that $\mathbf{I}_a = \mathbf{I}_1$, $\Delta_a = \Delta_1$, and $N=2$, the formal solution (16) and (18) for the angle coordinates is reduced to

$$\begin{pmatrix} \varphi_1 \\ \varphi_2 \end{pmatrix} = \Delta_1 \left[ \tan(\omega_1/\omega) + (\omega_2/\omega) \tan(\omega_2/\omega) \right]^{-1} \\ \times \left[ \begin{pmatrix} 1 \\ \mu \end{pmatrix} \frac{\sin\left[(\omega_1/\omega)\tau(\omega t + \varphi)\right]}{\cos(\omega_1/\omega)} + \frac{\omega_2}{\omega_1} \right. \\ \left. \times \begin{pmatrix} 1 \\ -\mu \end{pmatrix} \frac{\sin\left[(\omega_2/\omega)\tau(\omega t + \varphi)\right]}{\cos(\omega_2/\omega)} \right].$$

Let us take the parameter $\mu$ such that

$$\frac{\omega_2}{\omega_1} = \sqrt{\frac{\mu+1}{\mu-1}} = 2k-1.$$

In this case, the parameter $\lambda$ given by expression (20) is evaluated exactly as $\lambda = 1$. This value is relatively small (see appendix) so that the existence of impact modes is not affected, if condition (19) holds.

It is important to note that the first impact mode basic frequency $\Omega$ must be sufficiently close to the first linear frequency $\omega_1$, so that it is not approaching the left neighborhood of the next frequency, $\omega_2$. Since there are no more frequencies after $\omega_2$, hence its right “neighborhood” has no upper boundary
and can be extended to infinity, $\omega_2 < \Omega < \infty$. As a result, the highest impact mode becomes spatially localized when the parameter $\Omega$ is growing. The localization can be estimated explicitly by asymptotic ($\Omega^{-1} \to 0$) expansion

$$ \begin{aligned} \left. \frac{\varphi_2}{\varphi_1} \right|_{\tau=1} &= \mu \left( \tan \frac{\pi \omega_1}{2\Omega} - \frac{\omega_2}{\omega_1} \tan \frac{\pi \omega_2}{2\Omega} \right) / \left( \tan \frac{\pi \omega_1}{2\Omega} + \frac{\omega_2}{\omega_1} \tan \frac{\pi \omega_2}{2\Omega} \right) \\ &= \mu \frac{\omega_1^2 - \omega_2^2}{\omega_1^2 + \omega_2^2} \left[ 1 + \frac{\pi^2 \omega_1^2 \omega_2^2}{6(\omega_1^2 + \omega_2^2)} \Omega^{-2} + O(\Omega^{-4}) \right] \\ &= -1 - \frac{\pi^2}{12} \Omega^{-2} - O(\Omega^{-4}). \end{aligned} $$

It is seen that $(\varphi_2/\varphi_1)|_{\tau=1} \to -1$ as $\Omega \to \infty$, i.e., the horizontal displacement of the bottom mass becomes negligibly small whereas the impacting mass oscillates between the barriers. Fig. 5 represents the related trend of the mode shape.

**Example 2.** Let us consider a mass-spring chain of $N$ identical particles

$$ \ddot{x}_n + \frac{k}{m}(-x_{n-1} + 2x_n - x_{n+1}) = 0, \quad (21) $$

$$ n = 1, \dots, N; \quad x_0 = x_{N+1} = 0 $$

under the bilateral constraints condition (see Fig. 3),

$$ |x_a| \le \Delta_a \quad (22) $$

where $k$ and $m$ are stiffness of the springs and mass of the particles, respectively.

A set of no impact modes and frequencies can be found explicitly, since the system can be viewed as a finite element model of a continuous string. The $j$th mode vector and related frequency of the linear system are given by the expressions

$$ \mathbf{e}_j = \frac{4}{N+1} \left( \sin \frac{\pi j}{N+1}, \dots, \sin \frac{N \pi j}{N+1} \right)^{\mathrm{T}} $$

and

$$ \omega_j = \omega_{N+1} \sin \frac{\pi j}{2(N+1)} \quad (23) $$

where notation<fnref n="3" /> $\omega_{N+1} = 2\sqrt{k/m}$ has been introduced, the basis vectors are normalized such that $\mathbf{e}_j^\mathrm{T}\mathbf{e}_i = \delta_{ji}$.

The $a$th component of the $j$th normal mode vector is

$$ \mathbf{e}_j^\mathrm{T} \mathbf{I}_a = \frac{4}{N+1} \sin \frac{a \pi j}{N+1}. \quad (24) $$

Let us show that the impact mode periodic solution (16) and (18) becomes localized as the frequency parameter is infinitely growing, $\omega \to \infty$. First, replacing the trigonometric functions by their asymptotic expansions as $\sin(\omega_j \tau/\omega) \sim \omega_j \tau/\omega$, $\cos(\omega_j/\omega) \sim 1$, and $\tan(\omega_j/\omega) \sim \omega_j/\omega$, we obtain an estimate for an arbitrary $b$th particle impact mode motion in the form

$$ x_b \sim \Delta_a \frac{\sum_{j=1}^{N} (\mathbf{e}_j^\mathrm{T} \mathbf{I}_a)(\mathbf{e}_j^\mathrm{T} \mathbf{I}_b)}{\sum_{j=1}^{N} (\mathbf{e}_j^\mathrm{T} \mathbf{I}_a)^2} \tau(\omega t + \varphi) \quad \text{as } \omega \to \infty. \quad (25) $$

In this particular example, the sums can easily be evaluated. Taking into account expression (24) and the standard trigonometric sums [15], one obtains

$$ \begin{aligned} \sum_{j=1}^{N} (\mathbf{e}_j^\mathrm{T} \mathbf{I}_a)(\mathbf{e}_j^\mathrm{T} \mathbf{I}_b) &= \frac{16}{(N+1)^2} \sum_{j=1}^{N} \sin \frac{a \pi j}{N+1} \sin \frac{b \pi j}{N+1} \\ &= \frac{8}{N+1} \delta_{ab}. \end{aligned} \quad (26) $$

Later on, we will need also the sum

$$ \sum_{j=1}^{N} (\mathbf{e}_j^\mathrm{T} \mathbf{I}_a)(\mathbf{e}_j^\mathrm{T} \mathbf{I}_b) \omega_j^2 = \frac{4\omega_{N+1}^2}{N+1} \left( -\delta_{a,b-1} + 2\delta_{a,b} - \delta_{a,b+1} \right). \quad (27) $$

Substituting (26) into the asymptotic expression (25), one obtains

$$ x_b \sim \Delta_a \delta_{ab} \tau(\omega t + \varphi) \quad \text{as } \omega \to \infty. \quad (28) $$

This expression shows that the $a$th particle of the chain vibrates according to the saw-tooth temporal mode shape with an infinitely large frequency, whereas all other particles remains at rest. Thus, the impact mode becomes spatially localized as $\omega \to \infty$. The localization tendency of the impact modes is illustrated by Fig. 5a.

<footnote n="3">³ Note that this is just a suitable notation since the $(N+1)$th frequency does not physically exist.</footnote>
<figure><img src="image_8.png" /></figure>

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

<figure><img src="image_10.png" /><figcaption>Fig. 5. Localized impact modes of the mass-spring chain ($n = 25$) for (a) single constrained particle (b) two constrained particle and the same frequency $\Omega = \omega_{25} + 0.01$. The same case of two constrained particles when increasing the frequency is illustrated by (b) and (c) $\Omega = \omega_{25} + 0.2$.</figcaption></figure>

Fig. 6 provides Propositions 1 and 2 with some additional illustration on the spectral axes of the chain consisting of $N = 5$ particles. The figure represents diagrams of the impulsive force parameter,

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

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

<figure><img src="image_13.png" /><figcaption>Fig. 6. The impulsive force parameter, $p$, as a function of the trigonometric frequency, $\Omega$, for the first (a), second (b), and third (c) mass under the constraint condition. The chain consists of $N = 5$ masses. For all cases the clearance parameter is $\Delta = 0.5$, and $\omega_{N+1} = \omega_6 = 1$.</figcaption></figure>

$p$, as a function of the basic circular frequency of vibration, $\Omega = (\pi/2)\omega$, for different masses under the constraints, $a = 1, 2$, and $3$. It is seen that each of the frequencies of the linear spectrum possesses its right neighborhood where the parameter $p$ is positive, and hence the impact modes can exist (see fragments (a)-(c)). At the same time, the parameter $p$ is not necessarily negative on the left of the
frequencies $\omega_1, \dots, \omega_5$ as it is seen from the fragments (b) and (c) (see, respectively, the frequency $\omega_3$ and the frequencies $\omega_2$ and $\omega_4$). The diagrams structure depends on the particle under the constraint condition. An interesting feature of the diagrams is that they are more complicated around the lowest frequency. The parameter $p$ can change its sign many times in the domain between the first two frequencies, $\omega_1 < \Omega < \omega_2$.

### 4. Systems with multiple impacting particles

Let us consider the case of two constrained particles. It will be seen that a generalization on multiple constrained particles is a quite trivial task. For example, let us put the constraints condition on the $a$th and $b$th coordinates ($a \neq b$) of the general linear system (14). These conditions are $|x_a| \le \Delta_a$ and $|x_b| \le \Delta_b$ or in the vector notations

$$|\mathbf{I}_a^T \mathbf{x}| \le \Delta_a \quad \text{and} \quad |\mathbf{I}_b^T \mathbf{x}| \le \Delta_b. \qquad (29)$$

In this case, the impulsive excitation on the right-hand side of auxiliary equation must act on the both $a$th and $b$th particles, so that the equation takes the form

$$M\ddot{\mathbf{x}} + K\mathbf{x} = (p_a\mathbf{I}_a + p_b\mathbf{I}_b)\tau''(\omega t + \varphi) \qquad (30)$$

where $p_a$ and $p_b$ are parameters to be determined. The related solution includes some items related to $p_a$ and $p_b$, and can be represented in the form

$$\mathbf{x} = \sum_{j=1}^{N} \frac{p_a(\mathbf{e}_j^T \mathbf{I}_a)\mathbf{e}_j + p_b(\mathbf{e}_j^T \mathbf{I}_b)\mathbf{e}_j}{\omega \omega_j} \frac{\sin(\omega_j \tau / \omega)}{\cos(\omega_j \tau / \omega)} \qquad (31)$$

Following the normal modes ideology, we assume that the impact mode periodic regime is accompanied by synchronous impact interactions of the both impacting particles with the constraints according to conditions

$$\mathbf{I}_a^T \mathbf{x} = \pm \Delta_a \quad \text{and} \quad \mathbf{I}_b^T \mathbf{x} = \pm \Delta_b \quad \text{when } \tau = \pm 1. \qquad (32)$$

Substitution (31) into (32) gives linear algebraic equations with respect to $p_a$ and $p_b$ in the form

$$\begin{aligned} k_{aa}p_a + k_{ab}p_b &= \Delta_a, \\ k_{ba}p_a + k_{bb}p_b &= \Delta_b, \end{aligned} \qquad (33)$$

where

$$k_{ab} = \sum_{j=1}^{N} \frac{(\mathbf{e}_j^T \mathbf{I}_a)(\mathbf{e}_j^T \mathbf{I}_b)}{\omega \omega_j} \tan \frac{\omega_j}{\omega}. \qquad (34)$$

Expressions (31) and (33) give a formal impact mode solution. The impact mode (31) exists for those frequencies $\omega$ at which the determinant of system (33) is non-zero, and also condition (29) holds. Solution (31) can be viewed as a strongly non-linear superposition of the two simplest impact modes with a single impact pair. The result of the superposition is shown in Figs. 5b and c for two different magnitudes of the frequency, $\Omega$.

Let us consider the higher frequency domain for the above example of mass-spring system, when $\Omega \gg \omega_N$. In this case, one has $\tan(\omega_j/\omega) > 0$ for all $j = 1, \dots, N$, and hence, the coefficients $k_{ab}$ (34) create the Gram matrix with non-zero determinant [16]. The asymptotic estimation below confirms this conclusion.

Indeed, for $\Omega \gg \omega_N$ all ratios $\omega_j/\omega$ are small and thus an asymptotic expansion holds as

$$\frac{\sin[(\omega_j/\omega)\tau]}{\cos(\omega_j/\omega)} = \frac{\omega_j \tau}{\omega} + \frac{1}{2} \left(\frac{\omega_j}{\omega}\right)^3 \left(\tau - \frac{\tau^3}{3}\right) + O(\omega^{-5}). \qquad (35)$$

Substituting this expression into solution (31) and coefficients (34), one obtains the asymptotic solution for arbitrary $c$th particle in the form

$$x_c = (P_a \delta_{ac} + P_b \delta_{bc}) \tau - \frac{\pi^2}{16} \left(\frac{\omega_{N+1}}{\Omega}\right)^2 \left(\tau - \frac{\tau^3}{3}\right) \times \left[ P_a(\delta_{a,c-1} - 2\delta_{ac} + \delta_{a,c+1}) + P_b(\delta_{b,c-1} - 2\delta_{bc} + \delta_{b,c+1}) \right] + O(\Omega^{-4}) \qquad (36)$$

where $x_c = \mathbf{x}^T \mathbf{I}_c$ and $\tau = \tau(2\Omega t / \pi + \varphi)$.

The impact interaction parameters, $P_a$ and $P_b$, are given by the linear algebraic system

$$\begin{aligned} K_{aa}P_a + K_{ab}P_b &= \Delta_a, \\ K_{ba}P_a + K_{bb}P_b &= \Delta_b, \end{aligned} \qquad (37)$$

with the coefficients

$$K_{ab} = \delta_{ab} - \frac{\pi^2}{24} \left(\frac{\omega_{N+1}}{\Omega}\right)^2 (\delta_{a,b-1} - 2\delta_{ab} + \delta_{a,b+1}) + O(\Omega^{-4}). \qquad (38)$$
The parameters and the coefficients have been
rescaled as follows:

$$
\begin{align*}
\{P_a, P_b\} &= \frac{2\pi^2}{(N+1)\Omega^2} \{p_a, p_b\}, \\
\{K_{aa}, K_{ab}, K_{ba}, K_{bb}\} &= \frac{(N+1)\Omega^2}{2\pi^2} \{k_{aa}, k_{ab}, k_{ba}, k_{bb}\}.
\end{align*}
$$

It is seen that $K_{ab} \rightarrow \delta_{ab}$ as $\Omega \rightarrow \infty$, and Eqs. (37) give the solution $P_a = \Delta_a$ and $P_b = \Delta_b$. In this limit, according to (36), the vibration becomes localized on the two particles which vibrate between the barriers performing the saw-tooth temporal shape,

$$
x_c \sim (\Delta_a \delta_{ac} + \Delta_b \delta_{bc}) \tau \left( \frac{2\Omega}{\pi} t + \varphi \right) \quad \text{as } \Omega \to \infty.
$$

This is the leading-order asymptotic approach.
All other particles, $c \neq a$ and $c \neq b$, are at rest.
However, they oscillate with small amplitudes of
different orders of $\Omega^{-2}$ when the parameter $\Omega$ is
bounded. Asymptotic expansion (36) includes the
terms of order $\Omega^{-2}$. It is seen that the particles
numbered $c = a \pm 1$ and $c = b \pm 1$ have the ampli-
tudes of the lowest order, $\Omega^{-2}$. These particles are
the nearest to the impacting ones. Some analysis of
expansion (36) shows that the temporal mode
shapes of the impacting particles, $c = a$ and $c = b$,
are non-smooth and close to the saw-tooth sine
wave as it is expected to be due to the interactions
with the constraints. On the other hand, the tem-
poral mode shapes of the nearest particles,
$c = a \pm 1$ and $c = b \pm 1$, are estimated by the com-
bination of the saw-tooth sine, $\tau - \tau^3/3$, with the
first two continuous derivatives with respect to
time, $t$. These are

$$
\frac{d}{dt}\left(\tau-\frac{\tau^{3}}{3}\right)=\frac{2 \Omega}{\pi}\left(1-\tau^{2}\right) \tau^{\prime} \in C^{1}\left(R\right)
$$

and

$$
\begin{align*}
\frac{d^2}{dt^2}\left(\tau - \frac{\tau^3}{3}\right) &= \left(\frac{2\Omega}{\pi}\right)^2 \left[-2\tau(\tau')^2 + \underbrace{(1-\tau^2)\tau''}_{\phantom{2\tau(\tau')^2}}\right] \\
&= -\left(\frac{2\Omega}{\pi}\right)^2 2\tau \in C(R),
\end{align*}
$$

respectively. The prime denotes differentiation with
respect to the whole argument of the saw-tooth

sine, $\tau' \equiv d\tau/d(2\Omega t/\pi + \varphi)$.

In the above expressions, one should take into
account that the δ-type singularities of second de-
rivative, τ", are "located" at those points {t}, where
τ = ± 1. As a result, the underlined above term is
zero.<fnref n="4" />

Thus the temporal mode shapes of the particles
with no constraints are smooth functions of time,
although they are expressed through the saw-tooth
sine.

5. Conclusion

A class of strongly non-linear many degrees of
freedom vibrating systems including bilateral rigid
barriers was considered. Exact analytical periodic
solutions were expressed through the saw-tooth
temporal argument. The conditions of existence
and non-existence for impact mode regimes were
formulated based on the exact unit-form solutions.
More precisely, Proposition 1 gives a sufficient
condition of existence of impact modes, whereas
Proposition 2 gives a sufficient condition of their
non-existence. These two propositions do not com-
pose any necessary and sufficient condition of exist-
ence for the impact modes. Namely, one cannot
indicate *all the cases* when impact modes exist, and
also *all the cases* when they do not. A certain
exemption is given only by the high frequency vi-
bration Ω > ω<sub>N</sub>, when the impact modes become
localized spatially and always exist. The examples
clearly showed that a complete investigation pres-
ents a global geometrical problem. This, however,
can be done at least numerically, since the exact
expressions have been obtained in terms of trigon-
ometric sums. The problem of stability can be con-
sidered by evaluating the eigenvalues of the
Jacobian matrix for known periodic solutions as it
was done, for example, in [17]. The systems above
considered are strongly non-linear due to barriers
only. In some cases, the system displacements must
be significant in order to provide an interac-
tion with barriers, and hence, linearization of the

<footnote n="4">⁴I.e., it gives zero contribution into the related integrals of the theory of distributions.</footnote>
equations of motion between the barriers may be-
come doubtful. In general terms, however, the pro-
posed formulation remains applicable. For
example, Eq. (2) can be generalized as

$$
\ddot{x} + f(x, \dot{x}, t) = p\tau''(\omega t)
$$

provided that the function $f(x, \dot{x}, t)$ is periodic ($T = 4/\omega$) with respect to time, and the condition $|x| \le \Delta$ holds. In this case, the saw-tooth temporal argument can be introduced by using the two-component representation for periodic solutions [14,18].

$$
x(t) = X(\tau(\omega t)) + Y(\tau(\omega t))\tau'(\omega t).
$$

This results into the coupled non-homogeneous
boundary-value problem with respect to unknown
functions $X(\tau)$ and $Y(\tau)$,

$$
\omega^2 X''(\tau) + R_f(X, Y, X', Y', \tau) = 0,
$$

$$
X'(\tau)|_{\tau=\pm 1} = p\omega^{-2},
$$

$$
\omega^2 Y''(\tau) + I_f(X, Y, X', Y', \tau) = 0, \quad Y|_{\tau=\pm 1} = 0,
$$

where functions $R_f$ and $I_f$ are obtained by substi-
tution of the above representation for periodic
solutions into the equation of motion. Such a kind
of boundary value problems admit approximate
analytical solutions in the power series form with
respect to the new temporal argument, $\tau$, since it
automatically accounts for the periodicity of the
motion. Unfortunately, the related solutions may
appear to be so complicated that it will be either
difficult or even impossible to separate real and
"phantom" solutions.

To this end, the above manipulations with the
sawtooth temporal variable, should not be confused
with the method of non-smooth transformation of
the configuration space coordinates for the systems
with rigid barriers [10]. These two methods are
mutually antithetic in its physical meaning as well
as in its mathematical implementation. For
example, being applied to system (1), the method of
work [10] would give a non-linear oscillator (com-
pare with (5)),

$$
\ddot{s} + \omega_0^2 \tau'(s) \tau(s) = 0
$$

where $s(t)$ is a new positional variable, such that
$x = \tau(s)\Delta$; the result is written in the notations of
present work for the reason of comparison. The

later transformation does not require the motions
to be periodic. Unfortunately, this advantage
brings essential non-linearity of the transformed
differential equations. Further comparison of these
two techniques as well as their composition can be
found in [19].

Acknowledgements

The author would like to thank Professor A.A.
Zevin from TRANSMAG Research Institute of
Academy of Sciences of Ukraine for complete text
of reference [12], and a very helpful discussion.

**Appendix A**

**Proof of Proposition 1.** On the first stage, taking
into account the expansion [15, p. 43],

$$
\frac{\tan \lambda}{\lambda} = \frac{8}{\pi^2} \sum_{n=1}^{\infty} [(2n-1)^2 - (2\lambda/\pi)^2]^{-1}
$$

for $\lambda = \omega_j/\omega$, one can represent solution (16) and
(18) in a more convenient for our purposes form

$$
\begin{equation}
\begin{split}
\mathbf{x} = \Delta_a \left[ \sum_{j,n} \frac{(\mathbf{e}_j^\mathrm{T} \mathbf{I}_a)^2}{(2n-1)^2 - (\omega_j/\Omega)^2} \right]^{-1} \\
\times \sum_{j,n} \frac{(\mathbf{e}_j^\mathrm{T} \mathbf{I}_a) \mathbf{e}_j}{(2n-1)^2 - (\omega_j/\Omega)^2 \sin[(\pi/2)\omega_j/\Omega]}.
\end{split}
\tag{A.1}
\end{equation}
$$

Generally speaking, it is difficult to conclude
whether this function of τ is monotonic or not, since
it is composed of many oscillating terms. However,
when ωⱼ/Ω ~ 2n - 1 for certain Ω, ωⱼ and n, then
the only one term is expected to be dominating in
each of the two sums involved in solution (A.1). In
this case, the result can be explicitly estimated by
considering the neighborhood of the related pole of
the numerator and denominator

$$
\frac{\omega_j}{\Omega} = 2n - 1 - \sigma, \quad |\sigma| \ll 1 \qquad (A.2)
$$

where the parameter $\sigma$ shows how the both parts of
ratio (A.1) are close to their pole.
Substituting (A.2) into (A.1) and taking the *ath* component of the vector **x**, one obtains an estimation in the form

$$x_a = \mathbf{x}^T \mathbf{I}_a = \Delta_a \frac{\sin[\tau(2n - 1 - \sigma)\pi/2]}{\sin[(2n - 1 - \sigma)\pi/2]} + O(\sigma). \quad (A.3)$$

In this expression, the leading term is monotonic on the interval $-1 \le \tau \le 1$ and its absolute value does not exceed $\Delta_a$ if only $n=1$ and the parameter $\sigma$ is positive and sufficiently small. In this case expression (A.2) can be rewritten in the form (19), corresponding to the first condition of Proposition 1. To complete this part of the proof, however, one has to specially investigate a small neighborhood of the boundary points $\tau = \pm 1$. The neglected terms $O(\sigma)$ in (A.3) may become significant when $\tau$ is close to any of the two boundaries and the impacting mass coordinate in the range $|x_a - \Delta_a| \sim \sigma$ or $|x_a + \Delta_a| \sim \sigma$. Let us consider a small neighborhood of the boundary $\tau = 1$, were $\tau = 1 - \Delta\tau$ ($0 \le \Delta\tau \le 1$). Substituting this expression for $\tau$ and also the frequency $\Omega$ given by (19) into the exact solution (A.1), and keeping the leading order terms with respect to $\sigma$ and $\Delta\tau$, gives the estimation

$$x_a = \mathbf{x}^T \mathbf{I}_a = \Delta_a \left[ 1 - \frac{\pi^2 \sum_{j=1}^N (\mathbf{e}_j^T \mathbf{I}_a)^2}{4 (\mathbf{e}_j^T \mathbf{I}_a)^2} \sigma \Delta \tau \right].$$

Since the coefficient of the product $\sigma\Delta\tau$ is negative, it is seen that condition (13) holds, i.e., the coordinate $x_a$ remains inside the domain $|x_a| \le \Delta_a$ provided that the values $\sigma$ and $\Delta\tau$ are sufficiently small.

Now, let us check whether or not the sums involved in the exact solution (A.1) possess any poles different than (19) as $\sigma \to 0$. First, taking into account expression (19), one obtains

$$\begin{aligned} (2k - 1)^2 - \left(\frac{\omega_i}{\Omega}\right)^2 &= (2k - 1)^2 - \left(\frac{\omega_i}{\omega_j}\right)^2 (1 - \sigma)^2 \\ &= (2k - 1)^2 - \left(\frac{\omega_i}{\omega_j}\right)^2 \\ &\quad + (2\sigma - \sigma^2) \left(\frac{\omega_i}{\omega_j}\right)^2. \end{aligned} \quad (A.4)$$

It is seen that this value goes to zero as $\sigma \to 0$ when $k=1$ and $i=j$. This is the above considered case. However, one may have another zero limits when the frequencies ratio is an odd number, i.e., when $\omega_i/\omega_j = 2k-1$. In this cases the sums in (A.1) will have more the one term, and one cannot guarantee the monotonicity of $x_a$ as a function of $\tau$. So we eliminate these cases from consideration by the second of the above two sufficient conditions. The proof is completed.

The above expression “one cannot guarantee” does not determine those cases when the second condition of Proposition 1 is really significant. To investigate this point, let us suppose that the frequency ratio is odd number, i.e., the resonance ratio

$$\frac{\omega_i}{\omega_j} = 2k - 1 \quad \text{for certain } k \ge 2 \quad (A.5)$$

holds. In this case, the terms of order $O(1)$ on the right-hand side of expression (A.4) disappear. Then, substituting (19) and (A.4) into the exact solution (A.1), one obtains estimation

$$x_a = \mathbf{x}^T \mathbf{I}_a = \frac{\Delta_a}{1 + \lambda} \left[ \frac{\sin[\tau(1 - \sigma)\pi/2]}{\sin[(1 - \sigma)\pi/2]} + \lambda \frac{\sin[\tau(2k - 1)(1 - \sigma)\pi/2]}{\sin[(2k - 1)(1 - \sigma)\pi/2]} \right] + O(\sigma) \quad (A.6)$$

where the parameter $\lambda$ (20) depends on the system parameters and also numbers $i, j, k$ and $a$.

The first term of expression is of the same structure as (A.3) and appears to be monotonic for $-1 \le \tau \le 1$ and $0 < \sigma \le 1$. The second term arises due to condition (A.5), and it is not a monotonic function of $\tau$ for $k \ge 2$. Thus, under the resonance condition (A.5), it is not known whether the function $x_a(\tau)$ reaches its absolute maximum and minimum at the ends of the interval $-1 \le \tau \le 1$. Fig. 7 shows, however, that the answer essentially depends on the parameter $\lambda$.

The three fragments of this figure represent solution (A.6) for the impacting mass as a function of the two arguments, $x_a = x_a(\tau, \lambda^{-1})$, for two different odd frequency ratios, $\omega_2/\omega_1 = 3$ and 5. An inverse scale has been chosen for the parameter $\lambda$ in order to better visualize the regimes of the
surface, intersect the horizontal (non-meshed) planes, $x_a/\Delta_a = \pm 1$. These planes represent the two barriers. As follows from expressions (A.6) and (20), the parameter $\lambda$ determines a contribution of the $i$th (resonating) mode into the solutions and thus is an important characteristic of the system. Namely, the $i$th mode is responsible for violation of the condition $|x_a/\Delta_a| \le 1$. Fig. 7 also shows that increasing the frequency ratio $\omega_2/\omega_1$ and the $\sigma$-distance from the basic frequency $\omega_1$ will expand the $\lambda^{-1}$-region of non-existing “phantom” solutions.

<figure><img src="image_14.png" /><figcaption>(a) $\frac{\omega_2}{\omega_1} = 5, \sigma = 0.1$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) $\frac{\omega_2}{\omega_1} = 5, \sigma = 0.15$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) $\frac{\omega_2}{\omega_1} = 3, \sigma = 0.2$</figcaption></figure>

Fig. 7. Family of solutions (shown by meshed surfaces) for different odd frequency ratios, $\omega_2/\omega_1$, and $\sigma$-distances from the basic frequency, $\omega_1$. The horizontal planes represent the barriers, $\lambda^{-1}$ is a parameter of the family.

“phantom” solutions associated with large magnitudes of this parameter.

In Fig. 7, the “phantom” solutions correspond to those $\lambda^{-1}$ at which solutions, shown by the meshed

## References

[1] V.I. Babitsky, Theory of Vibroimpact Systems and Applications, Springer, Berlin, 1998.

[2] M.F. Dimentberg, Statistical Dynamics of Nonlinear and Time-Varying Systems, Wiley, New York, 1988.

[3] F. Peterka, Introduction to Oscillations of Mechanical Systems with Internal Impacts, Academia Prague, 1981 (in Czech).

[4] S.W. Shaw, R.H. Rand, The transition to chaos in a simple mechanical system, Int. J. Nonlinear Mech. 24 (1) (1989) 41-56.

[5] S. Lenci, G. Rega, Controlling nonlinear dynamics in a two-well impact system. Parts I and II, Int. J. Bifurc. Chaos 8 (12) (1998) 2387-2424.

[6] A.F. Vakakis, L.I. Manevitch, Yu.V. Mikhlin, V.N. Pilipchuk, A.A. Zevin, Normal Modes and Localization in Non-linear Systems, Wiley-Interscience, New York, 1996.

[7] C. Jaffe, P. Brumer, Local and normal modes: a classical perspective, J. Chem. Phys. 73 (11) (1980) 5646-5658.

[8] L.I. Manevitch, Ju.V. Mikhlin, V.N. Pilipchuk, Method of Normal Vibrations for Essentially Non-linear Systems Moscow, Nauka, 1989 (in Russian).

[9] E.G. Vedenova, L.I. Manevitch, V.N. Pilipchuk, Normal oscillations of a string with concentrated masses on nonlinearly supports, Prikl. Mat. Mekh. 49 (2) (1985) 572-578.

[10] V.F. Zhuravlev, A method for analyzing vibration-impact systems by means of special functions, Izv. Akad. Nauk SSSR Mekh. Tverd. Tela (Mech. Solids) 11 (2) (1976) 30-34.

[11] Yu. Mikhlin, A.F. Vakakis, G. Salenger, Direct and inverse problems encountered in vibro-impact oscillations of a discrete system, J. Sound Vib. 216 (2) (1998) 227-250.

[12] A.A. Zevin, Localization of periodic oscillations in vibroimpact systems, XXXV Symposium “Modeling in Mechanics” Gliwice, Poland, Vol. 2, 1996, pp. 261-266.

[13] L.I. Manevitch, M.A.F. Azeez, A.F. Vakakis, Exact solutions for a discrete systems undergoing free vibroimpact oscillations, in: V.I. Babitsky (Ed.), Dynamics of Vibro-Impact Systems, 386 EUROMECH Colloquium, Springer, Berlin, 1999.
[14] V.N. Pilipchuk, Application of special non-smooth temporal transformations to linear and nonlinear systems under discontinuous and impulsive excitation, Nonlinear Dyn. 18 (1999) 203–234.

[15] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, 5th edition, Academic Press, New York, 1994.

[16] R. Bellman, Introduction to Matrix Analysis, McGraw-Hill Company, New York, 1960.

[17] A.A. Zevin, On Stability of Periodic Motions with a Presence of Inter-Impacting Masses, Vol. 4, USSR Academy of Sciences, Mashinovedenie, 1968, pp. 28–36.

[18] V.N. Pilipchuk, A transformation of vibrating systems based on a non-smooth periodic pair of functions. Dokl. Akad. Nauk UkrSSR (Ukrainian Acad. Sci. Rep.) A (4) (1988) 37–40 (in Russian).

[19] V.N. Pilipchuk, Non-smooth spatio-temporal transformation for impulsively forced oscillators with rigid barriers, J. Sound Vib. (2000), in press.

[20] A.F. Filippov, Differential Equations with Discontinuous Right-hand Sides, Kluwer Academic Publishers, Dordrecht, 1988.