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

HAL
open science

# Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators

Mats Fredriksson, Arne Nordmark

► To cite this version:

Mats Fredriksson, Arne Nordmark. Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1997, 10.1098/rspa.1997.0069. hal-01297285

HAL Id: hal-01297285

https://hal.science/hal-01297285

Submitted on 4 Apr 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 - ShareAlike 4.0 International License
# Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators

BY MATS H. FREDRIKSSON AND ARNE B. NORDMARK

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

The transition from stable periodic non-impacting motion to impacting motion is analysed for a mechanical oscillator. By using local methods, it is shown that a grazing impact leads to an almost one-dimensional stretching in state space. A condition can then be formulated, such that a grazing trajectory will be stable if the condition is fulfilled. If this is the case, the bifurcation will be continuous and the motion after the bifurcation can be understood by a one-dimensional mapping. This mapping is known to exhibit chaotic solutions as well as arbitrary long stable cycles, depending on parameters.

## 1. Introduction

Impact oscillators are examples of systems that are inherently strongly nonlinear. As such, they have proven to be a challenge to analyse and understand. This and the many practical applications of impact problems, have led to extensive studies, using the progress in dynamical systems theory in recent years.

One peculiarity in the dynamics of impact oscillators comes from grazing or zero velocity impacts. Previous investigations have mainly dealt with one degree of freedom systems. Studying such a system, Shaw & Holmes (1983) noted that a zero velocity impact gives a singularity in the derivative of the Poincaré mapping. This lack of smoothness affects the dynamics in a profound way. Whiston (1992) analysed consequences for global dynamics. Nordmark (1991) studied the bifurcation of a stable periodic orbit when it is brought to a grazing impact by a change in a single parameter. In the generic case stability is lost, but Ivanov (1993) showed that there are cases when it is possible for the periodic orbit to be stable after the bifurcation. Further analysis of grazing bifurcations can be found in Foale & Bishop (1992), Budd & Dux (1994) and in Ivanov (1994). Effects similar to those found in Nordmark (1991) have been observed in systems with more than one degree of freedom, and not necessarily in systems with explicit time dependence. Knudsen et al. (1992) studied the dynamics of a model of a railway wheelset, where impact between the wheelset and the rail was included. On p. 461, figure 4, a bifurcation diagram is shown where stable periodic motion bifurcates to chaotic motion as the wheelset starts to impact. Close to the bifurcation point no periodic windows are visible, which is one of the possible scenarios described in Nordmark (1991).

The purpose of this paper is to study the grazing bifurcation of a stable periodic motion in a quite general class of mechanical systems. The first section of the paper reviews the formulation of equations of motion for constrained systems of rigid bodies
according to Kane. It is briefly discussed how this method can be used in impact problems where the duration of the impact is negligible. The system to investigate is then taken to be described by a differential equation, with generalized coordinates and speeds as coordinates on the state space, and the motion is restricted to take place on one side of a rigid boundary. When a part of the system impacts with the boundary, an instantaneous change of generalized speeds is assumed to take place. We emphasize a geometrical viewpoint, where the boundary is treated as a submanifold of the state space. The analysis then focuses on the problem of understanding the dynamics close to a grazing trajectory, close to the boundary. By series expansions an algebraic expression is obtained, which can be combined with the local dynamics of the non-impacting motion to give a stability criterion. Consequences for the Poincaré mapping technique are analysed. If the stability criterion is fulfilled, the motion after bifurcation can be understood. The dynamics will be similar to the one degree of freedom case. Finally a numerical example is studied.

## 2. Equations of motion and impacts in complex systems

Different methods to obtain equations of motion for complex mechanical systems have been proposed. The fundamental equation of mechanics, Newton's second law, has the disadvantage that, for a constrained system, the constraint forces have to be eliminated to get a minimal set of equations that describes the system. By far the most well known method that does not suffer from this problem is to use the Lagrangian equations. Here the method of Kane & Levinson (1985) will be briefly outlined, using the interpretation by Lesser (1995). The impact problem can be treated in the same spirit, formally giving an impact law.

### (a) Equations of motion

The positions and attitudes of a set of $k$ constrained bodies are given by $n$ generalized coordinates $q_i$, $i = 1, \dots, n$, and we denote

$$ q = \begin{pmatrix} q_1 \\ \vdots \\ q_n \end{pmatrix}. \qquad (2.1) $$

Differentiating $q$ with respect to time, $dq/dt = \dot{q}$, we get the generalized velocity. By the transformation

$$ u = A(q, t)\dot{q} + a(q, t), \qquad (2.2) $$

where $A$ is an invertible $n \times n$ matrix, we introduce a new velocity $u$. This velocity is called the generalized speed in Kane & Levinson (1985). By solving for $\dot{q}$ one obtains

$$ \dot{q} = B(q, t)u + b(q, t), \qquad (2.3) $$

where $B = A^{-1}$ and $b = -A^{-1}a$. Equation (2.3) is called a kinematic differential equation. The centre of mass velocity and the angular velocity of a body can then be written as a linear combination of vectors with either $\dot{q}_i$ or $u_i$ as coefficients (plus another term if there is an explicit time dependence or if $b \neq 0$). Lesser (1995) interpreted these vectors as tangent vectors to the instantaneous configuration manifold and showed that Kane's method is equivalent to a projection procedure of the
Newton–Euler equations, using the tangent vectors. This automatically eliminates the constraint forces and the resulting equations are on the form

$$ M \dot{u} = f, \qquad (2.4) $$

where $M = M(q, t)$ is an $n \times n$ matrix and $f = f(q, u, t)$.

There are several advantages of this method. Compared to that of Lagrange, there is no need to form the kinetic energy and then proceed by differentiation. The introduction of generalized speeds can sometimes significantly reduce difficulties in problem solving. This method also handles (simple) non-holonomic problems with ease, as they can be viewed as systems with a knowledge of some $u_j: u_j = 0, j = m+1, \dots, n$. Indeed, there will not be any need to exclude non-holonomic systems. We will thus assume that the generalized speed is

$$ u = \begin{pmatrix} u_1 \\ \vdots \\ u_m \end{pmatrix}, \qquad (2.5) $$

where $m \le n$. The matrix $B$ of the kinematic differential equation (2.3) will then be an $n \times m$ matrix.

### (b) Impacts

Collisions between stiff bodies involve a change of velocities in a short period of time. Thus, in the limit of infinitely stiff bodies one might anticipate that an instantaneous change in velocity will give a reasonable model. The impact approximation is then to introduce a discontinuity in velocities (speeds) at the time $t^*$ of impact, while the generalized coordinates vary continuously. By labelling the generalized speed just after and just before impact $u_a$ and $u_b$, we assume that an impact law can be written as $u_a = G(u_b, q^*)$, where $q^* = q(t^*)$. We can formally obtain the function $G$ by integrating equation (2.4),

$$ \int_{t^* - \Delta t}^{t^* + \Delta t} M \dot{u} \, dt = \int_{t^* - \Delta t}^{t^* + \Delta t} f \, dt. \qquad (2.6) $$

$M$ will be slowly varying during the impact. We let $\Delta u = u_a - u_b$, and by introducing the impulse

$$ I = \int_{t^* - \Delta t}^{t^* + \Delta t} f \, dt, \qquad (2.7) $$

we have

$$ M(q^*, t^*) \Delta u = I, \qquad (2.8) $$

and thus

$$ u_a = u_b + M^{-1}(q^*, t^*)I. \qquad (2.9) $$

Of course, the impulse must be determined by physical reasoning and by some kind of model of the impact (in finite time). This paper will not deal with this problem, it will be assumed that this can be done.

## 3. System, state space and local geometry

The system to be considered is a set of constrained rigid bodies, described by a choice of generalized coordinates, $q$, and generalized speeds, $u$. As seen in the previous
section, the motion is then given by equation (2.3) and (2.4). These equations might have an explicit time dependence. Only periodic time dependence will be of interest here. Whether the system is autonomous or periodically driven is of no importance for the derivations and conclusions below. In the case of a system without explicit time dependence the state space is $n+m$-dimensional, and we let

$$x = \begin{pmatrix} q \\ u \end{pmatrix}. \qquad (3.1)$$

The smooth evolution of the system is determined by a differential equation,

$$\dot{x} = F(x). \qquad (3.2)$$

In the periodically driven case we convert the system to the autonomous form (3.2) in the standard way, by extending $x$ with a phase angle proportional to time. The phase angle is unaffected by impacts, so we can view it as a generalized coordinate. Hence, there is no reason to treat this case separately. The formal solution of the differential equation is $X = X(x_0, t)$, which satisfies

$$\frac{\partial X}{\partial t}(x_0, t) = F(X(x_0, t)) \qquad (3.3)$$

and

$$X(x_0, 0) = x_0. \qquad (3.4)$$

The formal solution will define a one parameter mapping $\Phi_t$ by

$$\Phi_t : x_0 \mapsto X(x_0, t). \qquad (3.5)$$

In the model we also include impacts by assuming that the motion of the bodies of the system is restricted to take place on one side of a rigid boundary. An analysis of the global dynamics will have to deal with the possibility of several contact conditions, which almost certainly will lead to a complex behaviour. However, our main concern here will be local dynamics, which allows us to avoid some of the problems involved. We introduce a bifurcation parameter $\mu$, and we assume that for an interval of negative $\mu$ values there is a stable periodic non-impacting solution. The dependence of $\mu$ is such that when $\mu$ is increased, the periodic solution will come closer to the boundary. When $\mu = 0$ a body will make contact with the boundary somewhere along the trajectory of the periodic solution, and since this contact will take place at certain values of the generalized coordinates and speeds, it corresponds to a point $O^*$ in state space. This is assumed to be the only point where the separation between the periodic solution and the boundary vanishes. It is also assumed that if it were not for the possibility of impacts, the stable periodic solution would exist for a parameter interval of positive $\mu$ values.

We define a function $h = h(x)$ in a neighbourhood $\mathcal{B}$ of $O^*$ to be positive and equal to the minimal distance from the boundary to the impacting body when $x$ is on the side where motion takes place. For points which corresponds to a penetration of the boundary we take $h(x) < 0$. This function only depends on the generalized coordinates. One will have contact when $h(x) = 0$, and this can be viewed as an equation for a local surface in state space, dividing $\mathcal{B}$ in two parts. This surface, the impact surface, is clearly of crucial importance and we use the notation

$$\Sigma = \{x \in \mathcal{B} : h(x) = 0\}. \qquad (3.6)$$
<figure><img src="image_3.png" /><figcaption>Figure 1. Grazing trajectory.</figcaption></figure>

If the formal solution is inserted in $h$, we will have a new function

$$H(x_0, t) = h(\Phi_t(x_0)), \qquad (3.7)$$

which describes the minimum distance to the boundary as a function of time and initial condition. It will be convenient to define a function $v = v(x)$ on $\Sigma$:

$$v = \frac{\partial H}{\partial t}(x_0, 0), \quad x_0 \in \Sigma. \qquad (3.8)$$

This allows us to divide the impact surface into different subsets. We define

$$\Sigma^+ = \{x \in \Sigma : v > 0\}, \qquad (3.9)$$

$$\Sigma^0 = \{x \in \Sigma : v = 0\}, \qquad (3.10)$$

$$\Sigma^- = \{x \in \Sigma : v < 0\}. \qquad (3.11)$$

The impact of the periodic solution when $\mu = 0$ is assumed to be a grazing impact, which we define as an impact with $v = 0$ and $\partial^2 H/\partial t^2 = A_g$, where $A_g$ is positive and of order $O(1)$. Thus part of the system is brought to contact with the boundary, but is immediately pulled back.

Finally an impact law is used for the impact process. With the definitions above it can be viewed as a mapping $G: \Sigma^- \to \Sigma^+$, that leaves the position coordinate unchanged. For a grazing impact the impact law is assumed to be reduced to the identity mapping. This assumption is valid for simple models that involve motion in only one dimension, but for more complicated systems this might not be true.

In figure 1 the grazing trajectory is indicated, the point $O_1$ is mapped by the flow to $O^* \in \Sigma^0$. Here another geometrical feature is apparent: choosing a point near $O_1$ will either result in an impact or not when the trajectory is followed. Thus a neighbourhood of $O_1$ will be divided in two parts, separated by points that will be mapped to $\Sigma^0$. As we will see, this will be of significant importance for the dynamics.

## 4. Local mappings

The standard approach to the stability problem is to linearize the map $\Phi$. For a point that undergoes a grazing impact this does not work for the obvious reason that
<figure><img src="image_4.png" /><figcaption>Figure 2. The C mapping, $P \mapsto C(P)$.</figcaption></figure>

only for non-impacting points will the result be correct. For impacting points an error will be introduced due to the fact that the dynamics of the impact is absent. To solve this problem a discontinuity bypass mapping is introduced. This will be a crucial step towards an understanding of the grazing bifurcation, since local flow mappings are well understood and derivatives are easily calculated from the variational equations.

As in figure 1, let $O^* \in \Sigma^0$ be the image of $O_1$, using $\Phi_{t_1}$, where $t_1$ is the time of flight from $O_1$ to $O^*$. In the same way we let $O_2$ be the image of $O^*$, and the time of flight from $O^*$ to $O_2$ is $t_2$. We now construct a mapping $C$, defined for points near $O^*$, that will enable us to write a locally valid mapping from a neighbourhood of $O_1$ to a neighbourhood of $O_2$ as $\Phi_{t_2} \circ C \circ \Phi_{t_1}$.

If the full mapping is to be a series of compositions starting with $\Phi_{t_1}$, one should note that we then must let the flow pass through $\Sigma$. This is because the mapping $\Phi_{t_1}$ disregards the boundary. By doing so the model is violated and the mapping $C$ must correct this. The points close to $O^*$ that will need this correction are the points lying on trajectories that cross $\Sigma^-$ if the flowline is followed backwards or forwards for a short time. Starting at a point with coordinates $x, H(x, t)$ will tell us how the distance to the boundary changes as $t$ is varied. By our assumptions it is clear that $H$ will have a local minimum for a (small) time $\tau = \tau(x)$. We can then introduce a function $\psi$ by

$$ \psi(x) = H(x, \tau(x)). \qquad (4.1) $$

This will help us to distinguish impacting points with $\psi(x) \le 0$, from non-impacting points. For impacting points we define the mapping $\Phi_{t_C}$, where $t_C$ is defined as the time of flight from the impacting point to the crossing point of $\Sigma^-$. By the definition of $\Sigma$ we have that $t_C$ must fulfil $H(x, t_C) = 0$. This map takes impacting points to the impact surface. Hence a natural next step is to use the impact mapping to include the dynamics of the impact. The mapping $G \circ \Phi_{t_C}$ will not do as a candidate for $C$, since it does not take place in zero time. This is achieved by composition with
$\Phi_{-t_C}$, so that the full discontinuity bypass mapping will become

$$C = \begin{cases} \Phi_{-t_C} \circ G \circ \Phi_{t_C}, & \text{if } \psi \le 0, \\ \text{identity}, & \text{if } \psi > 0. \end{cases} \qquad (4.2)$$

Finally we note that a point close to $O_1$ that will undergo an impact, will have a time of flight $t_1 + t_C$ to reach $\Sigma^-$. By the properties of $\Phi$ we have $\Phi_{t_1+t_C} = \Phi_{t_C} \circ \Phi_{t_1}$, thus

$$\Phi_{t_2-t_C} \circ G \circ \Phi_{t_1+t_C} = \Phi_{t_2} \circ C \circ \Phi_{t_1}, \qquad (4.3)$$

which is the desired form.

### (a) Algebraic expressions

We now consider series expansions of $C$. It will be assumed that local coordinates have been introduced to make $O^*$ correspond to $x^* = 0$. The notation $O(x, t)^p$ will be used to indicate terms of total order $p$ in $x$ and $t$.

Let us start by analysing the impact law. To generate a set of $n+m-1$ independent coordinates $\rho$ on $\Sigma$ we note that $h(x) = 0$ can be used to write, say, $q_1 = q_1(q_2, \dots, q_n)$ by using the implicit function theorem. We define our new coordinates so that $v$ is used instead of $u_1$ by letting $\rho_1 = v(q_1(q_2, \dots, q_n), q_2, \dots, u)$ and taking the remaining coordinates to be $q_2, \dots, q_n, u_2, \dots, u_m$. Equation (2.9) then indicates that

$$\Delta u = w(\rho). \qquad (4.4)$$

By assumption, $\rho_1 = 0$ will give the identity mapping, thus $\Delta u = 0$ when $\rho_1 = 0$. If we let

$$\xi_u = -\frac{\partial w}{\partial \rho_1}(0), \qquad (4.5)$$

the expansion of $w$ must be

$$\Delta u = -\rho_1(\xi_u + \mathcal{O}(\rho)). \qquad (4.6)$$

We extend $\xi_u$ to the vector $\xi$,

$$\xi = \begin{pmatrix} 0 \\ \xi_u \end{pmatrix}, \qquad (4.7)$$

to enable us to write the change of state given by the impact law as

$$\Delta x = -\rho_1(\xi + \mathcal{O}(\rho)). \qquad (4.8)$$

Focusing on how to determine $t_C$ as a function of $x$, we start by introducing a smooth coordinate transformation $\tilde{x} = \tilde{x}(x)$. This is done by using $\psi(x)$ and $\tau(x)$ to define new coordinates to use instead of, say, $q_1$ and $u_1$. Let $\tilde{x}_1 = \psi(x)$, $\tilde{x}_2 = \tau(x)$ and let the remaining $\tilde{x}$ coordinates be equal to $q_2, \dots, q_n, u_2, \dots, u_m$. This transformation lets us write $H$ on a simple form,

$$H(x, t) = \tilde{x}_1 + (t - \tilde{x}_2)^2 (A_g/2 + \mathcal{O}(\tilde{x}, t)), \qquad (4.9)$$

where $A_g$ is the acceleration of the grazing impact. Now introduce another coordinate transformation $\chi = \chi(\tilde{x})$, by taking $\chi_1 = \sqrt{-\tilde{x}_1}$, and leaving all other coordinates unchanged. Using this set of coordinates we have

$$H(x, t) = -\chi_1^2 + (t - \chi_2)^2 (A_g/2 + \mathcal{O}(\chi, t)). \qquad (4.10)$$
By rewriting $H(x, t_C) = 0$ as

$$ \chi_1 - (\chi_2 - t_C)\sqrt{A_g/2 + \mathcal{O}(\chi, t_C)} = 0, \quad (4.11) $$

where the minus sign is chosen since it is the time to reach $\Sigma^-$ we are interested in, we obtain an equation which can be used to determine $t_C$ as a smooth function of $\chi$ to desired order. We solve to find

$$ t_C = \chi_2 - \chi_1(\sqrt{2/A_g} + \mathcal{O}(\chi)). \quad (4.12) $$

In order to express $\Delta x$ as a function of $\chi$ we need to calculate $v$ as a function of $\chi$,

$$ v = \frac{\partial H}{\partial t}(x, t_C) = -\chi_1(\sqrt{2A_g} + \mathcal{O}(\chi)). \quad (4.13) $$

This leads to

$$ \Delta x = \chi_1(\sqrt{2A_g}\xi + \mathcal{O}(\chi)). \quad (4.14) $$

There remains to expand the flow near $\mathcal{O}^*$, and we find that

$$ \Phi_t(x) = x + t(F(x^*) + \mathcal{O}(x, t)). \quad (4.15) $$

It is now easy to plug $t_C$ into $\Phi_t$, use the impact law and then $\Phi_t$ again. One obtains

$$ \Phi_{t_C}(x) = x + t_C(F(x^*) + \mathcal{O}(\chi)), \quad (4.16) $$

and by using (4.14)

$$ G \circ \Phi_{t_C}(x) = x + t_C(F(x^*) + \mathcal{O}(\chi)) + \chi_1(\sqrt{2A_g}\xi + \mathcal{O}(\chi)). \quad (4.17) $$

We finally obtain

$$ \Phi_{-t_C} \circ G \circ \Phi_{t_C}(x) = x + \chi_1(\sqrt{2A_g}\xi + \mathcal{O}(\chi)). \quad (4.18) $$

Transforming back to the original set of coordinates we see that for impacting points, $C$ can be written as

$$ C(x) = x + y(x)(\sqrt{2A_g}\xi + \mathcal{O}(x, y(x))), \quad (4.19) $$

where $y = \sqrt{-\psi(x)}$.

## 5. Stability and behaviour after bifurcation

In the following we will analyse the stability issue and what kind of motion one might expect after bifurcation. Although we use the assumption that the dynamics is given by the flow, except for the low velocity impacts, this is merely a question of convenience. It should be clear that impacts or any other discontinuity in forces could be allowed as long as they can be handled by differentiable local mappings.

### (a) The stability problem

When $\mu$ is increased to zero, $x^*$ will become a grazing periodic point, with some period time $T$: $\Phi_T(x^*) = x^*$. We let $L$ be the linearized part of $\Phi_T$,

$$ L = \frac{\partial \Phi_T}{\partial x}(x^*). \quad (5.1) $$

If we disregard impacts, $\Phi_T$ gives all the dynamics. Stability is assumed, hence $L$ has one eigenvalue equal to one and all other eigenvalues are inside the unit circle.
When the possibility of impacts is to be included, we know from the previous section that this can be done in terms of the discontinuity bypass mapping. For impacting points we have $\psi(x) \le 0$. One finds the expansion of $\psi$ to be $\psi = \eta x + O(x)^2$, where $\eta$ is the covector

$$ \eta = \frac{\partial h}{\partial x}(x^*). \qquad (5.2) $$

We note that $h$ only depends on generalized coordinates, so $\eta\xi = 0$.

It follows that the stability can be examined by the mapping

$$ x \mapsto \begin{cases} \sqrt{-2A_g\eta x}L\xi + O(x), & \text{if } \eta x \le 0, \\ Lx + O(x)^2, & \text{if } \eta x > 0. \end{cases} \qquad (5.3) $$

If $\eta L\xi < 0$, we see that an impact will be followed by another impact. If this is the case, then the square root will be iterated and as a result the periodic motion will become unstable. By just changing $L$ to $L^i$ above, this will happen if an impact is followed by non-impacting motion for some iterations, but eventually impacts again. Thus, if $\eta L^i\xi < 0$ for any $i \ge 1$, stability will be lost.

If we have equality for some $i$, $\eta L^i\xi = 0$, no conclusions can be drawn. This is because an understanding of this case must involve higher order terms. Lastly, if $\eta L^i\xi > 0$ for all $i \ge 1$, an impact is followed by non-impacting motion. The non-impacting motion was assumed to be stable, so in the generic case stability is assured. However, non-generic cases are possible where trailing terms will be of importance. We will identify these cases when discussing Poincaré mappings.

### (b) Poincaré mappings

A common method to employ when analysing periodic motion in systems described by differential equations is the Poincaré mapping technique. This will be most helpful in our case, because it lets us disregard the action of the mapping $C$ in the direction of the flow.

Thus, we assume that a suitable Poincaré section $\Sigma_P$ can be chosen, which is crossed transversally by the flow. It is described by coordinates $z$ and a mapping $V$, such that $V(z) = x \in \Sigma_P$. When $\mu = 0$ the grazing flowline will cross the section at $\bar{x}$, and we will assume that $z$ are local coordinates, hence $V(0) = \bar{x}$. The time of flight from $\bar{x}$ to $x^*$ is $t_1$, $\Phi_{t_1}(\bar{x}) = x^*$, and we let $t_2 = T - t_1$.

If impacts are disregarded, one notes that $\Phi_T \circ V$ will take a point close to $z = 0$ to a point $x$ close to $\bar{x}$, but not necessarily in the section. We define the mapping $W$ by starting at $x$, following the flowline for a short time $T_W(x)$ until the section is reached, and then let $z = W(x)$ be the coordinate for the corresponding crossing point. It is clear that $W \circ V$ becomes the identity mapping. We then have that the Poincaré mapping can be written $W \circ \Phi_T \circ V$. We include impacts by writing the complete Poincaré mapping as

$$ P = \hat{C} \circ \hat{P}, \qquad (5.4) $$

where $\hat{P} = W \circ \Phi_T \circ V$. In this way the effects of impacts are separated out from the smooth dynamics. Obviously, $\hat{C}$ must be related to $C$ and reduce to the identity mapping for non-impacting points. One finds that $\hat{C}$ can be written as $\hat{C} = W \circ \Phi_{t_2} \circ C \circ \Phi_{-t_2} \circ V$.

To be able to identify impacting points in $\Sigma_P$ we define $\psi' = \psi \circ \Phi_{t_1} \circ V$, so if $\psi'(z) < 0$ then $z$ is an impacting point. This gives a natural way to define a
covector $\eta'$,

$$ \eta' = \frac{\partial}{\partial z}(0) = \eta \frac{\partial \Phi_{t_1}}{\partial x}(\bar{x}) \frac{\partial V}{\partial z}(0), \qquad (5.5) $$

such that if second order terms are neglected, we have impact if $\eta'z \le 0$. The mapping $\hat{P}$ can be written as $\hat{P}(z) = L'z + \mathcal{O}(z)^2$, where

$$ L' = \frac{\partial W}{\partial x}(\bar{x}) \frac{\partial \Phi_T}{\partial x}(\bar{x}) \frac{\partial V}{\partial z}(0). \qquad (5.6) $$

Since $W$ is implicitly defined by

$$ \Phi_{-T_W(x)}(V(W(x))) = x, \qquad (5.7) $$

one can calculate $\partial W/\partial x(\bar{x})$ by differentiating this expression. By expanding $\hat{C}$ one finds a square root stretching in $\Sigma_P$. The direction is given by $\xi'$,

$$ \xi' = \frac{\partial W}{\partial x}(\bar{x}) \frac{\partial \Phi_{t_2}}{\partial x}(x^*) \xi. \qquad (5.8) $$

From now on we will discuss all dynamics in terms of the Poincaré mapping, so there will be no confusion if we drop all prime signs and simply write $\psi, \eta, L$ and $\xi$.

Assembling all pieces now gives the full Poincaré mapping in the form

$$ P(z) = \begin{cases} \sqrt{-2A_g\eta z}\xi + \mathcal{O}(z), & \text{if } \psi(z) \le 0, \\ Lz + \mathcal{O}(z)^2, & \text{if } \psi(z) > 0. \end{cases} \qquad (5.9) $$

By similar reasoning as in the previous section, the periodic motion is unstable if $\eta L^i \xi < 0$ for any non-negative integer $i$. As indicated earlier, when formulating a criterion for stability some care must be taken. To see this, assume that $L$ has a single largest (by modulus) eigenvalue and assume that the corresponding eigenvector is annihilated by $\eta$. It might be possible to have $\eta L^i \xi > 0$ for all $i \ge 0$, but the asymptotic motion will be in the plane given by $\eta$. Hence, a complete understanding requires higher order terms. There are similar situations, and they can all be viewed as special cases; arbitrary small changes in parameters will make the periodic motion unstable. Keeping away from this situation, a sufficient criterion for stability can be formulated: the grazing periodic point will be stable if

$$ \eta L^i \xi > 0 \text{ for } i = 0, 1, 2, \dots \qquad (5.10) $$

and this is valid also for neighbouring parameter values.

Figure 3 shows the geometric interpretation of the stability criterion. The successive iterates of $\xi$ always point in the positive direction out from the tangent plane to $= 0$.

### (c) Motion when $\mu > 0$

When $\mu > 0$ the mapping $\hat{P}$ will have a stable periodic point $Z(\mu)$, and by assumption we have $\psi(Z) < 0$, $\partial Z/\partial \mu \ne 0$. We introduce a new local coordinate $\tilde{z}$ by a $\mathcal{O}(\mu)$ transformation, $\tilde{z} = z - Z(\mu)$, so that the periodic point of $\hat{P}$ will correspond to $\tilde{z} = 0$. Writing the Poincaré mapping for non-impacting points as a function of $\tilde{z}$ we have $\hat{P} = L\tilde{z} + \mathcal{O}(\tilde{z}, \mu)^2$. If we let $\psi(Z) = -d$, we find by expressing $\psi$ as function of $\tilde{z}$ that $\psi = \eta\tilde{z} - d + \mathcal{O}(\tilde{z}, \mu)^2$. By truncating terms of order greater than one, the following mapping is obtained,

$$ \tilde{P}(\tilde{z}, \mu) = \begin{cases} \sqrt{2A_g(d - \eta\tilde{z})} \xi + K(\tilde{z}, \mu), & \text{if } \eta\tilde{z} \le d, \\ L\tilde{z}, & \text{if } \eta\tilde{z} > d, \end{cases} \qquad (5.11) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. The stability criterion.</figcaption></figure>

where $K$ is $O(\tilde{z}, \mu)$ and $K(\tilde{z}, \mu) = L\tilde{z}$ when $\eta\tilde{z} = d$. In the following we will assume that the stability criterion is fulfilled, and to simplify notation we write $z$ instead of $\tilde{z}$.

Denote the region where $\eta z \le d$ by I, and the region where $\eta z > d$ by II. We can see that $z = 0$ will not be a fixed point any more. This is because $\sqrt{2A_g d\eta\xi} + \eta K(0, \mu) > d$ for small $d$ values, so the image of $z = 0$ will be in region II. Since the eigenvalues of $L$ are inside the unit circle, the iterates of $z = 0$ will approach region I. Investigating how many iterations of the linear mapping that will be needed for $z = 0$ to return to region I, we see from

$$ \sqrt{2A_g} \eta L^i \xi + \frac{\eta L^i K(0, \mu)}{\sqrt{d}} \le \sqrt{d}, \quad (5.12) $$

that $i \to \infty$ as $d \to 0$. Assume that $L$ has a single largest eigenvalue (by modulus) $\lambda$, and let $\zeta$ be the corresponding eigenvector. From the stability criterion we have that $\lambda$ must be real and $\lambda > 0$. We choose $\zeta$ such that $\eta\zeta = 1$ and by denoting all other eigenvalues by $\lambda_r$ we find

$$ L^i \tilde{P}(0, \mu) = \sqrt{2A_g d} \alpha \lambda^i \zeta + \sqrt{d} O(\lambda_r)^i + O(\lambda, \lambda_r)^i O(\mu), \quad (5.13) $$

where $\alpha > 0$ is a constant. When approaching region I, we find that the iterates will line up on the line segment given by $\zeta$. We expect this behaviour also for points with $\eta z < 0$, at least for points close to $z = 0$ and small $d$ values. The crossing point of this line segment and the surface $\eta z = d$ is $z_\zeta = d\zeta$, which gives $\eta L z_\zeta = \lambda d$. Thus, in the limit $d \to 0$, this is the maximal penetration for an orbit which has spent a large number of iterations in region II. It follows that the motion will be trapped in a region of size $O(\sqrt{d}) = O(\sqrt{\mu})$.

Since all orbits entering region II finally will return to region I, it makes sense to study the first return mapping from I to I. Let $N(z, \mu)$ be the number of iterations for $z$ to re-enter region I. Clearly, from a stability point of view it is most advantageous to spend many iterations in region II. If the return to region I takes place in just a few iterations, the derivatives of the return mapping will be large, so no stable periodic motion is to be expected. For large $N$, the dominating part of the return mapping is $\alpha\sqrt{2A_g}\sqrt{d - \eta z}\lambda^{N(z,\mu)}\zeta$. Due to the one dimensionality for high iterates, it is sufficient to study the projection of $L^{N(z,\mu)}\tilde{P}(z, \mu)$, where $\eta z \le d$, using $\eta$. If $\delta$
is defined as $\delta(z) = \eta z/d$, we obtain

$$ \delta(z) \mapsto \delta(L^{N(z,\mu)} \tilde{P}(z, \mu)) \approx \alpha \sqrt{2A_g} (\lambda^{N(z,\mu)}/\sqrt{d}) \sqrt{1-\delta(z)}. \quad (5.14) $$

Nordmark (1996) shows by writing $d$ as a function of two new parameters, that the dependence on $z$ in equation (5.14) as $d \to 0$ is only through $\delta(z)$. The resulting one-dimensional limit mapping is then analysed. Results for similar one-dimensional mappings can be found in Nusse et al. (1994), Lamda & Budd (1994) and in Foale & Bishop (1994). An analysis of a two-dimensional mapping with the properties here described can be found in Chin et al. (1994). We briefly summarize the results from these investigations. If $\lambda > 2/3$, the motion is found to be chaotic for an interval of $\mu$ values, $0 < \mu < \mu'$. If $1/4 < \lambda < 2/3$ there will be an infinite sequence of windows with stable periodic motion alternating with chaotic bands. The period increases by 1 if $\mu$ is decreased from one periodic window to the next. There is a scaling relationship for the periodic windows: if period $N$ is found for $\mu$, then period $N+1$ can be found close to $\lambda^2\mu$. When $0 < \lambda < 1/4$ no chaotic bands are found. The scaling relationship still holds, but the windows will overlap.

## 6. A numerical example

A model of a two degrees of freedom (autonomous) oscillator will have a three-dimensional Poincaré section. With this in mind we construct an example mapping based on the theory presented above. We let

$$ z = \begin{pmatrix} z_1 \\ z_2 \\ z_3 \end{pmatrix} \qquad (6.1) $$

and the covector $\eta$ is taken to be

$$ \eta = (1 \ 0 \ 0). \qquad (6.2) $$

When $\eta z > d$, the dynamics is given by a linear mapping $L$, $z \mapsto Lz$. We choose $L$ to be on the simple form,

$$ L = \begin{pmatrix} \lambda & 0 & 0 \\ 0 & 0 & -\omega \\ 0 & \omega & 0 \end{pmatrix}. \qquad (6.3) $$

The eigenvalues of $L$ are $\lambda$ and $\pm i\omega$, and it is assumed that $0 < |\omega| < \lambda < 1$. For the impact part of the mapping, when $\eta z \le d$, we use

$$ z \mapsto \sqrt{d - \eta z} \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} + \begin{pmatrix} \lambda d + 3(d - z_1) \\ 5(z_1 - d) - \omega z_3 \\ 7(d - z_1) + \omega z_2 \end{pmatrix}. \qquad (6.4) $$

Thus, the vector which gives the direction of the square root stretching is

$$ \xi = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}. \qquad (6.5) $$
<figure><img src="image_6.png" /><figcaption>Figure 4. Bifurcation diagram, $\lambda = 0.5$, $\omega = 0.3$.</figcaption></figure>

It is easy to see that this choice of $\eta$, $L$ and $\xi$ fulfils the stability criterion. The bifurcation should then be continuous and the phenomena mentioned are expected to show up in simulations as $d$ is increased from zero.

In figure 4 a bifurcation diagram is shown, where $z_1$ is plotted against $d$. Using $\lambda = 0.5$ and $\omega = 0.3$, the example mapping has been simulated for 600 different $d$ values. For each value of $d$, the last 1000 iterates are shown from a total of 10000 iterates. When $d = 0.03$, the mapping has a stable four-periodic solution. As $d$ is decreased, the four-periodic solution is interrupted by a chaotic band, which eventually disappears, leaving a five-periodic orbit. This seems to repeat as $d$ is decreased further. To investigate this more thoroughly and to check the scaling relationship, another bifurcation diagram was plotted. It is displayed in figure 5. By dividing $z_1$ with $\sqrt{d}$ the range of the motion should be constant and by plotting against $\ln \sqrt{d} / \ln \lambda$ the width of the chaotic bands should also be constant. This is clearly seen to be the case.

Choosing $d = 1.5 \times 10^{-4}$ we find that $\ln \sqrt{d} / \ln \lambda \approx 6.35$, so this corresponds to chaotic motion. This value can be used to show the appearance of the attractor. Figure 6 shows 2000 iterates from a simulation including 10000 iterations of the mapping. The iterates line up in a many fingered object. This is similar to the result obtained from a two-dimensional mapping in Nordmark (1991, p. 294, figure 9). The longest finger consists of points that were located on the impacting side in the previous step. The expansion given by the square root can be seen, the finger is less dense closer to zero. Since most points end up well inside the non-impacting region, the linear mapping comes into action, giving the other fingers. Finally, impact occurs again and the pattern is repeated in a non-periodic way.

## 7. Results and discussion

In this paper the bifurcation of a stable periodic orbit in a general mechanical oscillator as it grazes a rigid boundary has been considered. By using the assumption
<figure><img src="image_7.png" /><figcaption>Figure 5. Bifurcation diagram using logarithmic scale.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. 2000 iterations when $d = 1.5 \times 10^{-4}$.</figcaption></figure>

that impacts can be modelled by an instantaneous change of velocities, given by a smooth impact law, we have been able to analyse this situation to some extent. To understand the interaction between the flow and the impact surface close to a grazing impact, it has been helpful to use the concept of a discontinuity bypass mapping. By series expansion we have found that for points that undergo a low velocity impact, this mapping has a square root dependence of the penetration depth if trajectories were allowed to move inside the impact surface. The square root term acts only in
one direction, which is given by the impact law. Combining this with the dynamics of the non-impacting motion gives a criterion for stability.

The bifurcation will be discontinuous if the square root term of the mapping makes the grazing periodic orbit unstable. The square root term will be iterated, hence the motion quickly moves away from the formerly stable orbit. The situation then escapes local analysis, but it should be noted that the local theory gives the direction in which the system at first moves. With a knowledge of the dynamics in a larger region of state space, it might be possible to understand this case as well. The direction of the square root stretching is given by the impact law, so a change in the impact law might, for certain systems, make different regions of the state space reachable. This could be useful in applications, where a certain kind of motion usually is preferable.

If the stability criterion is fulfilled, then the bifurcation will be continuous. The motion after bifurcation is trapped close to the former fixed point. The dynamics is characterized by occasional low velocity impacts that are followed by non-impacting motion for a while. Depending on the largest eigenvalue of the linearized non-impacting part of the Poincaré mapping, the motion will be periodic, with arbitrary long periods, or non-periodic. This is due to the one-dimensional nature of the dynamics.

Impacts are possible causes of wear and noise, an important issue in many applications. The theory presented above gives tools to handle and understand some of the problems involved in bifurcations when a part of a mechanism starts to impact. The main underlying cause is the tangency of the flow to a boundary in state space, a situation that could occur in a variety of different systems. Analysing situations when higher order terms are of importance and the dependence of parameters in the stability criterion are topics for future work.

We thank Dr Hanno Essén for helpful comments. The support of this work by the Swedish Research Council for Engineering Sciences (TFR) is gratefully acknowledged.

## References

Budd, C. & Dux, F. 1994 Intermittency in impact oscillators close to resonance. *Nonlinearity* **7**, 1191–1224.

Chin, W., Ott, E., Nusse, H. E. & Grebogi, C. 1994 Grazing bifurcations in impact oscillators. *Phys. Rev. E* **50**, 4427–4444.

Foale, S. & Bishop, S. R. 1992 Dynamical complexities of forced impacting systems. *Phil. Trans. R. Soc. Lond.* **A 338** 547–556.

Foale, S. & Bishop, S. R. 1994 Bifurcations in impact oscillators. *Nonlinear Dynamics* **6**, 285–299.

Ivanov, A. P. 1993 Stabilization of an impact oscillator near grazing incidence owing to resonance. *J. Sound Vib.* **162**, 562–565.

Ivanov, A. P. 1994 Impact oscillations: linear theory of stability and bifurcations. *J. Sound Vib.* **178**, 361–378.

Kane, T. R. & Levinson, D. A. 1985 *Dynamics: theory and applications*. New York: McGraw Hill Publishing.

Knudsen, C., Feldberg, R. & True, H. 1992 Bifurcations and chaos in a model of a rolling railway wheelset. *Phil. Trans. R. Soc. Lond.* **A 338**, 455–469.

Lamba, H. & Budd, C. 1994 Scaling of lyapunov exponents at nonsmooth bifurcations. *Phys. Rev. E* **50**, 84–90.

Lesser, M. 1995 *The analysis of complex nonlinear mechanical systems*. Singapore: World Scientific Publishing.
Nordmark, A. B. 1991 Non-periodic motion caused by grazing incidence in an impact oscillator. *J. Sound Vib.* **145**, 279–297.

Nordmark, A. B. 1997 A universal limit mapping in grazing bifurcations. *Phys. Rev. E*. (In the press.)

Nusse, H. E., Ott, E. & Yorke, J. A. 1994 Border-collision bifurcations: a possible explanation for observed bifurcation phenomena. *Phys. Rev. E* **49**, 1073–1076.

Shaw, S. W. & Holmes, P. J. 1983 A periodically forced piecewise linear oscillator. *J. Sound Vib.* **90**, 129–155.

Whiston, G. S. 1992 Singularities in vibro-impact dynamics. *J. Sound Vib.* **152**, 427–460.