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

HAL
open science

# Direct and inverse problems encountered in vibro-impact oscillations of a discrete system

Yuri Mikhlin, Alexander Vakakis, Gary Salenger

► To cite this version:

Yuri Mikhlin, Alexander Vakakis, Gary Salenger. Direct and inverse problems encountered in vibro-impact oscillations of a discrete system. Journal of Sound and Vibration, 1998, 216 (2), pp.227-250. 10.1006/jsvi.1998.1613. hal-01352131

HAL Id: hal-01352131

https://hal.science/hal-01352131

Submitted on 5 Aug 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
DIRECT AND INVERSE PROBLEMS ENCOUNTERED IN VIBRO-IMPACT
OSCILLATIONS OF A DISCRETE SYSTEM

Y. V. MIKHLIN

*Department of Applied Mathematics, Kharkov Polytechnic University, Kharkov 310002, Ukraine*

AND

A. F. VAKAKIS AND G. SALENGER

*Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign, 1206 W. Green Street, Urbana, IL 61801, U.S.A.*

We study direct and inverse problems that arise in the vibro-impact oscillations of a discrete system. Specifically, we examine a class of systems with two coordinates undergoing single- or double-sided impacts; however, the presented techniques are sufficiently general to apply to systems with multiple impacts. The analytical methods employed are a nonlinear normal mode (NNM)-type analysis and a boundary value problem (BVP) formulation, and enable the computation of various branches of bifurcating periodic solutions with different impacting characteristics. Additional insight on the dynamics of these systems is obtained by direct integrations of the equations of motion and by numerical Poincaré maps. It is found that the vibro-impact systems considered possess rich nonlinear dynamics, including vibro-impact localized and nonlocalized time-periodic motions, complicated bifurcation structures giving rise to new types of single- and double-sided impacting motions, mode instabilities, and chaotic responses. We also formulate inverse vibro-impact problems, whereby, we seek the class of dynamical systems that produce specified orbits in the configuration plane. The solutions of the inverse problems are generally non-unique, since they can be reduced to underdetermined sets of algebraic equations with multiple infinities of unknowns. Numerical applications are provided to demonstrate the techniques and validate the analytical results.

# 1. INTRODUCTION

Vibro-impact oscillations are of considerable practical importance and occur in many
engineering applications, including systems with clearances, loose joints, or motion
confining rigid constraints. Vibro-impacts lead to strong and essential nonlinearities and
their analytical or even numerical study presents challenging technical difficulties [1-3].

Among various studies in the literature, Masri and Caughey [4] analyzed linear systems
with rigid constraints by matching linear solutions computed before and after the time
instants of impacts. Studies of piecewise linear and vibro-impact oscillations with
analytical/numerical Poincaré maps and geometrical techniques were performed in a series
of works by Shaw and Holmes [5], Moon and Shaw [2], Shaw [6], Shaw and Rand [3],
and Shaw and Shaw [7]. Ivanov [8] studied vibro-impact oscillations by introducing
auxiliary phase planes, and Zhuravlev [9] analyzed these oscillations by employing special
nonsmooth transformations of variables. In a related work, Vedenova et al. [10] studied
localized and nonlocalized vibro-impact periodic solutions of discrete oscillators using nonsmooth coordinate transformations (cf. also reference [11]). Van de Vorst et al. [12], studied periodic motions of a beam undergoing vibro-impacts by using finite elements and reducing the problem to numerically solving two-point boundary value problems. Numerical and experimental investigations of flexible elastic and rotordynamic systems with impact nonlinearities were carried out in references [13, 14].

Additionally [5, 15–17], local and global bifurcations of periodic orbits and chaotic motions in systems with vibro-impacts were examined. Chaotic motions of such systems were studied in these and other works [18–21]. Whiston [22] presented a study of the “shredding” of the invariant manifolds of a periodically excited system with impacts due to the non-differentiable nature of the vibro-impact dynamics, and Brogliato [23] discussed the dynamics and controls of systems with impact nonlinearities.

In this work we study direct and inverse problems related to vibro-impact oscillations of discrete linear or nonlinear systems with rigid constraints. The theoretical part of the work centers on the use of functional relations between coordinates to study periodic solutions and local dynamics, as well as, the formulation of two-point boundary value problems (BVPs) to study analytic continuations and bifurcations of solutions. The numerical part of the work focuses on the numerical solution of the formulated BVPs and on the construction of Poincaré maps to study the global dynamics. An inverse problem is also posed, whose solution leads to the determination of linear or nonlinear systems that are capable of producing a specified vibro-impact motion in configuration space. Bifurcating problems associated with the inverse problem are discussed and numerical simulations that validate the theoretical findings are presented.

## 2. THE DIRECT PROBLEM

### 2.1. ANALYTICAL RESULTS

We analyze the dynamics of the system depicted in Figure 1(a). It possesses two degrees of freedom, but the methodology to be developed can be extended to the multi-degree-of-freedom case. The motion of each particle is restricted by two rigid barriers, and takes place in the interval $[-e, e]$. As a result, for sufficiently high energy of motion impacts occur and the equations of motion are expressed in the form

$$ \ddot{x} + \frac{\partial \Pi(x, y)}{\partial x} + P(x) = 0, \quad \ddot{y} + \frac{\partial \Pi(x, y)}{\partial y} + P(y) = 0, \quad (1) $$

where $\Pi(x, y)$ denotes the (smooth) potential energy, and $P(u)$ models purely elastic impacts of the two particles with their rigid constraints. The elastic impacts are mathematically modelled by the relationship [cf. Figure 1(b)]

$$ P(u) = \lim_{n \to \infty} c(n) \left(\frac{u}{e}\right)^{2n-1}, \quad (2) $$

where $c(n)$ is a scaling constant depending on $n$.

In the direct problem, we seek analytic approximations and numerical solutions to vibro-impact, time-periodic orbits of the system, with a specified number of impacts per cycle. Whereas this problem has been addressed in the past by other authors (especially for the case of single degree of freedom oscillators), the technique that we propose is new since it relies on the concept of a ‘nonlinear normal mode’ (NNM) [11] to formulate a functional relation between the two coordinates of the system during the vibro-impact
<figure><img src="image_3.png" /><figcaption>Figure 1. The vibro-impact oscillator: (a) system configuration, (b) the limiting process defining the elastic impact forces as $n \to \infty$.</figcaption></figure>

motion. Hence, we assume that during the vibro-impact oscillation the following expression is satisfied at all time instants:

$$y(t) = \hat{y}[x(t)], \quad t \in \mathbb{R}. \qquad (3)$$

The function $y[\cdot]$ is determined by substituting equation (3) into equation (1) and eliminating the time variable from the problem. This leads us to the following quasi-linear ordinary differential equation governing $\hat{y}[\cdot]$ in open intervals of $x$ and $y$ between consecutive impacts,

$$2\hat{y}'' \frac{h - \Pi(x, \hat{y})}{1 + \hat{y}^2} + \hat{y}' \left[ -\frac{\partial \Pi(x, \hat{y})}{\partial x} \right] = -\frac{\partial \Pi(x, \hat{y})}{\partial y}, \quad |x| < e, \quad |\hat{y}(x)| < e, \quad (4)$$

where $h$ denotes the (conserved) energy of the system, and primes denote differentiation with respect to $x$. Assuming that $h - \Pi(x, \hat{y}) > 0$ during the motion, i.e. that the vibro-impacts occur before the system reaches its maximum potential energy value, the equation above is nonsingular, and must be complemented by appropriate boundary conditions at the elastic impacts.

We suppose at this point that the periodic motion under consideration involves only impacts in the $x$-coordinate, and is composed of two branches $y = \hat{y}^1(x)$ and $y = \hat{y}^2(x)$, as depicted in Figure 2(a). Between impacts the functions $\hat{y}^{1,2}[\cdot]$ are governed by equation (4), whereas at each instant of impact they satisfy the compatibility conditions

$$\hat{y}^1(\pm e) = \hat{y}^2(\pm e), \quad \hat{y}^1(\pm e) = -\hat{y}^2(\pm e). \qquad (5)$$

That is, at the time instant $t = t^*$ when the elastic impact occurs, there is $C^0$ and $C^1$ continuity in $y(t)$, but only $C^0$ continuity in $x(t)$. The $C^1$ discontinuity in $x(t)$ satisfies the
<figure><img src="image_4.png" /><figcaption>Figure 2. Vibro-impact motion with impacts only in the x-coordinate: (a) motion in the configuration plane, (b) joining of branches at the point of elastic impact.</figcaption></figure>

jump condition, $\dot{x}(t^* - ) = -\dot{x}(t^* + )$, which in turn leads to the conditions (5). Hence, the direct vibro-impact problem consists of solving two nonlinear ordinary differential equations governing the two branches $\hat{y}^{1,2}[\cdot]$ [each similar to expression (4)] subject to the compatibility conditions (5). Since no exact analytical solutions exist for the general nonlinear problem, it will be necessary to resort to asymptotic approximations. We note that more complicated vibro-impact oscillations corresponding to an arbitrarily large number of impacts in the x- and/or y-coordinates can also exist, but these will not be considered herein. For such motions, a similar asymptotic analysis can be performed by considering vibro-impact motions possessing more than two branches and/or multiple impacts per period in the x- and y-coordinates.

Before we consider the asymptotic solution of the direct problem, we remark that the previous discussion regarding the vibro-impact motion is valid only if the condition $h - \Pi(x, \hat{y}) > 0$ is satisfied. If the system reaches at any time instant its maximum equipotential surface [corresponding to the condition $\Pi(x, \hat{y}) = h$], the coefficient of the highest derivative in equation (4) vanishes and the problem becomes singular. Such problems with singularities at the maximum equipotential surface are encountered in calculations of nonlinear normal modes (NNMs) of discrete oscillators [11]; their solution requires imposing the following boundary orthogonality condition that is valid when $\Pi(x, \hat{y}) = h$:

$$ \left\{ \hat{y}' \left[ -\frac{\partial \Pi(x, \hat{y})}{\partial x} \right] + \frac{\partial \Pi(x, \hat{y})}{\partial y} \right\}_{\Pi(x, \hat{y}) = h} = 0. \quad (6) $$

This condition must be imposed whenever the orbit of the vibro-impact oscillator reaches the maximum equipotential surface.

We now consider the analytic approximation of the two-branch time-periodic orbit depicted in Figure 2(a). It should be clear that such a vibro-impact response can only occur for sufficiently large values of the total energy $h$. Indeed, for sufficiently small $h$ the two particles cannot reach their rigid constraints and the sought time-periodic motions are
merely NNMs of the system with no impacts [11]. To be able to perform analytic computations we assume the following specific form for the potential energy:

$$ \Pi(x, y) = \frac{\alpha}{2}(x^2 + y^2) + \frac{\delta}{3}(x^3 + y^3) + \frac{\beta}{4}(x^4 + y^4) + \frac{\gamma}{2}(x - y)^2. \quad (7) $$

In the above expression, the scalars $\alpha, \beta, \delta$ are the coefficients of the quadratic, quartic and cubic parts of the potential functions of the oscillators, whereas $\gamma$ is the coefficient of the weak linear coupling stiffness.

At low energies and sufficiently small coupling $\gamma$, this system possesses two localized NNMs, corresponding to periodic oscillations of the two masses with amplitudes of $O(1)$ and $O(\gamma)$, respectively [27]. We are interested in examining the behavior of these localized NNMs as the energy increases and impacts start to occur. Hence, in what follows we will assume that all coefficients in equation (7) are of $O(1)$ with the exception of the coupling coefficient which will be assumed small, $\gamma \ll 1$. Considering the vibro-impact orbit of Figure 2(a), the two branches $\hat{y}^{1,2}(x)$ must satisfy the relations (4) and (5). For a potential function given by equation (7), and taking into account that the problem possesses a small (perturbation) parameter $\gamma$, we express the two branches in the following series of successive approximations:

$$ \hat{y}^j(x) = \sum_{k=1} \gamma^k \hat{y}_k^j(x), \quad j = 1, 2. \quad (8) $$

With equation (8) we seek vibro-impact motions that localize to the particle with coordinate $x$ [note that with $x = O(1)$, we seek $y = O(\gamma) \ll x$. Substituting equation (8) into the equations governing $\hat{y}^{1,2}(x)$, we obtain

$$ \begin{aligned} & 2[\gamma\hat{y}_1^j + \gamma^2\hat{y}_2^j + O(\gamma^2)]\left[h - \frac{\alpha x^2}{2} - \frac{\delta x^3}{3} - \frac{\beta x^4}{4} - \frac{\gamma}{2}(x^2 - 2\gamma x \hat{y}_1^j) + O(\gamma^2)\right][1 + O(\gamma^2)] \\ & \quad + [\gamma\hat{y}_1^j + \gamma^2\hat{y}_2^j + O(\gamma^3)][-\alpha x - \delta x^2 - \beta x^3 - \gamma(x - \gamma\hat{y}_1^j) + O(\gamma^3)] \\ & = -\alpha\gamma\hat{y}_1^j - \alpha\gamma^2\hat{y}_2^j - \delta\gamma^2(\hat{y}_1^j)^2 - \gamma(\gamma\hat{y}_1^j - x) + O(\gamma^3) \end{aligned} \quad (9) $$

for $j = 1, 2$. This expression is valid for $|x| < e$, i.e. between consecutive impacts, and for sufficiently large values of $h$ so that the condition $h - \Pi(x, \hat{y}) > 0$ is satisfied pointwise for each of the two branches of the vibro-impact motion. When $|x| = e$ the elastic impact conditions (5) are imposed and matching of the two branches takes place.

Considering terms in equation (9) proportional to $\gamma$ we obtain the following quasi-linear ordinary differential equation governing $\hat{y}_1^j$:

$$ 2\hat{y}_1^j \left[ h - \frac{\alpha x^2}{2} - \frac{\delta x^3}{3} - \frac{\beta x^4}{4} \right] + \hat{y}_1^j (-\alpha x - \delta x^2 - \beta x^3) + \alpha \hat{y}_1^j = x, \quad |x| < e, \quad j = 1, 2. \quad (10) $$

Since by assumption the coefficient of the leading derivative is a nonzero quantity, we can express the solution in the power series form

$$ \hat{y}_1^j = \sum_{p=0}^{\infty} c_{jp} x^p, \quad (11) $$
with the coefficients of the series determined by substituting equation (11) into equation (10) and matching corresponding powers of $x$. Following this procedure we obtain the approximation

$$ 
\begin{aligned}
\hat{y}_1' &= c_{j0} \left[ 1 - \frac{\alpha x^2}{4h} - \frac{\alpha^2 x^4}{32h^2} - \frac{\alpha \delta x^5}{48h^2} \right] + c_{j1} \left[ x + \frac{\delta x^4}{24h} + \frac{\beta x^5}{40h} \right] + \left[ \frac{x^3}{12h} + \frac{\alpha x^5}{60h^2} \right] + O(x^6) \\
&\equiv c_{j0} f_1(x, h) + c_{j1} f_2(x, h) + f_3(x, h) + O(x^6), \quad j = 1, 2,
\end{aligned} 
\quad (12) $$

where the leading coefficients of the series remain undetermined up to this point. A similar expression for the second (or higher) order approximations $\hat{y}_2'$ can be derived in the form, $\hat{y}_2' \equiv d_{j0}g_1(x, h) + d_{j1}g_2(x, h) + g_3(x, h) + O(x^6)$, but this higher order correction will not be pursued in the present work.

Imposing the elastic impact conditions (5) and taking into account only terms of $O(\gamma)$ [expression (12)], we obtain the following matrix equation governing the leading order coefficients $c_{jp}$, $j = 1, 2, p = 0, 1$:

$$ 
\begin{bmatrix} f_1(e, h) & f_2(e, h) & -f_1(e, h) & -f_2(e, h) \\ f'_1(e, h) & f'_2(e, h) & f'_1(e, h) & f'_2(e, h) \\ f_1(-e, h) & f_2(-e, h) & -f_1(-e, h) & -f_2(-e, h) \\ f'_1(-e, h) & f'_2(-e, h) & f'_1(-e, h) & f'_2(-e, h) \end{bmatrix} \begin{Bmatrix} c_{10} \\ c_{11} \\ c_{20} \\ c_{21} \end{Bmatrix} = \begin{Bmatrix} 0 \\ -2f_3(e, h) \\ 0 \\ -2f_3(-e, h) \end{Bmatrix}. \quad (13) $$

Provided that the matrix of the determinant of coefficients is nonsingular, the above equation yields a unique solution for the coefficients, which is derived symbolically using *Mathematica*, as follows:

$$ 
\begin{aligned}
c_{10} &= c_{20} = \frac{8\delta e^4 (\alpha e^2 + 3h)}{\alpha(-9\alpha\beta e^6 + 10\delta^2 e^6 - 72\alpha e^2 h - 36\beta e^4 h - 288h^2)}, \\
c_{11} &= c_{21} = \frac{6e^2 (\alpha^2 e^4 + 7\alpha e^2 h + 12h^2)}{h(-9\alpha\beta e^6 + 10\delta^2 e^6 - 72\alpha e^2 h - 36\beta e^4 h - 288h^2)}.
\end{aligned} 
\quad (14) $$

This solution shows that leading order approximations for the two branches coincide, $\hat{y}_1' = \hat{y}_2'$, and, furthermore, that the slope of the vibro-impact trajectory in the configuration plane at the points of impact is equal to zero, $\hat{y}_1'(\pm e) = \hat{y}_2'(\pm e) = 0$. In Figure 3 we depict the asymptotic approximation (11)–(14) for parameter values $\alpha = \beta = \delta = 1, e = 1, \gamma = 0.1$, and $h = 2$. In the same figure we show the maximum equipotential energy level corresponding to $h = 2$, which is sufficiently wide to allow double-sided vibro-impact motions to occur.

Numerical simulations were also performed, integrating the equations of motion (1) with potential energy given by equation (7). The numerical results validated the analytical predictions, indicating that the double-sided vibro-impact localized motion is orbitally stable, and, hence, physically realizable. Moreover, using the numerical simulations it was established that the localized motion persists as a stable solution for arbitrarily high levels of the energy $h$, and that no bifurcations associated with this motion take place in the system for large $h$. The above theoretical and numerical results prove the existence of localized vibro-impact oscillations, where the energy of the motion is mainly confined to one of the two particles. Moreover, the motion in the configuration plane is represented by a line (and not by a closed loop), which indicates synchronous motion of the two particles (purely in-phase or out-of-phase). We will show later that this is not the case for systems
<figure><img src="image_5.png" /><figcaption>Figure 3. First order asymptotic approximation $\hat{y}_1^1(x) \equiv \hat{y}_1^2(x)$ to a localized double-sided vibro-impact motion with parameters $\alpha = \beta = \delta = 1$, $e = 1$, $\gamma = 0.1$ and $h = 2$.</figcaption></figure>

with other types of potential function, which can support periodic motions possessing nontrivial phase differences between particles.

An interesting question concerns the behavior of the previous localized solution as $h$ decreases, and we now briefly discuss this issue. For the system with (asymmetric) potential function (7), we expect that as the energy decreases below a critical value the particle characterized by the coordinate $x$ will make the transition from double-sided to single-sided impacts, and (after further decrease of energy), eventually, to no impacting motions. Indeed, assuming small $\gamma$ and a localized periodic motion of the form $\hat{y}(x) = \gamma \hat{y}_1(x) + O(\gamma^2)$, the function $\hat{y}[\cdot]$ is still governed by equation (4). The first critical value of the energy $h_{c1}$ separating the regimes of single- and double-sided impacts is determined approximately by the condition

$$ \Pi(e, \hat{y}(e)) = \frac{\alpha e^2}{2} + \frac{\delta e^3}{3} + \frac{\beta e^4}{4} + \frac{\gamma e^2}{2} + O(\gamma^2) = h_{c1}, \quad (15) $$

which corresponds to the energy level for which the particle characterized by $x$ reaches its maximum amplitude when $x = e$, and, hence, the relation $h_{c1} - \Pi(e, \hat{y}) = 0$ is satisfied. Since the potential energy of the system is asymmetric, at $x = -e$ the condition $h_{c1} - \Pi(-e, \hat{y}) > 0$ is still satisfied and impact of the particle with its rigid constraint occurs. For $\alpha = \beta = \delta = 1$, $e = 1$ and $\gamma = 0.1$ we find that $h_{c1} = 1.1333$. As a result, for a sufficiently small neighborhood of $h$ bounded from above by $h_{c1}$ the system undergoes single-sided localized vibro-impact oscillations, and the motion in the configuration plane $\hat{y}(x) = \gamma \hat{y}_1(x) + O(\gamma^2)$ is determined by solving equation (4) subject to the boundary conditions

$$ \hat{y}'(-e) = 0, \quad (\text{elastic impact at } x = -e), \quad (16a) $$

$$ \left\{ \hat{y}' \left[ -\frac{\partial \Pi(x, \hat{y})}{\partial x} \right] + \frac{\partial \Pi(x, \hat{y})}{\partial y} \right\}_{x=x_1} = 0, \quad (\text{right boundary of the motion}), \quad (16b) $$

where $X_1$ is the maximum amplitude attained by $x$ at the right boundary of the single-sided vibro-impact motion (cf. Figure 4). Performing a perturbation analysis similar to the
<figure><img src="image_6.png" /><figcaption>Figure 4. First order asymptotic approximation $\hat{y}_1(x)$ to the localized single-sided vibro-impact motion with parameters $\alpha = \beta = \delta = 1$, $e = 1$, $\gamma = 0.1$ and $h = 0.933333$.</figcaption></figure>

double-sided impact case, we express $\hat{y}_1(x)$ as, $\hat{y}_1(x) = c_0f_1(x, h) + c_1f_2(x, h) + f_3(x, h) + O(x^6)$, where the functions $f_k(x, h)$, $k = 1, 2, 3$, are defined by equation (12). Imposing the conditions (16) we then determine the values of the coefficients $c_0$ and $c_1$. For the above values of the system parameters and $h = 0.93333 < h_{c1}$, we find $c_0 = 0.20042$ and $c_1 = -0.49791$. In Figure 4 the motion in the configuration plane is presented, superimposed to the maximum equipotential energy level. The left boundary curve possesses zero slope since it encounters an impact, whereas the right boundary terminates at the maximum equipotential energy level; at that point the system possesses zero kinetic and purely potential energy. The localized motion of Figure 4 is the analytic continuation of the double-sided vibro-impact motion discussed earlier when the energy is decreased a small distance below the critical value $h_{c1}$.

With a further decrease of the energy below a second critical value $h_{c2}$, the system ceases to impact with its rigid barriers and performs localized NNM oscillations (Vakakis et al. [11]). The value $h_{c2}$ is determined from the relation

$$ \Pi(-e, \hat{y}(-e)) = \frac{\alpha e^2}{2} - \frac{\delta e^3}{3} + \frac{\beta e^4}{4} + \frac{\gamma e^2}{2} + O(\gamma^2) = h_{c2}, \quad (17) $$

which for the above values of the system parameters is equal to 0.46666. For $h < h_{c2}$ the localized NNM is computed by solving equation (4) subject to the boundary conditions (6). Expressing the mode as $\hat{y}(x) = \gamma\hat{y}_1(x) + O(\gamma^2)$, a perturbation analysis can be used to compute the first order approximation. The resulting analytical approximation for $h = 0.26666 < h_{c2}$ is depicted in Figure 5, superimposed to the maximum equipotential energy level.

The previous analysis proves the existence of a localized vibro-impact motion. As the energy decreases this localized motion is preserved in the single-sided impact regime, and ultimately degenerates to a localized NNM of the (smooth) dynamical system with no impacts. Hence, the localized vibro-impact periodic motions discussed in this section can be viewed as analytical continuations of strongly localized NNMs that are encountered in many classes of smooth dynamical systems [11]. As mentioned earlier, the localized vibro-impact motion is orbitaly stable at all levels of energy. This result, along with
<figure><img src="image_7.png" /><figcaption>Figure 5. First order asymptotic approximation $\hat{y}_1(x)$ to the localized nonlinear normal mode of the system with parameters $\alpha = \beta = \delta = 1$, $e = 1$, $\gamma = 0.1$ and $h = 0.26666$.</figcaption></figure>

additional ones concerning bifurcations of other (nonlocalized) types of vibro-impact motions, was deduced by means of direct numerical simulations and construction of Poincaré maps.

## 2.2. NUMERICAL RESULTS

The analytical results of the previous section were *local* in nature, establishing the existence of localized vibro-impact motions and providing asymptotic approximations for sufficiently small values of the coupling parameter $\gamma$. The *global* vibro-impact response of the system will now be considered by numerically integrating the equations of motion (1) and constructing numerical Poincaré maps [24–26].

To construct the Poincaré maps we note that the motion of equation (1) generally takes place in the four dimensional phase space $(x, \dot{x}, y, \dot{y})$. Since the impacts are assumed to be purely elastic the energy $h$ is conserved throughout the oscillation, and by fixing it to a specific level we can restrict the flow of the dynamical system to an isoen energetic three-dimensional space. If, in addition, we intersect transversely the three-dimensional isoen energetic flow by a two-dimensional cut-plane, the cut-plane defines a two-dimensional Poincaré map which can be used to study the global dynamics of the system. Choosing the cut plane as $T: \{x = 0\}$, the Poincaré map $P(\cdot)$ is defined as

$$P: \Sigma \to \Sigma, (y, \dot{y}) \to P(y, \dot{y}),$$

where the Poincaré section is given by $\Sigma = \{x = 0, \dot{x} > 0\} \cap \{\text{energy } h\}$. The additional restriction of positive $\dot{x}$ at the point of intersection was imposed in order that the Poincaré map be orientation preserving [25]. A period-one vibro-impact orbit pierces the cut-section $\Sigma$ only once and corresponds to a fixed point of the map. The stability of the periodic orbit can be determined by examining near-by orbits of the Poincaré map: if the fixed point appears as a center and is surrounded by closed curves, the corresponding orbit is stable; otherwise it is unstable. Subharmonic vibro-impact orbits of higher periods pierce $\Sigma$ in more than one point, whereas quasiperiodic orbits are represented by closed loops. Chaotic motions appear as randomly distributed points that occupy a definite region of the map.
The numerical simulations were carried out by integrating equations (1) between impacts using the fifth-order Runge–Kutta SLATEC routine DDERKF. The energy-conserving elastic impacts were modelled by imposing new initial conditions after each impact, to reflect the abrupt sign change in velocity after the x- or y-particle collide with their rigid constraint. The exact time instant when the impact occurs is of great importance in this numerical scheme, since it determines the points where the new initial conditions must be imposed. The time instants of impact were determined numerically by applying Newton’s method to the numerical results of the integration close to each instance of impact. Special care was taken to choose a sufficiently small time steps in order to ensure accurate convergence of Newton’s method and precise calculation of the instant of impact. This is

<figure><img src="image_8.png" /><figcaption>Figure 6. Poincaré map of the neighborhood of the localized vibro-impact motion of a system with $\alpha = \beta = \delta = 1$, $e = 1$, $\gamma = 0.1$, and energy (a) $h = 18$, and (b) $h = 99$.</figcaption></figure>
particularly important for grazing (near impact) oscillations where the accuracy of the numerical integrations is important in order to construct correct bifurcation diagrams and study mode instabilities initiated by impact.

### 2.2.1. Nonlinear potential function

In Figures 6(a, b) we present the neighborhoods of the Poincaré maps close to the localized double-sided vibro-impact motions of equation (1) with potential energy (7), parameters $α = β = δ = 1$, $e = 1$, $γ = 0.1$, and energy $h = 18$ and $99$. The Poincaré maps consist of centers surrounded by closed curves, indicating that the vibro-impact localized motions are stable even at high values of the energy. The global dynamics of the system are more clearly discerned from the Poincaré maps of Figure 7. In Figure 7(a) we depict the global map of the system for energy $h = 2$. We note the large central region that is dominated by the localized vibro-impact periodic motion. In addition, the system possesses a stable in-phase periodic orbit represented by the center at the top of the map; this motion corresponds to in-phase synchronous oscillations of the two particles satisfying $x = y$. We also note the large chaotic ring, composed of chaotic vibro-impact motions, that surrounds the central region. This ‘stochastic sea’ possess islands of regular subharmonic motions. For comparison purposes, in Figure 7(b) we present the Poincaré map of the same system but with no impacts. The system with no impacts possess a stable localized NNM that also dominates the central part of the map, a smaller chaotic region, and a stable in-phase NNM at the top of the map. Islands of subharmonic motion can also be detected in this case.

### 2.2.2. Linear potential function

To perform a systematic study of the effects of the vibro-impacts on the global dynamics of system (1) we performed a detailed computation of Poincaré maps for the case of the quadratic potential,

$$ \Pi(x, y) = \frac{\alpha}{2}(x^2 + y^2) + \frac{\gamma}{2}(x - y)^2, \qquad (18) $$

with parameter values and varying values of the clearance $e$. In Figure 8(a–j) we present the results of the numerical simulations, ranging from a system with no vibro-impacts [Figure 8(a)] to a system with vibro-impacts and clearance $e = 1$.

We note that as vibro-impacts start to occur, a region of chaotic motions starts to develop at the boundaries of the domains of influence of the stable in-phase and out-of-phase modes [cf. Figures 8(a–e)]. This central chaotic region expands with decreasing clearance and surrounds the domains of influence of the two modes, where regular quasiperiodic motions continue to occur. For $e \approx 1.78$ two saddle-node bifurcations occur and two pairs of stable-unstable vibro-impact modes are generated. In the Poincaré map of Figure 8(f) only the one pair of newly created modes can be observed, since the other pair is very close to the lower boundary curve of the map, and is difficult to observe. The small chaotic region surrounding one of the newly generated stable vibro-impact mode close to the origin of the axes is due to transverse homoclinic intersections of the invariant manifolds of one of the now unstable out-of-phase mode, and, hence, is a direct result of the mode bifurcation. A similar bifurcation of nonlinear normal modes was previously observed in a system with smooth nonlinearities [26].

With the clearance further reduced, the newly created mode becomes localized (it approaches the origin of the axes), and its domain of influence expands. At the same time, the regular regions surrounding the in-phase and out-of-phase vibro-impact modes contract [cf. Figure 8(g)], until the in-phase mode becomes unstable close to $e = 1.4$.
<figure><img src="image_9.png" /><figcaption>Figure 7. Global Poincaré map of the system with $\alpha = \beta = \delta = 1$, $\gamma = 0.1$ and energy $h = 2$, when (a) vibro-impacts occur at $e = 1$, and (b) no vibro-impacts exist.</figcaption></figure>

[cf. Figure 8(b)] and the out-of-phase mode becomes unstable below $e = \sqrt{5/3}$. We note that the in-phase mode regains stability for $e \le 1.4$; the bifurcations associated with the small instability region of this mode is discussed in detail below. The out-of-phase mode loses stability in a pitchfork bifurcation and remains unstable for smaller values of $e$. This behavior of the out-of-phase mode is consistent with the findings of reference [26] for a similar system with smooth nonlinearities. Further decreases of the clearance results in an expansion of the regular domain of influence of the localized modes, with the in-phase mode continuing to be stable [cf. Figure 8(i, j)]. As a general observation we note that the region of chaotic motions contracts as the clearance is further reduced and the stable localized vibro-impact modes dominate the global dynamics of the system. In the following
<figure><img src="image_10.png" /><figcaption>Figure 8—(Caption on next page)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8. Global Poincaré map of the system with $\alpha = 1, \beta = \delta = 0, \gamma = 0.1$ and energy $h = 2$, for clearances: (a) $e = 2.0$ (no impacts), (b) $e = 1.9$, (c) $e = 1.88$, (d) $e = 1.85$, (e) $e = 1.8$, (f) $e = 1.78$, (g) $e = 1.5$, (h) $e = 1.4$, (i) $e = 1.2$, (j) $e = 1.0$.</figcaption></figure>

discussion we examine more closely the separate sequences of bifurcations associated with the instability of the in-phase mode, and the generation of the localized mode in the vibro-impact system.

We start with the bifurcations giving rise to the instability of the in-plane vibro-impact mode. These are complicated and generate additional periodic vibro-impact motions that have no counterparts in the corresponding smooth dynamical systems with no impacts. To investigate the series of bifurcations associated with the in-phase mode instability we performed detailed Poincaré map calculations in the neighborhood of that mode, as well as, analytical/numerical boundary value problem (BVP) computations to precisely determine the bifurcation points. We start with a discussion of the Poincaré maps in the vicinity of the in-phase mode, as depicted in Figure 9(a–f).

In Figure 9(a) (corresponding to $e = 1.418$) the in-phase mode (labelled ‘A’) is orbitally stable and non-impacting. That is, for the level of energy $h = 2$ the in-phase synchronous oscillations of the two particles do not possess enough amplitude to reach the rigid barriers, and are, in fact, in-phase NNMs of the system. We also note the large region of chaotic vibro-impact motions surrounding the in-phase mode, and the plethora of co-existing regular quasi-periodic motions. As $e$ is reduced to $1.417$ [cf. Figure 9(b)], the in-phase motion ‘A’ is still stable and non-impacting, and an additional stable periodic motion appears (labelled as ‘B’) on the right of the in-phase mode. As shown below, this new orbit corresponds to a double-impacting synchronous motion where both x- and y-coordinates undergo single-sided impacts; clearly this type of orbit cannot be realized in the corresponding smooth dynamical system. Moreover, it will be shown that this new motion is generated through a Saddle-node bifurcation, and that there exists an additional unstable orbit that cannot be discerned in the Poincaré map; also, due to symmetry, there exists an additional stable-unstable pair of similar orbits lying on the left side of the in-phase mode, which is not depicted in Figure 9(b). At $e = \sqrt{2}$ [Figure 9(c)] the in-phase mode ‘A’ starts impacting and loses stability, whereas, the stable motion ‘B’ expands its surrounding domain of influence which is composed of regular quasi-periodic orbits. It turns out that at least four unstable branches of vibro-impacting periodic motions converge to the in-phase mode at its point of instability, including the two unstable counterparts of motion ‘A’. As discussed below, the in-phase impacting mode remains unstable until $e \approx 1.40099$, when it coalesces with two new unstable branches of vibro-impact motions and regains orbital stability. This can be realized from Figures 9(e) for $e = 1.398$, where both the in-phase impacting ‘A’ and the periodic motion ‘B’ are stable.
<figure><img src="image_12.png" /><figcaption>Figure 9. Poincaré maps in the neighborhood of the in-phase mode of the system with $\alpha = 1, \beta = \delta = 0, \gamma = 0.1$ and energy $h = 2$, for clearances: (a) $e = 1.418$, (b) $e = 1.417$, (c) $e = \sqrt{2}$, (d) $e = 1.402$, (e) $e = 1.398$, (f) $e = 1.390$</figcaption></figure>

A further reduction of $e$ to $1.390$ reduces the domain of influence of motion ‘B’ [cf. Figure 9(f)], and at $e \approx 1.38813$ the orbit ‘B’ ceases to exist when it coalesces with another unstable vibro-impact motion in a saddle-node bifurcation. Hence, the small region of instability of the in-phase motion is associated with a complicated series of bifurcations, which we now proceed to examine in detail using BVP formulations.

The BVP formulation is based on the realization that, since the motions of the two particles between instants of impact is linear, they can be explicitly determined as

$$ 
\begin{aligned}
x(t) &= C_1 \cos(\sqrt{\alpha}t - \phi_1) - C_2 \cos(\sqrt{\alpha + 2\gamma}t - \phi_2), & |x| \le e; \\
y(t) &= C_1 \cos(\sqrt{\alpha}t - \phi_1) + C_2 \cos(\sqrt{\alpha + 2\gamma}t - \phi_2), & |y| \le e.
\end{aligned} 
\quad (19) 
$$
The undetermined amplitudes and phases can be evaluated by imposing appropriate boundary conditions at instants of impact for each of the vibro-impact orbits. Then, by construction the solution (19) will be periodic in $t$. The first time-periodic solution to be examined is a double-sided impacting synchronous motion, where both x- and y-coordinates undergo single-sided impacts. For this type of motion the appropriate boundary conditions that must be imposed on equations (19) are

$$ \dot{x}(T/2) = 0, \quad x(0) = \mp e, \quad \dot{y}(0) = 0, \quad y(T/2) = \pm e, \qquad (20) $$

where $T$ is the period of the motion, and the origin of the $t$-axis is chosen at the time instant when the $x$-coordinate impacts. Combining equations (19) and (20), we obtain the set of algebraic equations

$$ \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & \omega_1 & 0 & \omega_2 \\ \cos \omega_1 T & \sin \omega_1 T & \cos \omega_2 T & \sin \omega_2 T \\ -\omega_1 \sin \omega_1 T & \omega_1 \cos \omega_1 T & \omega_2 \sin \omega_2 T & -\omega_2 \cos \omega_2 T \end{bmatrix} \begin{Bmatrix} C_1 \cos \phi_1 \\ C_1 \sin \phi_1 \\ C_2 \cos \phi_2 \\ C_2 \sin \phi_2 \end{Bmatrix} = \begin{Bmatrix} -e \\ 0 \\ e \\ 0 \end{Bmatrix}, \quad (21) $$

where $\omega_1 = \sqrt{\alpha}$ and $\omega_2 = \sqrt{\alpha + 2\gamma}$. Provided that the matrix of coefficients is nonsingular, the above relation provides a unique solution for the coefficients $C_1, C_2, \phi_1$, and $\phi_2$ and determines the solution (19). Once the solution is determined, the period $T$ of the motion is computed by imposing the additional condition that the total energy of the motion is conserved and equal to $h$. This completes the calculation. Singularities of the matrix of coefficients indicates points where analytic continuation of this particular branch of solutions cannot take place, i.e. bifurcation points of the solution.

Similar BVPs can be formulated for determining the branches of solutions and bifurcation points of other types of vibro-impact modes, such as, in-phase synchronous motions [solution (19) with boundary conditions $x(0) = -e$, $x(T/2) = e$, $y(0) = -e$, $y(T/2) = e$], out-of-phase synchronous motions [boundary conditions $x(0) = -e$, $x(T/2) = e$, $y(0) = e$, $y(T/2) = -e$], and ‘loop’ vibro-impact modes (see discussion and results below, corresponding to boundary conditions $x(0) = \mp e$, $\dot{x}(T/2) = 0$, $\dot{y}(0) = 0$, $y(T/2) = \mp e$).

The numerical solutions of the aforementioned BVPs enable us to fully understand the series of bifurcations that generate the instability of the in-phase mode, as depicted in the Poincaré maps of Figure 9. In Figure 10 we provide the complete bifurcation structure that gives rise to this instability for the system whose Poincaré maps were depicted previously. In the plots of this figure we graph the values of $y$ and $\dot{y}$ vs $e$ for each branch of periodic solutions at the time instant when $x = 0$ and $\dot{x} > 0$. There are four families of periodic solutions in this bifurcation structure. For $e \ge \sqrt{2}$ the in-phase mode ‘A’ is non-impacting and stable. At $e = \sqrt{2}$ the in-phase mode becomes where it coalesces with four unstable branches: two of loop’ vibro-impact modes ‘C’, and two of double-sided synchronous modes ‘B’. The mode ‘B’ (which were observed in the Poincaré plots of Figure 9) are generated in two saddle-node bifurcations and exist in stable-unstable pairs up to $e = \sqrt{2}$. At $e = (\sqrt{h}/\omega_1) \sin(\pi\omega_1/2\omega_2) \approx 1.40099$ the in-phase mode regains stability
<figure><img src="image_13.png" /><figcaption>Figure 10. The bifurcation structure associated with the instability of the in-phase mode. Graphs of (a) *y* and (b) *y*' vs *e* for each periodic solution at the time instant when *x* = 0 and *x*' > 0: ——, stable; ----, unstable motions.</figcaption></figure>

in a pitchfork bifurcation and generates two unstable branches of solutions 'D'. These are in-phase synchronous motions with coinciding endpoints with the in-phase mode, and cease to exist at $e = (1/\omega_1)\sqrt{h/\{1 + 2\tan^2[(\pi/2) + (\pi\omega_1/2\omega_2)]\}} \approx 1.38813$, when they coalesce with the stable modes 'B' in saddle-node bifurcations. The bifurcation structure of Figure 10 agrees completely with the numerical simulations of Figure 9, and leads to a full understanding of the dynamic phenomena associated with the instability of the in-phase mode.

The bifurcation results depicted in Figure 10 reveal the richness of the dynamics of the system. Indeed, in the region where the symmetric in-phase NNM loses stability it gives rise to stable or unstable motions with differing impacting patterns. Solutions C (cf. Figure 10) are of particular interest due to their loop-like appearance, indicating nontrivial phase differences between the two position coordinates; moreover, these motions
<figure><img src="image_14.png" /><figcaption>Figure 11. The bifurcation structure associated with the generation of the localized vibro-impact mode. Graphs of y vs e for each periodic solution at the time instant when x = -e: —, stable; —, unstable motions.</figcaption></figure>

correspond to anisochronous and single-impacts for x and y. Similar anisochronous and single-impacts occur in solutions B, but for these motions trivial phase differences occur between the coordinates. Nontrivial phase differences between x and y also occur in solutions D, but these represent unstable free oscillations. We point out that it would be of considerable interest to analytically study the degenerate bifurcations associated with the loss of stability of the in-phase mode at e = √2. These bifurcations indicate that the initiation of impacts (grazing) for the in-phase mode is destabilizing, although the mode regains stability when the non-linear impacting effects become stronger for e ≤ 1·40099. Additional bifurcations are expected to occur involving motions in the configuration plane of increased complexity. Such bifurcations can be studied using the techniques outlined herein, by imposing different boundary conditions and impacting restrictions for the motions to be studied.

Similar BVP computations help us understand the bifurcations giving rise to the localized vibro-impact mode that appears in the Poincaré maps of Figure 8. As mentioned previously, due to symmetry there exists two such localized modes, one localizing in the x-coordinate, and the other in the y-coordinate. Here we will analyze only the mode localizing in the x-coordinate, but a similar calculation can be performed for the other mode. To this end, we seek localized solutions governed by equation (19) between impacts, with only the x-coordinate impacting, and boundary conditions

$$x(0) = -e, \quad x(T/2) = e, \quad \dot{y}(0) = 0, \quad \dot{y}(T/2) = 0. \qquad (22)$$

In Figure 11 we depict the numerical solutions of the BVP. The solutions are depicted as graphs of y vs e for each periodic solution at the time instant when x = -e. We note that the localized mode is generated at e ≈ 1·79 in a saddle-node bifurcation. The unstable mode of the bifurcation has the same impacting properties with the localized mode and ceases to exist at e = √5/3 ≈ 1·291 when it coalesces with the out-of-phase vibro-impact mode in a pitchfork bifurcation. After this bifurcation the out-of-phase modes loses stability, and remains unstable for lower values of e. It is the transverse intersections of the invariant manifolds of this unstable mode that give rise to the small chaotic regime in the central region of the Poincaré plot of Figure 8(f).

In the previous sections we analyzed some of the dynamic phenomena occurring in the direct vibro-impact problem of a two-degree-of-freedom linear or nonlinear oscillator undergoing purely elastic impacts. We considered both analytical and numerical
techniques, which enabled us to study local and global dynamics, and bifurcations generating vibro-impact mode localization and mode instabilities. We now turn our attention to the inverse problem associated with the two degree-of-freedom system. This problem seeks the systems that under vibro-impact conditions generate an a priori specified trajectory in the configuration plane. The solutions of the inverse problem provide interesting insight on the class of dynamical systems that are capable of producing required pre-determined vibro-impact orbits.

### 3. THE INVERSE PROBLEM

The inverse problem associated with the vibro-impact oscillator of Figure 1 is posed as follows. Consider the equations of motion (1) and assume that the trajectory of the system in the configuration plane follows a specified orbit: What is the required potential function $\Pi(x, y)$ required for such an orbit to exist? Moreover, under what conditions does such a potential function exist and when is it unique? In what follows we will discuss in detail a method of solution to the inverse problem and will provide some numerical applications. Although the inverse problem considered herein addresses the two-degree-of-freedom vibro-impact oscillator, it can be easily generalized to higher dimensions. In addition, in the following exposition the inverse problem is formulated for a motion composed of two two-sided impacting branches. Generalization to other types of vibro-impact motion can be similarly performed.

Considering again the equations of motion (1), we assume that the orbit of the system in the configuration space is composed of two double-sided impacting branches. We further assume that each branch is sufficiently smooth to be expanded in power series in $x$ for the entire range of the vibro-impact motion:

$$y = \hat{y}^1(x) = \sum_{k=0} a_k x^k \quad \text{and} \quad y = \hat{y}^2(x) = \sum_{k=0} b_k x^k \quad \text{for} \quad -e \le x \le e. \quad (23)$$

Assuming that the coefficients $a_k$ and $b_k$ in the above series are specified, we seek the potential energy function $\Pi(x, y)$ that is required for such a motion. Expressing the potential function in the general form

$$\Pi(x, y) = \sum_{i+j \ge 2} c_{ij} x^i y^j \quad \text{for} \ i, j = 0, 1, 2, \dots, \quad (24)$$

the solution of the inverse problem reduces to the determination of the coefficients $c_{ij}$ in the above expression.

Since the branches (23) represent double-sided vibro-impact orbits, they must satisfy the differential equation (4) and the set of compatibility conditions (5). Substituting equations (23) and (24) into equations (4) and (5) we obtain the algebraic equations

$$\begin{aligned} & 2\left[\sum_{k=2} k(k-1)\zeta_k x^{k-2}\right]\left[h - \sum_{i+j \ge 2} c_{ij} x^i \left(\sum_{k=0} \zeta_k x^k\right)^j\right] \\ & - \left[1 + \left(\sum_{k=1} k\zeta_k x^{k-1}\right)^2\right]\left(\sum_{k=1} k\zeta_k x^{k-1}\right)\left[\sum_{i+j \ge 2} i c_{ij} x^{i-1} \left(\sum_{k=0} \zeta_k x^k\right)^j\right] \\ & + \left[\sum_{i+j \ge 2} j c_{ij} x^i \left(\sum_{k=0} \zeta_k x^k\right)^{j-1}\right]\left[1 + \left(\sum_{k=1} k\zeta_k x^{k-1}\right)^2\right] = 0, \quad \zeta_k = a_k \text{ or } b_k, \quad (25a) \end{aligned}$$
complemented by the compatibility conditions

$$ \sum_{k=0} a_k(\pm e)^k = \sum_{k=0} b_k(\pm e)^k \quad \text{and} \quad \sum_{k=1} ka_k(\pm e)^{k-1} = -\sum_{k=1} kb_k(\pm e)^{k-1}. \quad (25b) $$

These equations govern the sought coefficients $c_{ij}$, and their solution provides the answer to the inverse problem. We note that conditions (25b) do not involve the unknowns $c_{ij}$, and represent mere geometric restrictions on the assumed vibro-impact branches (23). Considering relations (25a), we make the initial observation that, if we truncate the infinite summations to only the first $N$ terms, by matching coefficients of respective powers $x^k$ we obtain an *underdetermined* set of $2N$ algebraic equations for $N^2$ unknowns. Hence, in general, the inverse problem appears to possess non-unique solutions.

Due to the difficulty of analyzing the general equations (25), it will be necessary to consider specific geometric forms for the branches (23) in order to proceed with concrete solutions of the inverse problem. Hence, as an application, we consider a vibro-impact time periodic solution composed of the two branches

$$ y = \hat{y}^1(x) = a - (a/e^2)x^2 \quad \text{and} \quad y = \hat{y}^2(x) = -\hat{y}^1(x) = -a + (a/e^2)x^2 \quad \text{for} \quad -e \le x \le e. \quad (26) $$

In terms of the previous notation we have, $a_0 = a$, $a_1 = 0$, $a_2 = -(a/e)^2$, $a_k = 0$, $k \ge 3$, and $b_0 = -a$, $b_1 = 0$, $b_2 = (a/e)^2$, $b_k = 0$, $k \ge 3$. Substituting these values into relations (25a), and considering terms only up to $O(x^3)$ we obtain the following matrix relation governing the leading coefficients of the potential function:

$$ \begin{bmatrix} K^+ & L_2^+ & L_3^+ & 0 & 0 & 0 & 0 & 0 \\ K^- & L_2^- & L_3^- & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & M_1^+ & M_2^+ & 0 & 0 & 0 \\ 0 & 0 & 0 & M_1^- & M_2^- & 0 & 0 & 0 \\ 0 & N_2^+ & N_3^+ & 0 & 0 & O_0^+ & O_1^+ & 0 \\ 0 & N_2^- & N_3^- & 0 & 0 & O_0^- & O_1^- & 0 \\ 0 & 0 & 0 & P_1^+ & P_2^+ & 0 & 0 & Q_0^+ \\ 0 & 0 & 0 & P_1^- & P_2^- & 0 & 0 & Q_0^- \end{bmatrix} \begin{Bmatrix} h \\ c_{02} \\ c_{03} \\ c_{11} \\ c_{12} \\ c_{20} \\ c_{21} \\ c_{30} \end{Bmatrix} = \begin{Bmatrix} D_1^+ \\ D_1^- \\ D_2^+ \\ D_2^- \\ D_3^+ \\ D_3^- \\ D_4^+ \\ D_4^- \end{Bmatrix}, \quad (27) $$

where the various coefficients are defined in Appendix A. Assuming that all coefficients $c_{ij}$, $i \ge 4$, are set equal to zero, the above matrix relation governs exactly the inverse problem and involves no approximation. The vector on the right-hand-side depends on an infinite number coefficients $c_{ij}$ (cf. Appendix A), and, hence, the inverse problem possesses non-unique solutions. Furthermore, it can be realized that problem (27) can be decomposed into two separate subproblems, one governing the coefficients $c_{1j}$ and $c_{3p}$ with $j \ge 1$, $p \ge 0$, and the other governing the coefficients $c_{0k}$ and $c_{2p}$ with $k \ge 2$, $p \ge 0$. The solutions of these subproblems are summarized below.

(1) Subproblem I governing $c_{1j}$ and $c_{3p}$ with $j \ge 1$, $p \ge 0$.

* If $c_{1j} = 0$, $j \ge 1$, then, either, $c_{3p} = 0$, $p \ge 0$

or, $c_{30} = D_4^+/Q_0^+$, and $c_{3j}$, $j \ge 1$ satisfy $\sum_{j \ge 1} c_{3j}(Q_j^+ + Q_j^-) = 0$.
• If not all $c_{1j}, j \ge 3$ are zero, then the coefficients are computed by the relations,

$$ \begin{cases} c_{11} \\ c_{12} \end{cases} = \begin{bmatrix} M_1^+ & M_2^+ \\ M_1^- & M_2^- \end{bmatrix}^{-1} \begin{Bmatrix} D_2^+ \\ D_2^- \end{Bmatrix}, \quad c_{30} = (D_4^+ - P_1^+ c_{11} - P_2^+ c_{12})/Q_0^+ $$

and $c_{3j}, j \ge 1$ satisfy the relation

$$ (D_4^+ - P_1^+ c_{11} - P_2^+ c_{12})/Q_0^+ = (D_4^- - P_1^- c_{11} - P_2^- c_{12})/Q_0^- $$

(2) Subproblem II governing $c_{0k}$ and $c_{2p}$ with $k \ge 2, p \ge 0$.

In this case we cannot have the trivial solution $c_{0j} = 0, j \ge 2$ and these coefficients satisfy

$$ \begin{cases} c_{02} \\ c_{03} \end{cases} = \begin{bmatrix} L_2^+ & L_3^+ \\ L_2^- & L_3^- \end{bmatrix}^{-1} \begin{Bmatrix} -hK^+ + D_1^+ \\ -hK^- + D_1^- \end{Bmatrix}. $$

the remaining coefficients satisfy the relation

$$ \begin{cases} c_{20} \\ c_{21} \end{cases} = \begin{bmatrix} O_0^+ & O_1^+ \\ O_0^- & O_1^- \end{bmatrix}^{-1} \left( \begin{bmatrix} N_2^+ & N_3^+ \\ N_2^- & N_3^- \end{bmatrix} \begin{Bmatrix} -c_{02} \\ -c_{03} \end{Bmatrix} + \begin{Bmatrix} D_3^+ \\ D_3^- \end{Bmatrix} \right). $$

All matrices in the solutions above are invertible, provided that the parameter $a$ of the trajectory (26) is nonzero. This eliminates the possibility of a trivial trajectory in the configuration plane.

A solution to the inverse problem can be formed by combining one of the two solutions of subproblem 1 with the solution of subproblem 2. We note that since the above solutions involve an infinite number of coefficients, one can obtain an infinity of solutions to the inverse problem by assigning arbitrary values to all the coefficients except the ones required to satisfy the above relations. Perhaps the simplest of these solutions is obtained by setting $c_{1j} = 0, j \ge 1$ and $c_{3p} = 0, p \ge 0$ (solution of subproblem 1), and $c_{02} = 2h/(2a^2 + e^2)$, $c_{0j} = 0, j \ge 3$, $c_{20} = h/2(2a^2 + e^2)$, $c_{2j} = 0, j \ge 1$ (solution of subproblem 2). Hence, we obtain a linear system with clearance nonlinearities, which is predicted to possess the vibro-impact periodic solution (26). In Figure 12 we superimpose the theoretical orbit (26) and the numerical orbit obtained by direct numerical integration of the vibro-impact equations of

<figure><img src="image_15.png" /><figcaption>Figure 12. Comparison between theoretical and numerical solutions of the inverse problem.</figcaption></figure>
motion. The agreement between theory and numerical integration validates the previous inverse analysis.

Summarizing, in this section we formulated an inverse vibro-impact problem by computing the required potential function of a system in order to produce a specified trajectory in the configuration plane. Although the analysis herein was performed for a double-branch orbit of parabolic shape, the technique can be applied to a more general class of orbits with single- or double-sided impacts. An interesting aspect of the inverse problem analyzed is the degeneracy of its solutions: indeed, it was found that an infinity of vibro-impact systems can produce the specified trajectory, a result which is due, perhaps, to its relatively simple parabolic shape. It is conjectured that the degeneracy of solutions of the inverse problem decreases as the shape of the trajectory becomes more complicated. This is an interesting topic for some future research.

## 4. DISCUSSION

We studied direct and inverse problems associated with vibro-impact oscillations of discrete systems. The usefulness of the presented analytical/numerical techniques stems from the fact that the majority of previous vibro-impact works were direct problems and centered on dynamical systems with a single impacting coordinate. By contrast, in this work we examined systems with two impacting coordinates, undergoing single- or double-sided impacts. The analytical methods employed were a NNM-type analysis and a BVP formulation. In the NNM-type analysis the displacement of one of the particles was expressed as a function of the displacement of the other particle, and the functional relation of the motions of the two particles was asymptotically approximated. The BVP enabled the computation of various branches of bifurcating time-periodic solutions with different impacting properties. Numerical results were obtained by direct integrations of the equations of motion and then used to construct Poincaré maps.

The class of vibro-impact systems considered possesses rich nonlinear dynamical behavior, including vibro-impact localized and nonlocalized time-periodic motions (which can be considered as analytic continuations in the vibro-impact regime of NNMs of corresponding smooth dynamical systems), as well as, complicated bifurcation sequences giving rise to new types of single- or double-sided impacting motions, instabilities, and chaotic responses. These bifurcations can be studied using appropriate BVPs and numerical Poincaré maps.

We also formulated inverse vibro-impact problems, whereby we sought systems that would produce specified orbits in the configuration plane. As expected from experience with similar types of inverse problems for other engineering areas, the solution of the vibro-impact inverse problem is generally non-unique, since it can be reduced to an underdetermined set of algebraic equations with multiple infinities of unknowns. The inverse problems discussed in this work can find applications to tracking problems were the design task calls for the responses of the coordinates of a system to track specified paths. Such problems arise often in kinematics and dynamics of mechanisms, and in other engineering applications.

The presented techniques can be extended to vibro-impact oscillators with more than two degrees-of-freedom, although in such cases the computational and numerical effort is expected to increase. Moreover, the outlined analysis can be used to study bifurcations associated with the inverse problem, i.e., to investigate degeneracies in the class of dynamical vibro-impact systems that produce a specified orbit in configuration space. That would be an interesting problem in inverse nonlinear dynamics. In addition, the mentioned analytical and numerical techniques can be applied to study the bifurcation structure of
the periodic solutions of multi-degree-of-freedom vibro-impact systems, and to establish
regions of regular and chaotic responses.

ACKNOWLEDGMENTS

This work was supported in part by NSF Young Investigator Award CMS-94-57750 (A.F.V.). Dr Devendra Garg is the Grant monitor.

REFERENCES

1. J. M. T. THOMPSON and R. GHAFFARI 1982 *Physics Letters* **91A**, 5-8. Chaos after period-doubling bifurcations in the resonance of an impact oscillator.

2. F. MOON and S. W. SHAW 1983 *International Journal of Non-linear Mechanics* **18**(6), 465-477. Chaotic vibrations of a beam with nonlinear boundary conditions.

3. S. W. SHAW and R. H. RAND 1989 *International Journal of Non-linear Mechanics* **24**(1), 41-56. The transition to chaos in a simple mechanical system.

4. S. F. MASRI and T. K. CAUGHEY 1966 *Journal of Applied Mechanics* **33**, 586-592. On the stability of the impact damper.

5. S. W. SHAW and P. HOLMES 1989 *Physical Review Letters* **51**, 623-626. A periodically forced linear oscillator with impacts: chaos and long period motions.

6. S. W. SHAW 1985 *Journal of Applied Mechanics* **52**, 453-464. The dynamics of a harmonically excited system having rigid amplitude constraints, Part I: subharmonic motions and local bifurcations, and Part II: chaotic motions and global bifurcations.

7. J. SHAW and S. W. SHAW 1989 *Journal of Applied Mechanics* **56**, 168-174. The onset of chaos in a two-degree-of-freedom impacting system.

8. A. P. IVANOV 1993 *Prikl. Mat. Mekh.* **57**(2), 5-21. Analytical methods in the theory of vibro-impact systems.

9. V. F. ZHURAVLEV 1977 *Iz. AN SSR, Mekhanica Tverdogo Tela* **12**(6), 24-28. Investigation of certain vibration-impact systems by the method of nonsmooth transformations.

10. E. G. VEDENOVA, L. I. MANEVITCH and V. N. PILIPCHUK 1985 *Prikl. Mat. Mekh.* **49**(2), 203-211. Normal oscillations of a string with concentrated masses on nonlinearly elastic supports.

11. A. F. VAKAKIS, L. I. MANEVITCH, Y. V. MIKHLIN, V. N. PILIPCHUK and A. A. ZEVIN 1996 *Normal Modes and Localization in Continuous Systems*. New York: Wiley Interscience.

12. E. L. B. VAN DER VORST, D. H. VAN CAMPEN and A. DE KRAKER 1996 *Journal of Sound and Vibration* **192**(5), 913-925. Periodic solutions of a multi-dof beam system with impact.

13. E. EMACI, T. A. NAYFEH and A. F. VAKAKIS 1997 *Zeitschrift fur Angewandte Mathematik und Mechanik*—ZAMM **77**(7), 527-541. Numerical and experimental study of nonlinear localization in a flexible structure with vibro-impacts.

14. M. A. F. AZEEZ and A. F. VAKAKIS 1997 16th ASME Conference on *Mechanical Vibration and Noise, Symposium on Impact*, Sacramento, CA, 14-17 September. Numerical and experimental analysis of the nonlinear dynamics due to impacts of a continuous overhung rotor.

15. C. BUDD, F. DUX and A. CLIFFE 1995 *Journal of Sound and Vibration* **184**(3), 475-502. The effect of frequency and clearance variations on single-degree-of-freedom impact oscillators.

16. F. MOON, W. HOLMES and S. KHOURY 1991 *Chaos* **1**(1), 65-68. Symbol dynamic maps of spatial-temporal chaotic vibrations in a string of impact oscillators.

17. A. P. IVANOV 1994 *Journal of Sound and Vibration* **178**(3), 361-378. Impact oscillations: linear theory of stability and bifurcations.

18. F. PETERKA and J. VACIK 1992 *Journal of Sound and Vibration* **154**(1), 95-115. Transition to chaotic motion in mechanical systems with impacts.

19. T. P. VALKERING 1994 *Journal of Sound and Vibration* **175**(3), 397-422. Non-trivial dynamics in a driven string with impact nonlinearity.

20. F. FARAHANCHI and S. H. SHAW 1994 *Journal of Sound and Vibration* **177**(3), 307-324. Chaotic and periodic dynamics of a slider-crank mechanism with slider clearance.

21. R. P. S. HAN and A. C. J. LUO 1995 *Journal of Sound and Vibration* **181**(2), 231-250. Chaotic motion of a horizontal impact pair.

22. G. S. WHISTON 1992 *Journal of Sound and Vibration* **152**(3), 427-460. Singularities in vibro-impact dynamics.
23. B. BROGLIATO 1996 *Nonsmooth Impact Mechanics, Lecture Notes in Control and Information Sciences* 220. Berlin: Springer.

24. J. GUCKENHEIMER and P. HOLMES 1984 *Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*. Berlin: Springer.

25. S. WIGGINS 1990 *Introduction to Applied Nonlinear dynamical Systems and Chaos*. Berlin: Springer.

26. A. F. VAKAKIS and R. H. RAND 1992 *International Journal of Nonlinear Mechanics* **27**(5), 861–888. Normal modes and global dynamics of a two-degree-of-freedom nonlinear system I: Low energies, and II: High energies.

27. A. F. VAKAKIS 1992 *Journal of Sound and Vibration* **158**(2), 341–361. Nonsimilar normal oscillations in a strongly nonlinear discrete system.

## APPENDIX A: COEFFICIENTS OF THE MATRIX EQUATION (23)

The coefficients are expressed in terms of the geometric parameters of the trajectory (22) as follows:

$$K^{\pm} = (\mp 4a/e^2),$$

$$L_j^{\pm} = -4(\mp a/e^2)(\pm a)^j + j(\pm a)^{j-1}, \quad j = 2, 3, \dots,$$

$$M_j^{\pm} = -6(\mp a/e^2)(\pm a)^j + j(\pm a)^{j-1}, \quad j = 1, 2, \dots,$$

$$N_j^{\pm} = j(j-1)(\pm a)^{j-2}(\mp a/e^2), \quad j = 2, 3, \dots,$$

$$O_j^{\pm} = -8(\mp a/e^2)(\pm a)^j + j(\pm a)^{j-1}, \quad j = 0, 1, \dots,$$

$$P_j^{\pm} = -2j(a^2/e^4)(\pm a)^{j-1} - 8(\mp a/e^2)(a^2/e^4)(\pm a)^j + j(j-1)(\pm a)^{j-2}(\mp a/e^2)$$

$$j = 1, 2, \dots,$$

$$Q_j^{\pm} = -10(\mp a/e^2)(\pm a)^j + j(\pm a)^{j-1} \quad j = 0, 1, \dots,$$

$$D_1^{\pm} = -\sum_{j \ge 4} c_{0j} L_j^{\pm},$$

$$D_2^{\pm} = -\sum_{j \ge 3} c_{1j} M_j^{\pm},$$

$$D_3^{\pm} = -\sum_{j \ge 4} c_{0j} N_j^{\pm} - \sum_{j \ge 2} c_{2j} O_j^{\pm},$$

$$D_4^{\pm} = -\sum_{j \ge 3} c_{1j} P_j^{\pm} - \sum_{j \ge 1} c_{3j} Q_j^{\pm}.$$