# Non-periodic motion caused by grazing incidence in an impact oscillator

Arne Nordmark

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

► To cite this version:

Arne Nordmark. Non-periodic motion caused by grazing incidence in an impact oscillator. Journal of Sound and Vibration, 1991, 145 (2), pp.279-297. 10.1016/0022-460X(91)90592-8. hal-02546986

HAL Id: hal-02546986

https://hal.science/hal-02546986

Submitted on 19 Apr 2020

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>
# NON-PERIODIC MOTION CAUSED BY GRAZING INCIDENCE IN AN IMPACT OSCILLATOR

A. B. NORDMARK

Department of Mechanics, Royal Institute of Technology, S-100 44 Stockholm, Sweden

The motion of a single-degree-of-freedom, periodically forced oscillator subjected to a rigid amplitude constraint is considered. Using analytical methods, the singularities caused by grazing impact are studied. It is shown that as a stable periodic orbit comes to grazing impact under the control of a single parameter, a special type of bifurcation occurs. The motion after the bifurcation may be non-periodic, and a criterion for this based on orientation and eigenvalues is given.

## 1. INTRODUCTION

In recent years, single-degree-of-freedom impacting linear oscillators subjected to periodic external forcing have been intensively studied by using the methods of dynamical systems theory. Shaw and co-workers [1-4] have studied clearance systems and have analyzed subharmonic periodic motions and the bifurcations leading to them. In numerical simulations, their systems have also been found to exhibit symptoms of chaotic motion. A similar model with preload or clearance has been studied by Whiston [5, 6], who has provided an analysis of global bifurcations. There are also several experimental studies on impacting oscillators, for example, that of Moon [7].

As noted by Shaw and Holmes [1], the Poincaré mappings derived for impacting systems have singularities in the derivative. These singularities can clearly be seen in the one-dimensional mappings derived for a system with large dissipation [2]. Whiston has analyzed the singularities and shown their importance for the global dynamics of impact oscillators. The present work is focused on the effect of the singularities on a stable periodic orbit, showing by analytical methods the possibility of chaotic motion for these systems.

The oscillator model used in this paper is a general model for a single-degree-of-freedom, periodically forced impacting oscillator, where the amplitude is limited by a stop on one side. The impact process is idealized to take place in zero time, and with the velocity coming out of the stop being a function of incoming velocity only. The systems that can be described by this model include the linear impact oscillators and the problem of a ball bouncing on a vibrating table.

The first section of the paper consists of a presentation of the system studied, and a discussion of the geometry of the state space of the system and ways to introduce Poincaré sections and mappings. Locally valid, basic mappings between different sections of the state space are defined. These can then be composed to give locally defined Poincaré mappings near some orbit of the system.

In the Poincaré sections, there are curves of points where the mappings show a singularity in the Jacobian as the curves are approached from one side, while the Jacobians on the other side stays bounded and slowly varying. These curves represent points that
are brought to a zero velocity “impact” with the stop, here referred to as a *grazing impact*.
If a stable periodic point under the control of a single system parameter crosses one of
these curves from the non-singular side, the stability characteristics of the point change
suddenly. Because of the discontinuous and singular change in derivatives as the curve
is crossed, this bifurcation is not one of the well known bifurcations that frequently occur
in the study of smooth non-linear systems [8]. The later part of the paper consists of a
study of this bifurcation. At the bifurcation point two cases are distinguished: either the
periodic point suddenly loses its stability, causing the motion of the system to leave the
immediate neighborhood of the bifurcation point; or the point will retain its stability. In
the latter case, a trapping region is found to exist after the bifurcation. The size of the
trapping region grows continuously from the bifurcation point, and a simple analysis of
stability gives an inequality condition on the largest eigenvalue of the periodic point,
under which no stable periodic points are possible in the trapping region for an interval
of the bifurcation parameter. For other regions of the eigenvalues, the analysis shows the
possibility of a sequence of points of arbitrarily long periods, and the corresponding
sequence of bifurcation parameter values is found to obey a geometric scaling law.

## 2. THE SYSTEM

The system treated is a single-degree-of-freedom impact oscillator that is being driven
by a periodic force. As it is a single-degree-of-freedom system, being externally forced,
the acceleration of the mass is a function of displacement, velocity and time. The
periodicity of the forcing requires the acceleration function to be periodic in the time
variable. At a fixed value of the displacement co-ordinate, the mass strikes a stop, causing
it to rebound. The impact process is considered to be so rapid that the time spent in the
rebound can be taken to be zero. The velocity coming out of the stop is assumed to
depend only on the velocity coming in. The time co-ordinate is denoted by $\tau$ and the
displacement by $x$. With scaling and shifting the co-ordinates if necessary, it can be
assumed that the time period of the forcing is $2\pi$, the impact boundary is given by $x = \bar{x}_c$,
and the oscillator moves on the side $x \ge \bar{x}_c$.

### 2.1. FORMULATION OF DIFFERENTIAL EQUATION

The equation governing the motion is a second order differential equation:

$$ \frac{d^2 x}{d\tau^2} = A(x, \frac{dx}{d\tau}, \tau), \quad (1) $$

and the periodicity of the acceleration function is expressed by

$$ A(x, \frac{dx}{d\tau}, \tau + 2\pi) = A(x, \frac{dx}{d\tau}, \tau). $$

### 2.2. THE SYSTEM AS A THREE-DIMENSIONAL FLOW

Equation (1) describes a one-dimensional non-autonomous second order system, but
the system can equivalently be described as a three-dimensional autonomous first order
system in the following standard manner. Define a new time variable $t \stackrel{\text{def}}{=} \tau - \tau_0$, where $\tau_0$
can be freely chosen. Then define the three state variables: namely $x$ (still displacement),
$v \stackrel{\text{def}}{=} \frac{dx}{d\tau}$, and $\theta \stackrel{\text{def}}{=} \tau \pmod{2\pi}$. The flow equations for the new variables are

$$ \frac{d}{dt} \begin{pmatrix} x \\ v \\ \theta \end{pmatrix} = \begin{pmatrix} v \\ A(x, v, \theta) \\ 1 \end{pmatrix}. \quad (2) $$

This gives a picture of the system without reference to an absolute time, so all the dynamics
is in the geometry of the flow.
It is worth discussing the geometry of the state space and the flow in some detail. $x$ and $v$ are real numbers, and $\theta$ is a circular phase, so the state space is a subset of $R \times R \times S^1$. It can be visualized in $R^3$ by using cylindrical co-ordinates with the radius taken as $x-x_c+r_0$, $r_0$ a positive constant, the angle as $\theta$, and the height as $v$ (see Figure 1).

It is convenient to introduce notations for some subsets of the state space:

$$ 
\begin{aligned}
E^{\text{def}}(\bar{x}_c, \infty) \times R \times S^1, & \qquad \Theta^{\bar{\theta}}^{\text{def}}[\bar{x}_c, \infty) \times R \times \{\bar{\theta}\}, & \qquad \Sigma^{\text{def}}\{\bar{x}_c\} \times R \times S^1, \\
\Sigma^{+\text{def}}\{\bar{x}_c\} \times (0, \infty) \times S^1, & \qquad \Sigma^{0\text{def}}\{\bar{x}_c\} \times \{0\} \times S^1, & \qquad \Sigma^{-\text{def}}\{\bar{x}_c\} \times (-\infty, 0) \times S^1.
\end{aligned}
$$

The flow (2) is defined in $E$, $\Theta^{\bar{\theta}}$ is a plane of some constant phase $\bar{\theta}$, and $\Sigma$ is the impact boundary. On $\Sigma^-$ the oscillator hits the boundary with negative velocity, and on $\Sigma^+$ it leaves the boundary with positive velocity. $\Sigma^0$ represents the state of being at the boundary with zero velocity, and the possibility for a flow line to enter or leave a point in $\Sigma^0$ is determined by the value of the acceleration function $A$ (in $\Sigma^0$ a function of phase only).

If $A$ is positive (directed away from the impact boundary) for a phase interval around a point in $\Sigma^0$, then the flow can enter and leave that point. This is the situation of a *grazing impact*; the oscillator just touches the boundary, without being effected by it. A periodic orbit, projected onto the $x-v$ plane, with a grazing impact, is shown in Figure 2(a). As will be shown later, grazing impact has a large influence on the stability of the motion of the system.

<figure><img src="image_3.png" /><figcaption>Figure 1. The state space.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Periodic orbit with (a) a grazing impact, and (b) an infinite number of impacts.</figcaption></figure>
If A is negative in $\Sigma^0$, the flow cannot leave or enter that point. This happens when the oscillator is at rest at the boundary, being pressed against it by the force. The motion will thus stay in $\Sigma^0$ until the acceleration becomes positive again. Although such points are not at the end of any flow line in E, they could be reached by the motion. For example, if there is damping in the impact, they might be reached through an infinite number of impacts in finite time. This is illustrated in Figure 2(b). However, a motion coming through one of these points cannot in general be uniquely followed backwards in time.

The formal solution of the flow equations (2) for initial conditions $x_0$, $v_0$, $\theta_0$ given at $t=0$ can be written in terms of a displacement function, x, as

$$ \begin{pmatrix} x \\ v \\ \theta \end{pmatrix} = \begin{pmatrix} x(t, x_0, v_0, \theta_0) \\ \frac{\partial x}{\partial t}(t, x_0, v_0, \theta_0) \\ t + \theta_0 \pmod{2\pi} \end{pmatrix}, \qquad (3) $$

where

$$ x(0, x_0, v_0, \theta_0) = x_0, \quad \frac{\partial x}{\partial t}(0, x_0, v_0, \theta_0) = v_0, \quad \text{and} \quad \frac{\partial^2 x}{\partial t^2} = A \begin{pmatrix} x, \frac{\partial x}{\partial t}, t + \theta_0 \end{pmatrix}. $$

The impact law is described by a function $G, G:(-\infty, 0]\rightarrow[0, \infty)$, with the properties $v_{out}=G(v_{in})$, and $G(0)=0$. The impact law maps points in $\Sigma^-$ to points in $\Sigma^+$ with the displacement and phase unchanged.

## 2.3. POINCARÉ MAPPINGS

As an aid in describing and understanding the system one can introduce mappings between different sections of the state space (Poincaré mappings). There are two different kinds of sections that are natural for this system. One is the constant phase section $\Theta^{\bar{\theta}}$, where a point is determined by its displacement $x_{\bar{\theta}}$ and velocity $v_{\bar{\theta}}$ in the plane corresponding to a given, constant phase $\bar{\theta}$. This kind of section can be used in any problem with periodically forced oscillators. The other choice is connected to the impact process, and the section used is $\Sigma$, where a point is determined by its phase $\theta_c$ and velocity $v_c$. It should be noted that when using these sections as a base for defining a globally valid mapping, care has to be taken to make the definition precise, and in general the mapping will be discontinuous. The cause of discontinuity is, however, different for the two choices of section.

In the case of the constant phase plane mapping, the discontinuity is due to the fact that a point with $x = \bar{x}_c$ is really a pair of points, one with negative and one with positive velocity, since the impact process is considered to be instantaneous. The mapping will thus be discontinuous at points which are mapped to $x = \bar{x}_c$. However, the control over which constant phase plane is to be used as the section might be used to move the discontinuities out of a regions of interest. This discontinuity is then just an artifact of having an instantaneous impact process, and can be regarded as a limit of a steep, continuous mapping when the impact boundary stiffness goes to infinity. In the present work, treating only local regions, it will not be a cause of trouble.

In the case of mapping between impacts, there is a more fundamental problem due to tangential flow. At $v=0$, the flow is tangent to the Poincaré section, so that nearby points may impact with low velocity, or may miss and impact at a later time. This is an inherent problem when choosing a section where the flow is not everywhere transversal, and it will cause discontinuities to appear at points being mapped to the line $v=0$.

The local Poincaré mappings derived later in the paper will use a constant phase plane as the section.
## 2.4. LOCAL MAPPINGS

The careful definition of globally valid mappings will not be required for the present purpose. Instead, four basic, locally defined classes of mappings between different sections in state space will be studied. These mappings can then be composed to give a local Poincaré mapping in the neighborhood of a trajectory in state space. The primary goal is to obtain expressions for the derivative of the mappings, since this information is used for stability calculations, both in analytical work and in numerical simulations. The local decomposition of a Poincaré mapping into a sequence of mappings from four basic classes is convenient for calculation, and distinguishes between the contributions from the flow and those from the impact process.

The derivative of the final point with respect to the starting point when integrating flow equations such as (2) for a fixed interval of time is obtained by standard methods (the variational equations). The time taken to reach the impact boundary is, however, dependent on the starting point. Thus there is a need for an expression for the derivative when going from a plane of constant time to some other type of surface in state space. The following theorem, illustrated by Figure 3, can be used to determine the derivative of the projection mapping implied by the flow from one surface to another, at a point which lies on both surfaces.

<figure><img src="image_5.png" /><figcaption>Figure 3. Flow projection onto a surface.</figcaption></figure>

Suppose a flow in $R^n$ is given by $dx/dt = F(x)$ and a smooth surface $S$ in $R^n$ is given by $H(x) = 0$, and let $\bar{x}$ be a point in $S$. The points in $R^n$ and the functions $F$ are represented by column vectors and gradients by row vectors. Then the following holds. If

$$ (\partial H / \partial x) F \neq 0 \qquad (4) $$

at $\bar{x}$ (the flow is not tangent to $S$ at $\bar{x}$), then there is a neighborhood $B$ of $\bar{x}$ and a mapping $P: B \to S$, such that for $x \in B$, $x$ and $P(x)$ is connected by a part of a flow line, contained in $B$. The mapping $P$ is differentiable at $\bar{x}$, and the derivative is

$$ \partial P / \partial x = I - F(\partial H / \partial x) / (\partial H / \partial x) F, \qquad (5) $$

where $I$ is the $n \times n$ identity matrix.

To establish this result, one notes that the time of flight from $x$ to $P(x)$ is implicitly defined by the condition that $P(x) \in S$. A direct calculation with use of the implicit function theorem will then give the expression (5) for the derivative. It can be noted that the linear part (5) describes linear projection onto a plane tangent to $S$ at $\bar{x}$ in the direction of the flow at $\bar{x}$. When using this result to obtain a mapping from one surface ($T$) through $\bar{x}$ to another ($S$), local co-ordinates in $T$ are transformed to $R^n$, subjected to the mapping $P$ and then transformed to local co-ordinates in $S$.

The four classes of local mappings will now be defined. The position of the sections in state space near a typical trajectory is shown in Figure 4.
<figure><img src="image_6.png" /><figcaption>Figure 4. The local mappings.</figcaption></figure>

I. A flow mapping $P_I$ from a constant phase plane $\Theta^{\tilde{\theta}_1}$ to a constant phase plane $\Theta^{\tilde{\theta}_2}$, defined near a trajectory starting at $(\tilde{x}_{\tilde{\theta}_1}, \tilde{v}_{\tilde{\theta}_1}, \tilde{\theta}_1)$ and reaching $(\tilde{x}_{\tilde{\theta}_2}, \tilde{v}_{\tilde{\theta}_2}, \tilde{\theta}_2)$ with a time of flight $\tilde{t}$ and no impacts:

$$P_I : \Theta^{\tilde{\theta}_1} \to \Theta^{\tilde{\theta}_2}, \quad \begin{pmatrix} x_{\tilde{\theta}_1} \\ v_{\tilde{\theta}_1} \end{pmatrix} \mapsto \begin{pmatrix} x_{\tilde{\theta}_2} \\ v_{\tilde{\theta}_2} \end{pmatrix}, \quad \text{with} \quad P_I \begin{pmatrix} \tilde{x}_{\tilde{\theta}_1} \\ \tilde{v}_{\tilde{\theta}_1} \end{pmatrix} = \begin{pmatrix} \tilde{x}_{\tilde{\theta}_2} \\ \tilde{v}_{\tilde{\theta}_2} \end{pmatrix}.$$

The formal solution of the flow equations is given by equation (3), so $\tilde{x}_{\tilde{\theta}_2} = x(\tilde{t}, \tilde{x}_{\tilde{\theta}_1}, \tilde{v}_{\tilde{\theta}_1}, \tilde{\theta}_1)$, $\tilde{v}_{\tilde{\theta}_2} = (\partial x / \partial t)(\tilde{t}, \tilde{x}_{\tilde{\theta}_1}, \tilde{v}_{\tilde{\theta}_1}, \tilde{\theta}_1)$, and $\tilde{\theta}_2 = \tilde{t} + \tilde{\theta}_1 (\mod 2\pi)$. The Jacobian of $P_I$ at $(\tilde{x}_{\tilde{\theta}_1}, \tilde{v}_{\tilde{\theta}_1})$ is

$$JP_I = \begin{pmatrix} \partial x / \partial x_0 & \partial x / \partial v_0 \\ \partial^2 x / \partial t \partial x_0 & \partial^2 x / \partial t \partial v_0 \end{pmatrix}.$$

If the function $x$ is not known explicitly, the Jacobian can be found by integrating the variational equations numerically during the time $\tilde{t}$.

II. A flow projection mapping $P_{II}$ from a constant phase plane $\Theta^{\tilde{\theta}}$ to $\Sigma^{-}$, defined near a point $\tilde{p} = (\tilde{x}_c, \tilde{v}, \tilde{\theta})$ with $\tilde{v} < 0$. In the notation used in the above theorem, the vector field $F$ is given by (2), the surface $S$ is $\Sigma^{-}$, and the function $H$ is given as $x - \tilde{x}_c = 0$. The non-tangential condition (4) becomes $\tilde{v} \neq 0$, which is satisfied. Thus there is a local mapping from points near $\tilde{p}$ to $\Sigma^{-}$, and if restricted to points in $\Theta^{\tilde{\theta}}$, the mapping becomes:

$$P_{II} : \Theta^{\tilde{\theta}} \to \Sigma^{-}, \quad \begin{pmatrix} x_{\tilde{\theta}} \\ v_{\tilde{\theta}} \end{pmatrix} \mapsto \begin{pmatrix} \theta_c \\ v_c \end{pmatrix}, \quad \text{with} \quad P_{II} \begin{pmatrix} \tilde{x}_c \\ \tilde{v} \end{pmatrix} = \begin{pmatrix} \tilde{\theta} \\ \tilde{v} \end{pmatrix}.$$

The derivative according to equation (5) is

$$\begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} - \frac{1}{v} \begin{pmatrix} \tilde{v} & 0 & 0 \\ \tilde{A} & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix},$$

where $\tilde{A} \stackrel{\text{def}}{=} A(\tilde{x}_c, \tilde{v}, \tilde{\theta})$, and when this is restricted to the rows corresponding to $\theta$ and $v$, and the columns corresponding to $x$ and $v$, the Jacobian of the mapping $P_{II}$ at $(\tilde{x}_c, \tilde{v})$ becomes

$$JP_{II} = -\frac{1}{\tilde{v}} \begin{pmatrix} 1 & 0 \\ \tilde{A} & -\tilde{v} \end{pmatrix}.$$
### III. The impact mapping $P_{III}$ from $\Sigma^-$ to $\Sigma^+$, globally defined:

$$P_{III}: \Sigma^- \to \Sigma^+, \quad \begin{pmatrix} \theta_{c-} \\ v_{c-} \end{pmatrix} \mapsto \begin{pmatrix} \theta_{c+} \\ v_{c+} \end{pmatrix},$$

where $\theta_{c+} = \theta_{c-}$ and $v_{c+} = G(v_{c-})$. The Jacobian is

$$JP_{III} = \begin{pmatrix} 1 & 0 \\ 0 & dG/dv_{in} \end{pmatrix}.$$

### IV. A flow projection mapping from $\Sigma^+$ to a constant phase plane $\Theta^{\bar{\theta}}$, defined near a point $\bar{p} = (\bar{x}_c, \bar{v}, \bar{\theta})$ with $\bar{v} > 0$. The surface $S$ is $\mathcal{S}^{\bar{\theta}}$, given by $\theta - \bar{\theta} = 0$. The non-tangential condition (4) gives $1 \neq 0$, and restricted to $\Sigma^+$, the mapping becomes:

$$P_{IV}: \Sigma^+ \to \Theta^{\bar{\theta}}, \quad \begin{pmatrix} \theta_c \\ v_c \end{pmatrix} \mapsto \begin{pmatrix} x_{\bar{\theta}} \\ v_{\bar{\theta}} \end{pmatrix}, \quad \text{with} \quad P_{IV} \begin{pmatrix} \bar{\theta} \\ \bar{v} \end{pmatrix} = \begin{pmatrix} \bar{x}_c \\ \bar{v} \end{pmatrix}.$$

The Jacobian at $(\bar{\theta}, \bar{v})$ is found to be

$$JP_{IV} = \begin{pmatrix} -\bar{v} & 0 \\ -\bar{A} & 1 \end{pmatrix}.$$

A Poincaré mapping, either of the constant phase type or the impact type, can locally be constructed out of compositions of mappings from the four classes above. Mappings of class I are used to go from a constant phase plane through the point where the impact surface is left to a constant phase plane through the point of the next impact. A class II then goes from the constant phase plane to $\Sigma^-$, a class III to $\Sigma^+$, and a class IV to a constant phase plane through the point where the impact boundary is left, where the motion can again be continued with a class I mapping, and so on. A Poincaré mapping of the constant phase type thus locally consists of a series of mappings of class I, alternating with mappings of type $P_{IV} \circ P_{III} \circ P_{II}$. A Poincaré mapping of the impact type, mapping $\Sigma^-$ to $\Sigma^-$, would look like $P_{II} \circ P_I \circ P_{IV} \circ P_{III}$.

In this way these local mappings can be used to compute the stability of an orbit. It should be noted that the determinant of the Jacobian of a mapping of type $P_{IV} \circ P_{III} \circ P_{II}$ is $(v_{out}/v_{in})(dG/dv_{in})$. This expression is not large, even if $v_{in}$ is small. On the other hand, the Jacobian has elements that are very large, due to the factor $1/v_{in}$ in $JP_{II}$. Although the area of a constant phase region is changed by a factor of order one by a low-velocity impact, it will undergo a large stretching and compression. The effect of the mapping $P_{II}$ when $v_{in}$ is small can be compared with the shadow cast by a plane object onto the ground when the sun is standing low over the horizon, causing expansion in one direction, while the effect of the mapping $P_{IV}$ when $v_{out}$ is small can be compared to the shadow cast by a plane object which is held almost parallel to the light rays, causing contraction in one direction. Thus the projection mappings $P_{II}$ and $P_{IV}$ clearly plays an important role for stability, and cannot be neglected. This is also true for models in which the impact takes place over finite time, with velocity continuous across the boundary. If the acceleration is discontinuous across the boundary, Jacobians of type II and IV must be used as well as those obtained from the variational equations to calculate stability.

### 3. LOCAL MAPPINGS NEAR GRAZING IMPACT

In this section, low-order series expansions of mappings around an orbit having a grazing impact will be derived. In the derivation, it is always assumed that the generic case is at hand, for example that terms in series expansions are non-zero. The acceleration
function A and the impact law function G are considered to be sufficiently many times continuously differentiable whenever necessary, and the same for the inverse of G (which is assumed to exist).

The situation considered is one where there is a trajectory of the system that has a grazing impact at a point $O_c: (\bar{x}_c, 0, \bar{\theta}_c)$ in $\Sigma^0$ (see Figure 5). Two constant phase planes $\Theta^{\theta_1}$ and $\Theta^{\theta_2}$ are now chosen such that the trajectory first intersects $\Theta^{\theta_1}$ at a point $O_1: (\bar{x}_{\bar{\theta}_1}, \bar{y}_{\bar{\theta}_1}, \bar{\theta}_1)$, then has the grazing impact, and then intersects $\Theta^{\theta_2}$ at a point $O_2: (\bar{x}_{\bar{\theta}_2}, \bar{y}_{\bar{\theta}_2}, \bar{\theta}_2)$. It is assumed that there is only one grazing impact on the trajectory between the first and second planes. The task is to find representations for the local mapping from a neighborhood of $O_1$ in $\Theta^{\theta_1}$ to $\Theta^{\theta_2}$.

### 3.1. GEOMETRY NEAR GRAZING IMPACT

In Figure 6, the mapping of a region around $O_1$ in $\Theta^{\bar{\theta}_1}$, (a), to its image in $\Theta^{\bar{\theta}_2}$, (b), is shown. Mappings between constant phase planes are smooth except at points that lead to grazing impact. Points that impact with very low velocity acquire Jacobians that may have some very large elements, while points that just miss the impact will not be affected at all. At points leading to grazing impact the Jacobian is undefined. The locus of these points will form a curve separating impacting and non-impacting points. Likewise, the locus of the points coming from grazing impact will form a curve in the second constant phase plane.

At $O_c$, since the impact is grazing, the acceleration is positive. Upon looking at the Jacobian $JP_{IV}$, it is seen that a tangent vector in the positive $\theta$ direction in $\Sigma$ is mapped to a non-zero vector in the positive $v$ direction in $\Theta^{\bar{\theta}_c}$. Thus the mapping $P_{IV}$ maps a

<figure><img src="image_7.png" /><figcaption>Figure 5. A trajectory with a grazing impact.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Geometry of the mapping near grazing impact.</figcaption></figure>
segment of $\Sigma^0$ near $O_c$ to a smooth curve in $\Theta^{\tilde{\theta}_c}$, which is then mapped to a smooth curve $L_2$ in $\Theta^{\tilde{\theta}_2}$ through $O_2$. Moving backwards in time there is a smooth curve $L_1$ in $\Theta^{\tilde{\theta}_1}$ through $O_1$. On these curves a sense of direction is introduced in such a way that in moving in the positive direction on $L_1$, the time to reach grazing impact will decrease, corresponding to decreasing value of the phase at grazing impact. On $L_2$, the direction is the one implied by mapping $L_1$ to $L_2$, which means that, in moving in the positive direction along $L_2$, points coming from earlier grazing impacts are encountered. These directions makes it possible to say "to the right of a line", meaning on the right-hand side when facing in the positive direction. Using the fact that mappings of type I and II are orientation preserving, it is found that points to the left of $L_1$ are the impacting ones, while points to the right will miss. Likewise it is found that the impacting side of $L_2$ is again the left one. The mappings on the different sides of $L_1$ are totally different in character: on the non-impacting side the Jacobian is well behaved and slowly varying all the way up to $L_1$, while on the impacting side the Jacobian becomes singular as it approaches $L_1$. This means that a surface element near $L_1$ on the impacting side will be very distorted by the mapping. Furthermore, a smooth curve $A_1B_1C_1$ crossing $L_1$ at $B_1$ will be mapped to a curve $A_2B_2C_2$ with a corner at $B_2$. In fact, any curve segment directed from $L_1$ into the impacting side at a non-zero angle will be mapped to a segment tangent to and in the positive direction of $L_2$, and it will also be very stretched out along $L_2$.

## 3.2. ALGEBRAIC REPRESENTATIONS

The lines $L_1$ and $L_2$ are important for the mappings, and hence the co-ordinate systems should be chosen to simplify the form of these lines. In $\Theta^{\tilde{\theta}_1}$, new local co-ordinates $(\chi_1, \psi_1)$ are introduced by a smooth, orientation preserving co-ordinate transformation in such a way that $(\chi_1, \psi_1) = (0, 0)$ corresponds to $O_1$, $\chi_1 = 0$ to $L_1$, and $\chi_1 < 0$ to the impacting side of $L_1$. A similar co-ordinate system $(\chi_2, \psi_2)$ is introduced in $\Theta^{\tilde{\theta}_2}$. In $\Sigma$, local co-ordinates near $O_c$ are introduced by defining a local phase as $\phi_c \stackrel{\text{def}}{=} \theta_c - \tilde{\theta}_c$. The trailing part of a power series is represented by the symbol $\mathcal{O}$. For example, $\mathcal{O}(x, y, z)^3$ indicates terms that are of total order three or higher in the variables x, y and z, such as xyz or $x^2y$.

First look at the mapping of type $P_{IV}$ near $O_c$. The mapping is given by

$$x_{\tilde{\theta}_c} = x(-\phi_c, \bar{x}_c, v_c, \phi_c + \tilde{\theta}_c), \quad v_{\tilde{\theta}_c} = \frac{\partial x}{\partial t}(-\phi_c, \bar{x}_c, v_c, \phi_c + \tilde{\theta}_c).$$

Expanding in powers of $\phi_c$ and $v_c$ gives

$$x_{\tilde{\theta}_c} = \bar{x}_c + \bar{A}\phi_c^2/2 - \phi_c v_c + \mathcal{O}(\phi_c, v_c)^3, \quad v_{\tilde{\theta}_c} = -\bar{A}\phi_c + v_c + \mathcal{O}(\phi_c, v_c)^2,$$

where $\bar{A}$ is the acceleration at $O_c$. From this it is seen that the line segment $\phi_c = 0, v_c \ge 0$ is mapped tangent to the image of the line $v_c = 0$, in the direction of decreasing $\phi_c$, and that it curves away from it. Now the mapping from $\Theta^{\tilde{\theta}_c}$ to $\Theta^{\tilde{\theta}_2}$ is smooth, so it preserves these characteristics. By using these facts and the definition of the local co-ordinate system in $\Theta^{\tilde{\theta}_2}$, the mapping from $\Sigma^+$ to $\Theta^{\tilde{\theta}_2}$ is found to have the form:

$$\chi_2 = -\gamma v_c^2 [1 + \mathcal{O}(\phi_c, v_c)], \quad \psi_2 = -\beta \phi_c + \alpha v_c + \mathcal{O}(\phi_c, v_c)^2,$$

where $\alpha, \beta, \gamma > 0$.

Going backwards in time from $O_c$, the mapping from $\Sigma^-$ to $\Theta^{\tilde{\theta}_1}$ is found to have the same form:

$$\chi_1 = -c v_c^2 [1 + \mathcal{O}(\phi_c, v_c)], \quad \psi_1 = -b \phi_c + a v_c + \mathcal{O}(\phi_c, v_c)^2,$$

where $a, b, c > 0$.
Finally, the mapping from $\Sigma^-$ to $\Sigma^+$ has the form:

$$ \phi_{c+} = \phi_{c-}, \quad v_{c+} = -rv_{c-} + \mathcal{O}(v_{c-}^2), $$

where $r \stackrel{\text{def}}{=} -G'(0)$.

Putting these three mappings together to obtain the mapping from $\Theta^{\vec{\delta}_1}$ to $\Theta^{\vec{\delta}_2}$ (on the impact side $\chi_1 < 0$) gives

$$ \chi_2 = -\frac{\gamma r^2}{c} (-\chi_1)[1 + \mathcal{O}(\sqrt{-\chi_1}, \psi_1)], \quad \psi_2 = \frac{\alpha br + a\beta}{b\sqrt{c}} \sqrt{-\chi_1} + \frac{\beta}{b} \psi_1 + \mathcal{O}(\sqrt{-\chi_1}, \psi_1)^2. $$

On the non-impact side $\chi_1 > 0$ the mapping to first order is given by

$$ \begin{pmatrix} \chi_2 \\ \psi_2 \end{pmatrix} = A \begin{pmatrix} \chi_1 \\ \psi_1 \end{pmatrix} + \mathcal{O}(\chi_1, \psi_1)^2, $$

where $A$ is an invertible matrix that satisfies certain relations to make the full mapping continuous along the line $\chi_1 = 0$.

This is the lowest order form of the mapping from $\Theta^{\vec{\delta}_1}$ to $\Theta^{\vec{\delta}_2}$ near $O_1$, in which the singularities in derivative on the impact side are clearly seen.

### 4. GRAZING IMPACT IN A PERIODIC ORBIT, AND TRAPPING CONDITIONS

Now consider the case in which a stable periodic orbit of the system, with no low-velocity impacts, is close to a grazing impact somewhere in the orbit. If a system parameter is changed, then the shape of the periodic orbit changes with it, and it might happen that the orbit is brought to a grazing impact by this, so the situation looks like Figure 2(a). To study this phenomenon a Poincaré mapping $P$ is set up, mapping some constant phase plane $\Theta^{\vec{\delta}}$ back to itself after a time of flight equal to the period of the periodic orbit. This will make the mapping have fixed points where the orbit crosses the plane. One such fixed point $O^*$ is studied. Near this point there is a line $L_0$ of points leading to grazing impact and a line $L_1 (=P(L_0))$ of points coming from grazing impact. When the system parameter is changed, $O^*$ crosses $L_0$ from the non-impact side. It is here assumed that as the system parameter is changed, the fixed point crosses with non-zero speed, so the distance is to first order proportional to the offset of the system parameter. The offset of the system parameter is described by a normalized bifurcation parameter $\mu$, which is negative when $O^*$ is on the non-impact side, zero at the crossing, and positive when $O^*$ is on the impact side.

Since the Jacobian of $P$ changes discontinuously at $L_0$, with singularities on the impact side, the fixed point will in general loose stability as it crosses, and the dynamics will consequently change. It will be shown that under certain conditions the motion will *still be trapped near the former fixed point*, and the dynamics in this situation will be analyzed. Under other conditions, the motion will be mapped away from the region near the former fixed point, and in such cases a global analysis is needed to determine the resulting motion. No such analysis is presented here.

It should be noted, that there is no information in the stability of the periodic orbit predicting that it is close to a grazing impact situation. This is in contrast to the ordinary types of bifurcations, where some eigenvalues of the mapping approach the unit circle in the complex plane as the bifurcation is approached. In the case of grazing impact, the *bifurcation comes totally unexpectedly*, from a local point of view. From the local information of the mapping on the non-impact side, there is no way of predicting where the line $L_0$ is located.
### 4.1. THE MAPPING WHEN $\mu=0$

At $\mu=0$, the periodic orbit has a grazing impact, which means that $O^*$ is both coming from and leading to a grazing impact. Thus $L_0$ and $L_1$ must intersect at $O^*$. In the general case the intersection will be transversal, and since the lines are oriented there are two possible orientations for the crossing of the lines:

$L_0$ and $L_1$ form a right-hand system (the same orientation as the x- and v-axes), so the positive direction of $L_1$ points into the left-hand side of $L_0$. This will be called positive orientation of $L_0$ and $L_1$. In this case there are two local invariant half-manifolds in the impact side of $L_0$, both beginning at $O^*$. The unstable one is tangent to the positive direction of $L_1$ with the derivative along the manifold going to $+\infty$ as $O^*$ is approached, while the stable one is tangent to the negative direction of $L_0$ with derivative going to $0^+$. On the impact side, $O^*$ acts like an infinitely attracting/repelling saddle point; thus it loses its stability the moment $\mu$ has increased to zero. Since the derivative along the unstable manifold is positive, most points on the impact side will be mapped back to the impact side again, making nearly all points quickly vanish from the neighborhood of $O^*$. When the orientation is positive the system will quickly move away from the former stable periodic orbit.

$L_0$ and $L_1$ form a left-handed system, so the positive direction of $L_1$ points into the right-hand side of $L_0$. This is then called negative orientation. Now most points on the impact side are mapped to the non-impact side, so it is more likely that the motion stays close to $O^*$. A local co-ordinate system $(\chi, \psi)$ is chosen with origin at $O^*$, $L_1$ as the $\chi$-axis, and $L_0$ as the $\psi$-axis (keeping the directions). Now let the subscripts 0 and 1 denote a point and its image under the mapping $P$. Using the fact that $\chi=0$ (the curve $L_0$) is mapped to $\psi=0$ (the curve $L_1$), the mapping on the non-impact side is found to be

$$ \chi_1 = a_1\chi_0 + a_3\psi_0 + \mathcal{O}(\chi_0, \psi_0)^2, \quad \psi_1 = -a_2\chi_0[1 + \mathcal{O}(\chi_0, \psi_0)], \quad \chi_0 \ge 0. $$

By changing scales if necessary, the element $a_3$ can be made equal to 1. The elements $a_1$ and $a_2$ can then be expressed in terms of the two eigenvalues of the non-impact mapping, evaluated at $O^*$, $\lambda_1$ and $\lambda_2$: $a_1 = \lambda_1 + \lambda_2$, $a_2 = \lambda_1\lambda_2$. With this notation the mapping becomes

$$ \chi_1 = (\lambda_1 + \lambda_2)\chi_0 + \psi_0 + \mathcal{O}(\chi_0, \psi_0)^2, \quad \psi_1 = -\lambda_1\lambda_2\chi_0[1 + \mathcal{O}(\chi_0, \psi_0)], \quad \chi_0 \ge 0. \quad (6) $$

On the impact side, points will have gone through an extra low-velocity impact, described by a mapping of type $P_{IV} \circ P_{III} \circ P_{II}$. This will contribute an extra factor $(v_{out}/v_{in}) \times (dG/dv_{in})$ to the determinant of the Jacobian of the mapping. By l'Hôpital's rule,

$$ \lim_{v_{in} \to 0^{-}} (v_{out}/v_{in}) = -r, $$

so at $O^*$ the extra factor is $r^2$. Using this and the continuity of the mapping gives the mapping on the impact side as

$$ \chi_1 = B\sqrt{-\chi_0} + \psi_0 + \mathcal{O}(\sqrt{-\chi_0}, \psi_0)^2, \quad \psi_1 = \lambda_1\lambda_2 r^2 (-\chi_0) [1 + \mathcal{O}(\sqrt{-\chi_0}, \psi_0)], \quad \chi_0 \le 0, \quad (7) $$

where $B > 0$.

The next task is to establish whether $O^*$ is still stable. Since the mapping is orientation preserving, the eigenvalues $\lambda_1$ and $\lambda_2$ are either (a) both negative, (b) both complex or (c) both positive. Since the motion was assumed to be stable before coming to grazing impact, their absolute values are less than one. By taking a point on the impact side near $O^*$ with $\chi_0 = -\delta$ and iterating $P$ until the image returns to the impact side again, in cases (a) and (b) it is found that $\chi_n = \mathcal{O}(\sqrt{\delta})$ when returning, so points move away from $O^*$. In case (c) on the other hand, the iterates never return to the impact side, but instead converge to $O^*$ from the non-impact side. The former stable fixed point $O^*$ has thus become unstable except in case (c).
To summarize, a stable periodic orbit that comes to grazing impact will lose its stability unless the eigenvalues are positive and the orientation is negative. In this case, however, it will be shown that the motion stays near $O^*$ even when $\mu > 0$, and that the motion can be analyzed to a substantial extent.

## 4.2. THE MAPPING WHEN $\mu > 0$

When $\mu > 0$ there is no longer a fixed point on the non-impact side, but the non-impact side mapping is totally “unaware” of the existence of the line $L_0$ and the change in the mapping. This gives a way to define a virtual fixed point $O^*$ as the location of where a fixed point would have been, had not the change in the mapping been present. This virtual fixed point is located on the impact side of $L_0$. Denote the distance between $O^*$ and $L_0$ by $d$ (measured in $\chi$ units of the $(\chi, \psi)$ co-ordinate system defined when $\mu = 0$). Given that the crossing of $L_0$ by $O^*$ is transversal, $d$ is proportional to $\mu$ to the first order. The constants in the mappings (6) and (7) have changed by $\mathcal{O}(d)$ amounts, so $L_0$ is approximately described by a vertical line at a horizontal distance $d$ to the right $O^*$, and $L_1$ is approximately described by a horizontal line at a vertical distance $\lambda_1\lambda_2d$ below $O^*$. Now make an $\mathcal{O}(d)$ co-ordinate change from $(\chi, \psi)$ to $(\xi, \eta)$, making $O^*$ correspond to $(\xi, \eta) = (0, 0)$, $L_0$ to $\xi = d$, and $L_1$ to $\eta = -\lambda_1\lambda_2d$. In the new co-ordinates the mapping is

$$ \left\{ \begin{aligned} \xi_1 &= B\kappa_0 + \eta_0 + (\lambda_1 + \lambda_2)d + \mathcal{O}(d, \kappa_0, \eta_0)^2 \\ \eta_1 &= -\lambda_1\lambda_2d + \lambda_1\lambda_2r^2\kappa_0^2[1 + \mathcal{O}(d, \kappa_0, \eta_0)] \end{aligned} \right\}, \qquad \xi_0 \le d, \\ \left\{ \begin{aligned} \xi_1 &= (\lambda_1 + \lambda_2)\xi_0 + \eta_0 + \mathcal{O}(d, \xi_0, \eta_0)\mathcal{O}(\xi_0, \eta_0) \\ \eta_1 &= -\lambda_1\lambda_2d - \lambda_1\lambda_2(\xi_0 - d)[1 + \mathcal{O}(\xi_0, \eta_0)(1 + \mathcal{O}(d))] \end{aligned} \right\}, \qquad \xi_0 \ge d, $$

where $\kappa_0 \stackrel{\text{def}}{=} \sqrt{d - \xi_0}$. Upon noticing that $(0, 0)$ is mapped to a point with co-ordinates of $\mathcal{O}(\sqrt{d})$ it seems reasonable to make a uniform change of scale to make the co-ordinates describing the dynamics be of $\mathcal{O}(1)$. Therefore define

$$ \delta \stackrel{\text{def}}{=} \sqrt{d}/B, \quad X \stackrel{\text{def}}{=} \xi/B^2\delta, \quad Y \stackrel{\text{def}}{=} \eta/B^2\delta, \quad K \stackrel{\text{def}}{=} \sqrt{1 - X/\delta}. $$

Here, $\delta$ is a new parameter, proportional to $\sqrt{\mu}$ to first order, and $X$ and $Y$ are the new co-ordinates. The mapping becomes

$$ \left\{ \begin{aligned} X_1 &= K_0 + Y_0 + (\lambda_1 + \lambda_2)\delta + \mathcal{O}(\delta)\mathcal{O}(\delta, X_0, Y_0)^2 \\ Y_1 &= -\lambda_1\lambda_2\delta + \lambda_1\lambda_2r^2\delta K_0^2[1 + \mathcal{O}(\delta)\mathcal{O}(\delta, X_0, Y_0)] \end{aligned} \right\}, \qquad X_0 \le \delta, \\ \left\{ \begin{aligned} X_1 &= (\lambda_1 + \lambda_2)X_0 + Y_0 + \mathcal{O}(\delta)\mathcal{O}(\delta, X_0, Y_0)\mathcal{O}(X_0, Y_0) \\ Y_1 &= -\lambda_1\lambda_2\delta - \lambda_1\lambda_2(X_0 - \delta)[1 + \mathcal{O}(\delta)\mathcal{O}(X_0, Y_0)] \end{aligned} \right\}, \qquad X_0 \ge \delta. \quad (8) $$

In the case (c) it can be seen that if co-ordinates are of $\mathcal{O}(1)$ on the non-impact side, and of $\mathcal{O}(\delta)$ on the impact side, they will remain so under iteration of the mapping. This shows that the scale chosen is the correct one when describing the dynamics for $\mu > 0$. Neglecting higher order terms, a truncated mapping is obtained:

$$ \left\{ \begin{aligned} X_1 &= K_0 + Y_0 + (\lambda_1 + \lambda_2)\delta \\ Y_1 &= -\lambda_1\lambda_2\delta[1 - r^2 K_0^2] \end{aligned} \right\}, \qquad X_0 \le \delta, \\ \left\{ \begin{aligned} X_1 &= (\lambda_1 + \lambda_2)X_0 + Y_0 \\ Y_1 &= -\lambda_1\lambda_2 X_0 \end{aligned} \right\} \qquad X_0 \ge \delta. \quad (9) $$

When $\delta$ becomes smaller, the real mapping (8) will approach the truncated mapping (9), so conclusions drawn from the truncated mapping should apply to the real one for $\delta$ sufficiently small.
## 5. STUDY OF THE TRUNCATED MAPPING

In the following, lower case letters are used for the co-ordinates.

$$ egin{cases} x_1 = k_0 + y_0 + (\lambda_1 + \lambda_2)\delta \\ y_1 = -\lambda_1\lambda_2\delta[1 - r^2 k_0^2] \end{cases}, \quad x_0 \le \delta \quad (I), \\ 
\begin{cases} x_1 = (\lambda_1 + \lambda_2)x_0 + y_0 \\ y_1 = -\lambda_1\lambda_2x_0 \end{cases}, \quad x_0 \ge \delta \quad (II). \tag{10} $$

The non-linear (impact) region is denoted by I, and the linear region by II. By assumption, the eigenvalues satisfy $1 > \lambda_1 > \lambda_2 > 0$, and $\delta > 0$.

In region II the mapping is given by

$$ \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = L \begin{pmatrix} x_0 \\ y_0 \end{pmatrix}, \quad \text{where} \quad L = \begin{pmatrix} \lambda_1 + \lambda_2 & 1 \\ -\lambda_1\lambda_2 & 0 \end{pmatrix}. $$

The eigenvectors corresponding to $\lambda_1$ and $\lambda_2$ are

$$ \zeta_1 = \begin{pmatrix} 1 \\ -\lambda_2 \end{pmatrix} \quad \text{and} \quad \zeta_2 = \begin{pmatrix} 1 \\ -\lambda_1 \end{pmatrix}. $$

### 5.1. EXISTENCE OF A TRAPPING REGION

Consider the line segment $y = -\lambda_2 x$ for $x \ge 0$. This line is in the direction of $\zeta_1$, so the part of the line in region II is mapped onto the line. This makes the line act like a boundary in region II. The form of the mapping in region I shows that if $x_0 \ge 0$, then $y_1 \le 0$ (since $r \le 1$). Furthermore, in region II line $y=0$ is mapped below itself but above $y = -\lambda_2 x$. The sector between the lines $y = -\lambda_2 x$ for $x \ge 0$ and $y = 0$ for $x \ge 0$ is consequently mapped inside itself in region II. Examining the mapping in region I of the triangle bounded by these two lines and the line $x = \delta$, reveals that the corner $A_0=(0,0)$ maps to $A_1=(1+(\lambda_1+\lambda_2)\delta, -\lambda_1\lambda_2\delta(1-r^2))$, which is inside the sector, while the other two corners, by the continuity of the mapping and the results found for region II, also map inside the sector. This establishes a trapping sector of infinite size. By cutting it off with a line between the points $B_0=(c,0)$ and $C_0=(c\lambda_1/(1-\lambda_2), -c\lambda_1\lambda_2/(1-\lambda_2))$, where c is a constant slightly greater than one and chosen to make the image of $A_0$ fall inside the triangle, a trapping region of size $O(1)$ is obtained.

An example of the trapping region and its image at parameter values $\lambda_1 = 0.5$, $\lambda_2 = 0.2$, $r=1$ and $\delta = 0.1$ is given in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. The trapping region.</figcaption></figure>
## 5.2. STABILITY OF PERIODIC POINTS

The Jacobians in the two regions are

$$J_I = \begin{pmatrix} -1/(2\delta k_0) & 1 \\ -r^2\lambda_1\lambda_2 & 0 \end{pmatrix} \quad \text{and} \quad J_{II} = L = \begin{pmatrix} \lambda_1 + \lambda_2 & 1 \\ -\lambda_1\lambda_2 & 0 \end{pmatrix}.$$

All points in region II are eventually mapped into region I, so there can be no periodic point the orbit of which is entirely located in region II. In region II, where both eigenvalues are less than one, an area element will be compressed in both the eigenvalue directions by the mapping. In region I, the presence of the entry $1/(2\delta k_0)$ in the Jacobian will have the effect of a very large stretching and compression of the area element. The most favorable situation for stability in a periodic orbit is to have the iterates mainly stay in region II, with occasional visits to region I. If the stability acquired in region II outweighs the unstability acquired in region I, the point is stable. To determine whether stable periodic points are possible at all, an estimate of the eigenvalue of largest absolute value for any periodic point will be established.

For any steady motion in the trapping region, the value of $k$ obeys $0 \le k \le \sqrt{1 - \lambda_1} < 1$ when iterates are in region I. The dominating element in $J_I$ when $\delta$ is small is $1/(2\delta k_0)$, which is larger than $1/(2\delta)$. A point making $n$ consecutive visits in region II acquires a Jacobian $L^n$, that can be evaluated explicitly as

$$L^n = \frac{\lambda_1^n}{\lambda_1 - \lambda_2} \begin{pmatrix} \lambda_1[1 - (\lambda_2/\lambda_1)^{n+1}] & 1 - (\lambda_2/\lambda_1)^n \\ -\lambda_1\lambda_2[1 - (\lambda_2/\lambda_1)^n] & -\lambda_2[1 - (\lambda_2/\lambda_1)^{n-1}] \end{pmatrix}.$$

Now define a mapping $F$ taking points in region I to their first return to region I. The Jacobian $JF$ of this mapping will be $L^n J_I$, where $n$ is the number of iterations spent in region II. The structure of this matrix is

$$JF = L^n J_I = -A \begin{pmatrix} 1 & \beta \\ \alpha & \gamma \end{pmatrix},$$

where $\beta$ and $\gamma$ are $\mathcal{O}(\delta)$, and $\alpha$ is $\mathcal{O}(1)$; so, when seeking the largest eigenvalue, it can be approximated as

$$L^n J_I \approx -A \begin{pmatrix} 1 & 0 \\ \alpha & 0 \end{pmatrix}. \qquad (11)$$

If $k_0 = \mathcal{O}(\delta)$, then $x_1 = \mathcal{O}(\delta)$, so the point will return to region I in a few iterations, making $L^n = \mathcal{O}(1)$ and $A = \mathcal{O}(1/(\delta k_0)) = \mathcal{O}(\delta^{-2})$. If $k_0 = \mathcal{O}(1)$, then the number of iterations $n$ must be estimated.

Having $k_0 = \mathcal{O}(1)$ leads to

$$\begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} k_0 \\ 0 \end{pmatrix} + \mathcal{O}(\delta), \quad \text{so} \quad \begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} = L^n \begin{pmatrix} k_0 \\ 0 \end{pmatrix} + L^n \mathcal{O}(\delta).$$

By the definition of $n$, $\lambda_1 \delta \le x_{n+1} \le \delta$, so $L^n \mathcal{O}(\delta)$ is $\mathcal{O}(\delta)^2$. Using $x_{n+1} = (L^n)_{11} k_0 + \mathcal{O}(\delta)$, an estimate of $n$ is given by

$$\lambda_1 \le \frac{\lambda_1^{n+1} [1 - (\lambda_2/\lambda_1)^{n+1}]}{\delta(\lambda_1 - \lambda_2)} \le k_0 \le 1, \qquad (12)$$

to $\mathcal{O}(1)$. Evaluating (11) explicitly gives

$$A = \frac{\lambda_1^{n+1} [1 - (\lambda_2/\lambda_1)^{n+1}]}{2k_0 \delta(\lambda_1 - \lambda_2)}, \qquad (13)$$
and when combining expressions (12) and (13), $\lambda_1 \le 2k_0^2A \le 1$. Since $k_0^2 \le 1 - \lambda_1$, this gives

$$A \ge \lambda_1/2(1 - \lambda_1). \quad (14)$$

When multiplying two Jacobians of the form (11), the result is a matrix of the same form, with $A = -A_1A_2$. This shows that if every $A$ is greater than one, then there can be no stable periodic points for the mapping $F$, and consequently none for the truncated mapping itself. From (14), this happens when $\frac{2}{3} < \lambda_1$. The above calculations can be viewed as an estimate of a lower bound of the largest Lyapunov exponent, showing that when $\frac{2}{3} < \lambda_1$, the largest Lyapunov exponent must be positive. Since this estimate is independent of $\delta$ as long as $\delta$ is small, the conclusion is that if the largest eigenvalue is greater than $\frac{2}{3}$, there are no stable periodic points in the trapping region for an interval $0 < \delta < \delta'$.

### 5.3. DEPENDENCE ON THE PARAMETER $\mu$

As the parameter $\mu$ increases from zero, the dynamics of the system will in general undergo a complex series of bifurcations, accumulating on $\mu = 0$, which can be studied by the truncated mapping (10). The co-ordinates of first return are given to $O(\delta)$ by

$$\begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} \approx \frac{\lambda_1^{n+1}}{\lambda_1 - \lambda_2} \begin{pmatrix} 1 - (\lambda_2/\lambda_1)^{n+1} \\ -\lambda_2[1 - (\lambda_2/\lambda_1)^n] \end{pmatrix} k_0,$$

and if $n$ is large,

$$\begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} \approx \frac{\lambda_1^{n+1}}{\lambda_1 - \lambda_2} \begin{pmatrix} 1 \\ -\lambda_2 \end{pmatrix} k_0,$$

since $\lambda_2/\lambda_1 < 1$. If $n$ is large then (12) simplifies to

$$\lambda_1 \le \frac{\lambda_1^{n+1}}{\delta(\lambda_1 - \lambda_2)} k_0 \le 1,$$

so when $\delta$ is small, $n$ is indeed large, unless $k_0$ is $O(\delta)$. If $\delta$ changes by a factor $\lambda_1$ then $n$ changes by one. This means that except in the small region where $k_0$ is $O(\delta)$, the mapping $F$ taking points in region $I$ to their first returns will come back the same when $\delta$ has changed by a factor $\lambda_1$. The same goes for the dynamics, so if some type of motion is observed when $\delta = \delta^*$, then a similar type of motion (with the number of visits in region II between visits in region I increased by one) is observed when $\delta = \lambda_1\delta^*$. Thus similar types of motion tend to re-appear at a sequence of values $\delta_1, \delta_2, \dots$, where

$$\lim_{k \to \infty} \frac{\delta_{k+1} - \delta_k}{\delta_k - \delta_{k-1}} = \lambda_1,$$

or, since $\delta$ is proportional to $\sqrt{\mu}$ to first order,

$$\lim_{k \to \infty} \frac{\mu_{k+1} - \mu_k}{\mu_k - \mu_{k-1}} = \lambda_1^2. \quad (15)$$

In Figure 8 is shown a bifurcation diagram obtained by numerical simulation of the first return mapping. The parameters are $\lambda_1 = 0.5$, $\lambda_2 = 0.2$ and $r = 1$. Only the $x$ co-ordinate is shown. The similarity between motion at values of $\delta$ that differ by a factor $\lambda_1$ is clearly seen.

In the co-ordinates used in the truncated mapping (10) the dynamics move over a region of size $O(1)$. Since these co-ordinates are scaled by $1/\delta$, the size of the region in the usual co-ordinates $(x, v)$ is $O(\delta)$, or $O(\sqrt{\mu})$ (for positive $\mu$).
<figure><img src="image_10.png" /><figcaption>Figure 8. Bifurcation diagram for the first return mapping.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. One thousand iterations of the truncated mapping in a chaotic region.</figcaption></figure>

## 5.4. NUMERICAL EXAMPLE OF THE DYNAMICS

The dynamics of the truncated mapping for parameter values $\lambda_1=0.5$, $\lambda_2=0.2$, $r=1$ and $\delta=0.01$ is shown in Figure 9. $\delta=0.01$ corresponds to $\log(\delta)/\log(\lambda_1)=6.65$, so it is in one of the chaotic regions in Figure 8. A striking feature is that iterates line up in curve-like objects. This is due to the very strong expansion/compression in region I. In the figure showing the trapping region this effect is easy to see. This nearly one-dimensional feature of the dynamics becomes more pronounced as $\delta$ becomes smaller, and the same results regarding stability and $\delta$ dependence are found for a one-dimensional mapping consisting of a square root part and a linear part pasted together. Another feature which is typical of grazing impact is the corner of the attractor near (0,0). This is the same corner effect as was shown in Figure 6(b).

## 6. USING THE THEORY IN NUMERICAL SIMULATIONS

When applying the theory to a specific system, the system must first be put into the standard form with functions A and G and a fixed impact boundary. If the impact boundary is moving, this motion is subtracted from the displacement so the impact
boundary in the new co-ordinates is at rest. Next, one wants to locate periodic orbits
having a grazing impact. If a period and a plane of constant phase is chosen, this can
be achieved by taking a point in Σº with A > 0 and following the motion forwards and
backwards in time to the intersection with the constant phase plane. The points of
intersection for the forward and backward trajectories must then coincide, since the
motion is to be periodic. This condition gives two equations (equality in displacement
and velocity in the constant phase plane) for the starting phase in Σº and the system
parameters. If all but one of the system parameters are kept fixed, the two equations can
then be solved for the starting phase and the remaining system parameter by a numerical
root finding algorithm. If derivatives with respect to the system parameter are needed for
the algorithm, they can be calculated through the trick of augmenting the system (2) with
the system parameter as a new dependent variable, with a zero time derivative. The
Jacobians equivalent to $JP_{I-IV}$ for the augmented system are easily derived, by using the
same kind of derivation. The local mappings are composed to give a mapping from Σ to
the constant phase plane backwards in time, and a mapping from Σ to the constant phase
plane forwards in time. Subtracting one of the mappings from the other gives a mapping
which must be zero at a grazing fixed point. When a grazing fixed point is located, the
Poincaré mapping near it is given by the inverse of the backward mapping composed
with the forward mapping, and the two eigenvalues can thus be calculated. The directions
of the curves $L_0$ and $L_1$ are given by the image of a vector tangent to Σº in the direction
of decreasing phase under the backward and forward mappings respectively. The orienta-
tion condition is then determined by whether the two image vectors form a right-hand
pair or not. Thus all quantities needed in the theory are determined by the first derivative
information given in $JP_{I-IV}$, and can easily be computed.

7. RESULTS AND DISCUSSION

In this paper the motion of single-degree-of-freedom, periodically forced impact oscil-
lators has been discussed, and in particular the consequences when a stable periodic orbit
comes to a grazing impact under the control of a single parameter. Unless one orientation
condition (orientation being negative) and one condition on eigenvalues (eigenvalues
being positive) are satisfied, the orbit becomes unstable, so motion will go elsewhere.
However, if these conditions are satisfied, the motion will stay trapped near the periodic
orbit and the dynamics in this case include either orbits of arbitrarily long periods,
possibly interleaved by non-periodic motion, or non-periodic motion for an interval of
parameter values down to the point at which grazing impact appeared. If arbitrarily long
periods appear, the parameter values of similar bifurcations form a geometric sequence
such as expression (15). The linear size of the trapping region scales as the square root
of the normalized bifurcation parameter.

These results have been derived for a situation in which the impact is considered to
be instantaneous. In this case it was shown that a square root singularity existed at lines
of points leading to grazing impact, and a truncated map was derived that becomes a
better and better approximation as a bifurcation parameter goes to zero. If a model with
more structure in the impact is used, for example a stiff spring model, then the notion
of impact may become less distinct. There may be no true singularity at the lines leading
to grazing impact, and the lines themselves can be hard to define. Then the truncated
mapping will not be a good approximation for very small values of the parameter, since
impacts are so soft that the system will hardly be affected by them. On the other hand,
as the bifurcation parameter grows larger, the truncated mapping will lose accuracy, being
a truncated power series. For a more realistic model, one will have to investigate whether
these two regions of non-validity overlap, or if there is an interval where this theory is applicable. The shorter the time of impact is, the more favorable is the situation for finding such an interval.

A "fingered" structure of chaotic attractors in constant phase Poincaré sections has been observed earlier in impact oscillators: see, for example in Figure 10 of Shaw's paper [9]. The forms of the attractors are somewhat similar to the form of the attractor in Figure 9. In both cases only a minor part of the attractor impacts within a period of the force, forming a finger-like part of the attractor which is then mostly left undisturbed since it takes several periods until the points in that part impacts again. This phenomenon creates a sequence of fingers, each mapped onto the next during a period of the force. In the former case, however, when looking closely the structure of the attractor can be seen to be far more complex, showing clear fractal layers. The size of the attractor is not small, and it is not born out of a stable periodic point, but instead out of a period-doubling sequence. The number of fingers is also roughly proportional to the ratio of driving frequency to natural frequency of the oscillator, while in the grazing impact attractor the number of fingers grows steadily as the bifurcation parameter goes to zero. The theory presented here then does not claim to explain these types of attractors. There is, however, a geometrical feature of these attractors that can be attributed to low-velocity impact. In the "tips" of the fingers there is a curled-up structure that ends in a sharp corner, and which seems almost to be disconnected from the rest of the attractor. The curling is a result of a series of low-velocity impacts with negative acceleration, and when the acceleration subsequently becomes positive the structure leaves the neighborhood of the impact boundary through low-velocity impact with positive acceleration. A line of points will have zero-velocity impacts and at that line a corner is formed in the attractor, and the large stretching near this line is the cause of the apparent non-connectedness.

The non-periodic motion that can occur right after the bifurcation is not very dramatic in character. When the periodic orbit passes through the bifurcation, the orbit starts to jitter around the former periodic orbit, and the amplitude of the jitter grows after the bifurcation. This local character of the non-periodic motion is what makes a theory based on series expansions work and give precise results, but, as it is a local theory, the results are guaranteed to be valid only near the point of bifurcation. Nonetheless, the results give a mechanism for the onset of chaotic motion in impacting oscillator problems, and gives criteria that are easy to check for a given system.

The cause of the phenomena described here is the fact that the flow is tangent to a boundary in state space where the system characteristics change. Here a simple mechanical system has been studied, but in general, the tangency of a flow to a boundary in state space could be the cause of similar effects in different systems.

## ACKNOWLEDGMENTS

I would like to thank Professor Martin Lesser for helpful comments, and the Volvo Research Foundation for supporting this work.

## REFERENCES

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

2. S. W. SHAW and P. J. HOLMES 1983 *Journal of Applied Mechanics* 50, 849-857. A periodically forced impact oscillator with large dissipation.

3. S. W. SHAW 1985 *Journal of Applied Mechanics* 52, 453-458. Dynamics of harmonically excited systems having rigid amplitude constraints, Part 1—Subharmonic motions and local bifurcations.
4. S. W. SHAW 1985 *Journal of Applied Mechanics* **52**, 459-464. Dynamics of harmonically excited systems having rigid amplitude constraints, Part 2—Chaotic motions and global bifurcations.

5. G. S. WHISTON 1987 *Journal of Sound and Vibration* **115**, 303-319. The vibro-impact response of a harmonically excited and preloaded one-dimensional linear oscillator.

6. G. S. WHISTON 1987 *Journal of Sound and Vibration* **118**, 395-429. Global dynamics of a vibro-impacting linear oscillator.

7. F. C. MOON 1987 *Chaotic Vibrations*. New York: John Wiley.

8. J. GUCKENHEIMER and P. HOLMES 1986 *Non-linear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*. New York: Springer-Verlag, second edition.

9. S. W. SHAW 1985 *Journal of Sound and Vibration* **99**, 199-212. Forced vibrations of a beam with one-sided amplitude constraint: theory and experiment.