# On normal form calculations in impact oscillators

Mats Fredriksson, Arne Nordmark

► To cite this version:

Mats Fredriksson, Arne Nordmark. On normal form calculations in impact oscillators. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1999, 456 (1994), pp.315-329. 10.1098/rspa.2000.0519. hal-01510821

HAL Id: hal-01510821

https://hal.science/hal-01510821

Submitted on 20 Apr 2017

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On normal form calculations in impact oscillators

MATS H. FREDRIKSSON AND ARNE B. NORDMARK

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

Normal form calculations are useful for analysing the dynamics close to bifurcations. However, the application to non-smooth systems is a topic for current research. Here we consider a class of impact oscillators, where we allow systems with several degrees of freedom as well as nonlinear equations of motion. Impact is due to the motion of one body, constrained by a motion limiter. The velocities of the system are assumed to change instantaneously at impact. By defining a *discontinuity mapping*, we show how Poincaré mappings can be obtained as an expansion in a local coordinate. This gives the mapping the desired form, thus making it possible to employ standard techniques. All calculations are algorithmic in spirit, hence computer algebra routines can easily be developed.

Keywords: impact oscillations; normal forms; non-smooth systems; discontinuity mapping; computer algebra

## 1. Introduction

Modern dynamical systems theory supplies tools for local analysis of the dynamics of ordinary differential equations. Concepts such as Poincaré mappings, centre manifolds and normal forms are helpful in understanding motion and bifurcations. The underlying assumption for these techniques is that the vector field is sufficiently well behaved. However, in models of mechanical systems one can find that assumptions that are natural from the modelling point of view give vector fields that lack smoothness across boundaries in state space. Oscillating systems where impacts are possible is one example where it is not straightforward to employ standard methods.

In impact oscillators one or more motion limiters are present. Typically, a wall acts as a one-sided amplitude constraint. This situation is not uncommon in engineering systems, where models that incorporate this behaviour have been used to understand gears (Pfeiffer & Kunert 1990; Karagiannis & Pfeiffer 1991), railway wheelsets (Nordström Jensen & True 1997; Knudsen et al. 1992), mooring towers (Thompson 1983), printing devices (Tung & Shaw 1988), the mechanics of bipedal walking (Garcia et al. 1998; McGeer 1993; Adolfsson et al. 1999), and several other dynamical problems.

The motion during contact has a much shorter time-scale than the typical motion of the system; consequently, a very common model for the impact is to assume that velocities change at an instant. This is the model that will be used in this paper.
Many papers have dealt with this subject. Among the first investigations using a dynamical systems perspective we find Shaw & Holmes (1983), Shaw (1985a, b) and Thompson & Ghaffari (1983). Subsequent work has revealed dynamical features found in smooth systems, as well as special features of the model, such as the dynamics close to grazing impact (see Nordmark 1991), and the possibility of an infinite number of impacts in finite time (see Budd & Dux 1994). The main focus has been on one-degree-of-freedom systems, but systems with several degrees of freedom have been considered in Shaw & Shaw (1989), Aidanpää & Gupta (1993), Fredriksson (1997) and the thesis of Lee (1996). The study of such systems are, in general, complicated by the abundance of parameters, thus special cases are usually considered in order to reduce the number of parameters. To the authors' knowledge no review paper has been written, but the theme issue edited by Bishop (1994) is dedicated to impact oscillations.

Previous investigations usually assume that the dynamics between impacts is given by linear equations of motion. As the impact law relates the velocity after impact to the velocity before impact, a coefficient of restitution is used. The advantage of this approach is that a closed-form solution for the motion between impacts can conveniently be written down, at least as long as the system has one degree of freedom. For impact oscillators little attention has been paid to more general approaches to local stability calculations, where one cannot rely on special features of the system. A scheme for the derivation of local expressions would be desirable, giving means to analyse systems where the interesting dynamics is impossible to capture using linearized equations of motion. This would also be of interest for the engineering community. The aim of the present paper is to propose such an idea. Combined with standard techniques it gives a way to obtain a local Poincaré mapping $P$ on the form

$$P(z) = Lz + N(z), \quad (1.1)$$

where $z$ is a local coordinate in the Poincaré section, $L$ is the matrix that gives the linear part of the mapping and $N(z)$ are the higher-order terms. To obtain (1.1) is the starting point if one would like to use normal forms to study bifurcations.

The paper is organized as follows. Firstly, we discuss the type of systems that we have in mind. Local mappings are introduced, and an idea on how to handle the impact part of the local dynamics is presented. This is done by defining a discontinuity mapping. This mapping encapsulates the contribution to the dynamics coming from impacts. Three examples are studied. Finally, we discuss the results and how these can be extended.

## 2. Motion and impact

Our assumptions about the system are largely inspired by models based on connected rigid bodies. Such multibody models often offer sufficiently accurate and computationally efficient models. They also relate well to the impact approximation, since a truly rigid body must change velocities discontinuously upon impact with a rigid surface. We assume that impact is due to a single body colliding with a surface. Thus, only one condition for contact needs to be checked. The motion of the impacting body is coupled to a larger system (in an arbitrary way) between impacts.
(a) State space and flow

The system is assumed to have $n$ degrees of freedom. We denote coordinates on the configuration manifold by $q_1, q_2, \dots, q_n$. Similarly, the velocities of the system are $u_1, u_2, \dots, u_n$. The state $x$ is

$$x = \begin{bmatrix} q_1 \\ u_1 \\ \vdots \\ q_n \\ u_n \end{bmatrix} \qquad (2.1)$$

and the time derivative $\dot{x}$ of a motion is related to $x$ by a vector field $F$:

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

In many practical cases, time-periodic forcing is present. Indeed, all examples below have periodic forcing. In developing the general arguments, we do not need to treat this case separately. We can always rewrite the equations in the standard fashion, by extending $x$ with a phase angle $\theta \in S^1$, and add $\dot{\theta} = 1$ to the equations of motion.

We denote the state space where the motion limiter is removed (or thought to be penetrable) by $\mathcal{X}$. This is partly a notational convenience. The main concern here is how to obtain locally valid expressions, thus we will not have much to say about global features of the state space. We write $\Phi = \Phi(x, t)$ for the state space flow:

$$\Phi : \mathcal{X} \times \mathcal{R} \to \mathcal{X}, \qquad (2.3)$$

the mapping given by the solution of the differential equation (2.2). When writing the arguments explicitly we sometimes use an indexed time $t$, $\Phi_t(x)$. Thus, $\Phi_t(x)$ is the state reached by following a trajectory from an initial point $x$ during the time interval $t$.

(b) The impact law

We define the function $H = H(x)$ to be equal to the distance from the impacting body to the motion limiting surface when $H(x) > 0$. The equation $H(x) = 0$ is the condition for contact. We also allow $H(x) < 0$, meaning a state which breaks the geometry imposed by the model. The contact condition can be interpreted as a surface in $\mathcal{X}$ with an equation $H(x) = 0$. We denote the set of points that fulfils this equation as $\Sigma$:

$$\Sigma = \{x \in \mathcal{X} : H(x) = 0\}. \qquad (2.4)$$

When a trajectory reaches $\Sigma$, it is disrupted by the impact law. We view the impact law as a mapping $G : \mathcal{X} \to \mathcal{X}$, which is only of interest on $\Sigma$. The impact law leaves coordinates on the configuration manifold unchanged. The velocities are changed, where we allow the new velocities to depend on both velocities and the configuration coordinates of the impact. We write

$$x_a = G(x_b), \qquad (2.5)$$

where $x_a$ is the state immediately after impact, and $x_b$ is the state immediately before impact.
<figure><img src="image_2.png" /><figcaption>Figure 1. Impacting and non-impacting motion.</figcaption></figure>

### 3. The discontinuity mapping

We now develop the main idea of the paper. Related work can be found in Fredriksson & Nordmark (1997) and in Dankowicz & Nordmark (2000). For a derivation of the linearization, see also Müller (1995).

#### (a) Motivation

Away from $\Sigma$, all of the dynamics is in the flow. If we wish to analyse non-impacting motion close to a point $x'$, it is natural to consider the mapping $\Phi_T$, where $T > 0$ is the time-interval of interest (the period in the case of periodic motion). A trivial decomposition can be made,

$$ \Phi_T = \Phi_{t_2} \circ \Phi_{t_1}, \qquad (3.1) $$

where $t_1 + t_2 = T$.

Including impacts, we assume that $t_1$ is also the time of flight for an orbit starting at a point $\bar{x}$ to reach a point $\tilde{x} \in \Sigma$:

$$ \tilde{x} = \Phi_{t_1}(\bar{x}). \qquad (3.2) $$

Using the impact law, and then $\Phi_{t_2}$, the image of $\bar{x}$ is

$$ \Phi_{t_2} \circ G \circ \Phi_{t_1}(\bar{x}). \qquad (3.3) $$

This situation is indicated in figure 1.

However, one should note that the expression $\Phi_{t_2} \circ G \circ \Phi_{t_1}$ is valid as a mapping only for points having a flight time equal to $t_1$. In general the flight time for points in a neighbourhood of $\bar{x}$ will be different. A decomposition similar to that of flow maps would be convenient, which inspires the following idea: can we find a mapping $C$ such
<figure><img src="image_3.png" /><figcaption>Figure 2. The discontinuity mapping. Flowlines are dashed when $H < 0$.</figcaption></figure>

that the correct mapping for all points close to $\bar{x}$ can be written as $\Phi_{t_2} \circ C \circ \Phi_{t_1}$? This mapping $C$ must then incorporate the jump in velocities implied by the mapping $G$ as well as the differing flight times. The mapping $C$ will be referred to as the *discontinuity mapping*, since it takes care of all aspects of passing the discontinuity imposed by the impact law.

Using this viewpoint, the influence of impacts is in a natural way separated from the dynamics of the flow. Calculating the flow for a fixed time is also natural in numerical computations, where $\Phi$ can be locally expanded and the coefficients can be determined by integrating the variational equations for the fixed times $t_1$ and $t_2$.

If, firstly, we follow the flow for a fixed time $t_1$, the image of the starting point might have $H < 0$. This is not consistent with the model, and the discontinuity mapping $C$ will have to address this. The definition set for $C$ is thus a neighbourhood $N$ of $\tilde{x}$, where $N \subset X$. The first step in defining $C$ is to note that for $x \in N$ we can trace the flowline passing through $x$ to the intersection with $\Sigma$. We use $t_C = t_C(x)$ to denote the short time of flight from $x$ to $\Sigma$. If $H(x) < 0$ then $t_C(x) < 0$, and if $H(x) > 0$ the time of flight to $\Sigma$ is positive. The latter situation is shown in figure 2.

In this way we obtain a mapping from $N$ to $\Sigma$, and the image of $x$ is $\Phi_{t_C(x)}(x)$. Using the impact law we get a point $G \circ \Phi_{t_C(x)}(x)$ close to $G(\tilde{x})$. We now use this as an initial condition and we follow the motion (disregarding the wall) for a time $-t_C(x)$. The discontinuity mapping is then

$$C(x) = \Phi_{-t_C(x)}(G \circ \Phi_{t_C(x)}(x)). \qquad (3.4)$$

It follows that $C$ is a mapping from a state space neighbourhood of $\tilde{x}$ (disregarding the wall) to a state space neighbourhood of $G(\tilde{x})$ (again disregarding the wall). Choosing $x = \tilde{x}$ we have $t_C(\tilde{x}) = 0$; hence $C(\tilde{x}) = G(\tilde{x})$, as expected.
(b) Expansions

To write C in a useful form, we wish to express it as an expansion in $y = x - \tilde{x}$:

$$C(x) = G(\tilde{x}) + Ly + N(y), \qquad (3.5)$$

where L is the matrix giving the linearized mapping and N denotes nonlinear terms. It is not hard to translate the ideas in the previous section into a step-by-step description of how to obtain this expression.

(i) Firstly we calculate the flow close to $\tilde{x}$. We use the tilde to stress the expansion point, thus writing $\tilde{\Phi}$. It is sought as

$$\tilde{\Phi} = x + t\tilde{\Gamma}(y, t), \qquad (3.6)$$

where $\tilde{\Gamma}$ is an expansion in $y$ and $t$. We determine the unknown coefficients in $\tilde{\Gamma}(y, t)$ by using the fact that the flow satisfies

$$\frac{\partial \tilde{\Phi}}{\partial t}(x, t) = F(\tilde{\Phi}(x, t)) \qquad (3.7)$$

for all $x, t \in \mathcal{X} \times \mathbf{R}$.

(ii) The next step is to find the time $t_C(x)$ as an expansion in $y$. The unknown coefficients in $t_C$ are calculated by inserting the expression for the flow into an expansion of $H$:

$$H(\tilde{\Phi}(x, t_C)) = 0. \qquad (3.8)$$

(iii) We can now calculate the point of impact, $x_b$:

$$x_b = \tilde{\Phi}(x, t_C(x)). \qquad (3.9)$$

(iv) We now direct our attention towards the impact law $G$. Using the notation $\tilde{G}_x$ for the Jacobian of $G$ at $\tilde{x}$, we write the expansion as

$$G(x) = G(\tilde{x}) + \tilde{G}_x y + \mathcal{O}(y)^2, \qquad (3.10)$$

where $\mathcal{O}(y)^2$ denotes terms of order two and higher. We then have

$$x_a = G(x_b), \qquad (3.11)$$

where $x_b(x)$ is known.

(v) The flow close to $\hat{x} = G(\tilde{x})$ is now calculated. This is similar to the first step, with only a change in the expansion point. For obvious reasons it is convenient to denote the initial condition with a subscript ‘a’, and we use the circumflex as we have previously used the tilde:

$$\hat{\Phi} = x_a + t\hat{\Gamma}(x_a - \hat{x}, t). \qquad (3.12)$$

(vi) Lastly, we insert $t = -t_C$ in the expression for the flow:

$$C(x) = \hat{\Phi}(x_a, -t_C), \qquad (3.13)$$

where $x_a = x_a(x)$ has been calculated in previous steps.

Finding the linearization is straightforward. The zeroth-order term in $\tilde{\Gamma}$ is $\tilde{F} = F(\tilde{x})$, thus

$$\tilde{\Phi} = x + t(\tilde{F} + \mathcal{O}(y, t)). \qquad (3.14)$$
By writing

$$ \tilde{H}_x = \frac{\partial H}{\partial x}(\tilde{x}), \qquad (3.15) $$

the function $H(x)$ is

$$ H(x) = \tilde{H}_x y + \mathcal{O}(y)^2, \qquad (3.16) $$

where the gradient $\tilde{H}_x$ is a row vector. We use an unknown row vector $\tilde{t}_{C_x}$ to tentatively write $t_C$ as

$$ t_C = \tilde{t}_{C_x} y + \mathcal{O}(y)^2. \qquad (3.17) $$

By inserting this into (3.8) we get

$$ \tilde{H}_x y + \tilde{H}_x \tilde{F} \tilde{t}_{C_x} y + \mathcal{O}(y)^2 = 0, \qquad (3.18) $$

thus

$$ \tilde{t}_{C_x} = - \frac{\tilde{H}_x}{\tilde{H}_x \tilde{F}}. \qquad (3.19) $$

From (3.9) and (3.14), we obtain

$$ x_b = \tilde{x} + \left( I - \frac{\tilde{F} \tilde{H}_x}{\tilde{H}_x \tilde{F}} \right) y + \mathcal{O}(y)^2, \qquad (3.20) $$

where $I$ is the identity matrix. Using the impact law we find

$$ x_a = G(\tilde{x}) + \tilde{G}_x \left( I - \frac{\tilde{F} \tilde{H}_x}{\tilde{H}_x \tilde{F}} \right) y + \mathcal{O}(y)^2, \qquad (3.21) $$

and by substituting $x_a$ and $t = -t_C$ into $\tilde{\Phi}$ we have

$$ C(x) = G(\tilde{x}) + \left\{ \tilde{G}_x \left( I - \frac{\tilde{F} \tilde{H}_x}{\tilde{H}_x \tilde{F}} \right) + \frac{\tilde{F} \tilde{H}_x}{\tilde{H}_x \tilde{F}} \right\} y + \mathcal{O}(y)^2, \qquad (3.22) $$

thus

$$ L = \tilde{G}_x + \frac{(\tilde{F} - \tilde{G}_x \tilde{F}) \tilde{H}_x}{\tilde{H}_x \tilde{F}}. \qquad (3.23) $$

Note that this form of the linearized map is uniquely specified by demanding that vectors $v$ orthogonal to $\tilde{H}_x$ should be mapped to $\tilde{G}_x v$ (since $C$ is the same as $G$ on $\Sigma$) and that $\tilde{F}$ should map to $\tilde{F}$ (since $C$ maps a flowline through $\tilde{x}$ to a flowline through $\hat{x}$).

If we wish to derive the mapping to higher order the computations rapidly get complex. A computer algebra system is very helpful in order to automate the calculations and to minimize errors. For error checking it is also helpful to note some of the characteristics of the mapping. If we tentatively use the identity mapping for $G$, then $C$ is identity. If we take $x \in \Sigma$, then $C = G$. These observations can be used to check that complicated expressions evaluate as expected.
## 4. Examples

Let us investigate how the ideas above can be used in practical calculations. When writing out expansions we use brackets `[i, j, ..., k]` as indices to label coefficients. If $A$ is a function of the real variables $x_1, x_2, x_3$, we write

$$A = \sum_{\substack{i,j,k \ge 0 \\ i+j+k<m}} \frac{\tilde{A}_{[i,j,k]}}{i!j!k!} (x_1 - \tilde{x}_1)^i (x_2 - \tilde{x}_2)^j (x_3 - \tilde{x}_3)^k + \mathcal{O}(x - \tilde{x})^m. \quad (4.1)$$

An expansion that includes linear terms is then

$$A = \tilde{A}_{[0,0,0]} + \tilde{A}_{[1,0,0]}(x_1 - \tilde{x}_1) + \tilde{A}_{[0,1,0]}(x_2 - \tilde{x}_2) + \tilde{A}_{[0,0,1]}(x_3 - \tilde{x}_3) + \mathcal{O}(x - \tilde{x})^2. \quad (4.2)$$

### (a) A forced n-degree-of-freedom system

We introduce some assumptions to obtain a system with an arbitrary number of degrees of freedom, but with similarities to the familiar one-degree-of-freedom periodically forced case. The basis for the one-degree-of-freedom model is a particle moving along a line. This motion is now coupled to another system. Thus, the state space without motion limiter, $\mathcal{X}$, has the structure

$$\mathcal{X} = \mathbf{R}^2 \times \mathcal{M} \times \mathbf{S}^1, \quad (4.3)$$

where the factor $\mathbf{R}^2$ is due to the position and velocity of the particle, $\mathcal{M}$ denotes the submanifold of the state space which describes the system that is coupled to the particle and the factor $\mathbf{S}^1$ is included by assuming periodic forcing. Furthermore, we assume that the equations of motion are

$$\begin{bmatrix} \dot{q}_1 \\ \dot{u}_1 \\ \vdots \\ \dot{q}_n \\ \dot{u}_n \\ \dot{\theta} \end{bmatrix} = \begin{bmatrix} u_1 \\ A_1(x) \\ \vdots \\ u_n \\ A_n(x) \\ 1 \end{bmatrix}. \quad (4.4)$$

When the motion limiter is taken into the model we assume that the impact process only involves the particle model. The coordinate $q_1$ describes the position of the particle and the motion limiter is located at a critical coordinate $q_1 = q_c$, hence

$$H(x) = q_1 - q_c. \quad (4.5)$$

The impact law only changes the velocity of the particle, which we write as

$$u_{1a} = g(u_{1b}). \quad (4.6)$$

All other velocities are unaffected by the impact. Writing the impact law as

$$u_{1a} = \hat{u}_1 + g_{[1]}(u_{1b} - \tilde{u}_1) + \mathcal{O}(u_{1b} - \tilde{u}_1)^2, \quad (4.7)$$

and denoting

$$A_{i[0,...,0]} = A_{i[0]}, \quad (4.8)$$
we obtain from (3.23) the expression for the linearization L of C :

$$
L =
\begin{bmatrix}
\frac{\hat{u}_1}{\tilde{u}_1} & 0 & 0 & 0 & \cdots & 0 & 0 & 0 \\
\frac{\hat{A}_{1_{[0]}} - g_{[1]} \tilde{A}_{1_{[0]}}}{\tilde{u}_1} & g_{[1]} & 0 & 0 & \cdots & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & \cdots & 0 & 0 & 0 \\
\frac{\hat{A}_{2_{[0]}} - \tilde{A}_{2_{[0]}}}{\tilde{u}_1} & 0 & 0 & 1 & \cdots & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & \cdots & 1 & 0 & 0 \\
\frac{\hat{A}_{n_{[0]}} - \tilde{A}_{n_{[0]}}}{\tilde{u}_1} & 0 & 0 & 0 & \cdots & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1
\end{bmatrix}.
\quad (4.9)
$$

The local state-space volume change caused by an impact is given by the determinant of L . Since L is triangular we immediately have

$$
\det L = g_{[1]} \frac{\hat{u}_1}{\tilde{u}_1}. \qquad (4.10)
$$

This is independent of the dimension of the system. For the much-used model where the impact law is modelled by a coefficient of restitution r ,

$$
u_{1a} = -ru_{1b}, \tag{4.11}
$$

we have

$$
\det L = r^2. \tag{4.12}
$$

(i) Poincaré mappings

For periodically forced systems without impacts, the canonical candidate for a Poincaré section is Σ ¯ θ :

$$
\Sigma_{\bar{\theta}} = \{x \in \mathcal{X} : \theta = \bar{\theta}\}, \quad (4.13)
$$

for some choice of ¯ θ. The Poincaré mapping using this section is sometimes called the stroboscopic mapping and denoted as P$_{S}$ . For impacting systems the situation is more complicated. The impact velocity is often of major interest. A natural idea is then to use a subset of Σ ( u$_{1}$ > 0 or u$_{1}$ < 0) as a Poincaré section. The mapping is then referred to as the impact mapping, P$_{I}$ . Another choice is to use Σ ¯ θ (with q$_{1}$ > q$_{c}$ ) and use this set as a section. This choice is more in conformity with the current approach. Assume that we have a motion, making m impacts at the phases ˜ θ$_{1}$, ˜ θ$_{2}$, . . . , ˜ θ$_{m}$ before returning to Σ ¯ θ ; let P ˜ θ$_{i}$$_{→}$ ˜ θ$_{1}$ be the mapping, using only the flow, from a neighbourhood of the starting point in Σ ¯ θ to the section of fixed phase where the first impact occurs. Similarly we write P ˜ θ$_{1}$$_{→}$ ˜ θ$_{2}$, P ˜ θ$_{2}$$_{→}$ ˜ θ$_{3}$, . . . , P ˜ θ$_{m}$$_{→}$ ˜ θ . We use C | θ = ˜ θ to denote the restriction of C by fixing the phase angle to the impact phase. The stroboscopic mapping is then

$$
P_S = P_{\tilde{\theta}_m \to \tilde{\theta}} \circ C|_{\theta=\tilde{\theta}_m} \circ \dots \circ P_{\tilde{\theta}_1 \to \tilde{\theta}_2} \circ C|_{\theta=\tilde{\theta}_1} \circ P_{\tilde{\theta} \to \tilde{\theta}_1}. \quad (4.14)
$$
To find periodic motion, standard root-finding algorithms can be employed. The linearization of $C|_{\theta=\tilde{\theta}}$ is then useful. We obtain it simply by removing the last row and last column in (4.9). It is worth noting that the determinant is the same for the restricted mapping, thus the Poincaré section volume change due to impact is also given by (4.10).

(b) A general one-degree-of-freedom system

We now consider a one-degree-of-freedom system, and we calculate an expansion of $C$ to order two. We use (4.4) with $n = 1$, where for simplicity we write $q = q_1$, $u = u_1$ and $A = A_1$. The special structure of the system can be used to find shortcuts in the calculations. Using $Q$, $U$ and $\Theta$ for the components of the flow, we make the trial solution

$$Q = q + t(u + \frac{1}{2}t\alpha(x, t)), \qquad (4.15)$$

where we use $\tilde{\alpha}$ for expansions at $(\tilde{x}, 0)$, and similarly we write $\hat{\alpha}$ to stress that $\hat{x}$ is the point of expansion. By differentiating $Q$ with respect to time we obtain $U$:

$$U = \frac{\partial Q}{\partial t} = u + t \left( \alpha(x, t) + \frac{1}{2} t \frac{\partial \alpha(x, t)}{\partial t} \right). \qquad (4.16)$$

The unknown coefficients in $\alpha$ can be determined by differentiating $U$ with respect to time and inserting all expressions in an expansion of $A$. We obtain

$$\alpha + 2t \frac{\partial \alpha}{\partial t} + \frac{1}{2} t^2 \frac{\partial^2 \alpha}{\partial t^2} = A(Q, U, \Theta). \qquad (4.17)$$

Carrying out the first step, we wish to obtain $U$ to order two, hence we need $\tilde{\alpha}$ to order one:

$$\tilde{\alpha} = \tilde{\alpha}_{[0,0,0,0]} + \tilde{\alpha}_{[1,0,0,0]}(q - \tilde{q}) + \tilde{\alpha}_{[0,1,0,0]}(u - \tilde{u}) \\ + \tilde{\alpha}_{[0,0,1,0]}(\theta - \tilde{\theta}) + \tilde{\alpha}_{[0,0,0,1]}t + \mathcal{O}(x - \tilde{x}, t)^2. \qquad (4.18)$$

Expanding $A$ to order one and using (4.17), we get

$$\left.
\begin{aligned}
\tilde{\alpha}_{[0,0,0,0]} &= \tilde{A}_{[0,0,0]}, \\
\tilde{\alpha}_{[1,0,0,0]} &= \tilde{A}_{[1,0,0]}, \\
\tilde{\alpha}_{[0,1,0,0]} &= \tilde{A}_{[0,1,0]}, \\
\tilde{\alpha}_{[0,0,1,0]} &= \tilde{A}_{[0,0,1]}, \\
\tilde{\alpha}_{[0,0,0,1]} &= \frac{1}{3}(\tilde{A}_{[1,0,0]}\tilde{u} + \tilde{A}_{[0,0,0]}\tilde{A}_{[0,1,0]} + \tilde{A}_{[0,0,1]}).
\end{aligned}
\right\} \qquad (4.19)$$

To find $t_C$ we note that since $t_C = 0$ when $q = \tilde{q}$, it is natural to seek $t_C$ of the form

$$t_C(x) = (q - \tilde{q})\beta(x), \qquad (4.20)$$

where $\beta$ is an expansion in $x-\tilde{x}$. The unknown coefficients in $\beta$ are found by inserting $t = t_C$ into $Q(x, t) = \tilde{q}$, from which we obtain

$$1 + (u + (q - \tilde{q})\frac{1}{2}\beta\alpha)\beta = 0. \qquad (4.21)$$
This gives the first-order expansion of $\beta$:

$$ 
\left. 
\begin{aligned}
\tilde{\beta}_{[0,0,0]} &= -\frac{1}{\tilde{u}}, \\
\tilde{\beta}_{[1,0,0]} &= -\frac{1}{2\tilde{u}^3}\tilde{A}_{[0,0,0]}, \\
\tilde{\beta}_{[0,1,0]} &= \frac{1}{\tilde{u}^2}, \\
\tilde{\beta}_{[0,0,1]} &= 0.
\end{aligned}
\right\} \qquad (4.22)
$$

Inserting $t_C$ into $U$ gives $u_b$:

$$ 
\begin{aligned}
u_b = & u - \frac{1}{\tilde{u}}\tilde{A}_{[0,0,0]}(q - \tilde{q}) \\
& + \frac{1}{2\tilde{u}}\left\{-\tilde{A}_{[1,0,0]} + \frac{1}{\tilde{u}}\tilde{A}_{[0,0,1]} + \frac{1}{\tilde{u}}\tilde{A}_{[0,0,0]}\tilde{A}_{[0,1,0]} - \frac{1}{\tilde{u}^2}\tilde{A}_{[0,0,0]}^2\right\}(q - \tilde{q})^2 \\
& + \frac{1}{\tilde{u}}\left\{-\tilde{A}_{[0,1,0]} + \frac{1}{\tilde{u}}\tilde{A}_{[0,0,0]}\right\}(q - \tilde{q})(u - \tilde{u}) \\
& - \frac{1}{\tilde{u}}\tilde{A}_{[0,0,1]}(q - \tilde{q})(\theta - \tilde{\theta}) + \mathcal{O}(x - \tilde{x})^3.
\end{aligned}
\qquad (4.23)
$$

The expansion of the impact law is

$$ u_a = \hat{u} + g_{[1]}(u_b - \tilde{u}) + \frac{1}{2}g_{[2]}(u_b - \tilde{u})^2 + \mathcal{O}(u_b - \tilde{u})^3. \qquad (4.24) $$

To get the expansion of the flow close to $\hat{x}$ we just change the tilde to a circumflex. Inserting $t = -t_C$ and the initial conditions, we find that the position component $q_C$ of $C(x)$ is

$$ 
\begin{aligned}
q_C = \tilde{q} &+ \frac{\hat{u}}{\tilde{u}}(q - \tilde{q}) + \frac{1}{2\tilde{u}^2}\left\{\hat{A}_{[0,0,0]} + \left(\frac{\hat{u}}{\tilde{u}} - 2g_{[1]}\right)\tilde{A}_{[0,0,0]}\right\}(q - \tilde{q})^2 \\
& + \frac{1}{\tilde{u}}\left\{g_{[1]} - \frac{\hat{u}}{\tilde{u}}\right\}(q - \tilde{q})(u - \tilde{u}) + \mathcal{O}(x - \tilde{x})^3,
\end{aligned}
\qquad (4.25)
$$

and the velocity component $u_C$ is

$$ 
\begin{aligned}
u_C = & \hat{u} + \frac{1}{\tilde{u}}\{\hat{A}_{[0,0,0]} - g_{[1]}\tilde{A}_{[0,0,0]}\}(q - \tilde{q}) + g_{[1]}(u - \tilde{u}) \\
& + \left\{\frac{1}{2\tilde{u}}\left(\frac{\hat{u}}{\tilde{u}}\hat{A}_{[1,0,0]} - g_{[1]}\tilde{A}_{[1,0,0]}\right) - \frac{1}{2\tilde{u}^2}(\hat{A}_{[0,0,1]} - g_{[1]}\tilde{A}_{[0,0,1]})\right. \\
& + \frac{1}{2\tilde{u}^3}\tilde{A}_{[0,0,0]}(\hat{A}_{[0,0,0]} - g_{[1]}\tilde{A}_{[0,0,0]}) \\
& + \frac{1}{2\tilde{u}^2}(\hat{A}_{[0,0,0]}\hat{A}_{[0,1,0]} - 2g_{[1]}\hat{A}_{[0,1,0]}\tilde{A}_{[0,0,0]} + g_{[1]}\tilde{A}_{[0,0,0]}\tilde{A}_{[0,1,0]}) \\
& + \left.\frac{g_{[2]}}{2\tilde{u}^2}\tilde{A}_{[0,0,0]}^2\right\}(q - \tilde{q})^2 \\
& + \left\{\frac{g_{[1]}}{\tilde{u}}(\hat{A}_{[0,1,0]} - \tilde{A}_{[0,1,0]}) - \frac{1}{\tilde{u}^2}(\hat{A}_{[0,0,0]} - g_{[1]}\tilde{A}_{[0,0,0]}) - \frac{g_{[2]}}{\tilde{u}}\tilde{A}_{[0,0,0]}\right\}(q - \tilde{q})(u - \tilde{u}) \\
& + \frac{1}{\tilde{u}}\{\hat{A}_{[0,0,1]} - g_{[1]}\tilde{A}_{[0,0,1]}\}(q - \tilde{q})(\theta - \tilde{\theta}) + \frac{1}{2}g_{[2]}(u - \tilde{u})^2 + \mathcal{O}(x - \tilde{x})^3.
\end{aligned}
\qquad (4.26)
$$

The phase-angle component $\theta_C$ is simply

$$ \theta_C = \theta. \qquad (4.27) $$
(c) A saddle-node bifurcation

In the case where the equations of motion between impacts are linear, all nonlinearities are in $C$. Assuming that we study an orbit with one impact, we can use the section with $\theta = \tilde{\theta}$ to write the Poincaré mapping as

$$P_S(z) = P_{\tilde{\theta} \to \tilde{\theta}} \circ C|_{\theta=\tilde{\theta}}(z), \quad (4.28)$$

where $z$ is the Poincaré section coordinate and $P_{\tilde{\theta} \to \tilde{\theta}}$ is a linear mapping. This viewpoint can be practical both for (numerically) finding periodic orbits and for studying the effects of the nonlinear terms. As an example we take the one-degree-of-freedom system with the acceleration given by

$$A = -\omega^2 q - 2\delta\omega u + a\{(\omega^2 - 1)\cos(\theta) - 2\delta\omega \sin(\theta)\}, \quad (4.29)$$

and with the impact law

$$u_a = -ru_b. \quad (4.30)$$

Here the forcing is chosen so that if impacts are disregarded the motion will settle on $q = a \cos(\theta)$, $u = -a \sin(\theta)$. In the following $a$ is varied while the other parameters are kept fixed. Since we work with the restriction $C|_{\theta=\tilde{\theta}}$, the point of expansion is the triplet $\tilde{q}$, $\tilde{u}$, $a_c$. We use the notation $C_{[i,j;k]}$ for the coefficients in the expansion of $C|_{\theta=\tilde{\theta}}$. Hence $C_{[1,1;1]}$ is the coefficient in front of $(q - \tilde{q})(u - \tilde{u})(a - a_c)$.

We can obtain the coefficients $C_{[i,j;k]}$ by applying the result of the previous section. The Jacobian of (the restriction of) $C$ is $L_C$,

$$L_C = \begin{bmatrix} -r & 0 \\ \frac{(1+r)\gamma}{\tilde{u}} & -r \end{bmatrix}; \quad (4.31)$$

the $C_{[0,0;1]}$ coefficient is zero,

$$C_{[0,0;1]} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}; \quad (4.32)$$

and the coefficients for the quadratic terms are

$$C_{[2,0;0]} = \begin{bmatrix} \frac{(1+r)\gamma}{\tilde{u}^2} \\ \frac{(1+r)}{\tilde{u}^3} \end{bmatrix}, \quad (4.33)$$

$$C_{[1,1;0]} = \begin{bmatrix} 0 \\ -\frac{(1+r)\gamma}{\tilde{u}^2} \end{bmatrix}, \quad (4.34)$$

$$C_{[1,0;1]} = \begin{bmatrix} 0 \\ \frac{(1+r)\{(\omega^2 - 1)\cos(\tilde{\theta}) - 2\delta\omega \sin(\tilde{\theta})\}}{\tilde{u}} \end{bmatrix}, \quad (4.35)$$

$$C_{[0,2;0]} = C_{[0,0;2]} = C_{[0,1;1]} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \quad (4.36)$$

where

$$\gamma = a_c\{(\omega^2 - 1)\cos(\tilde{\theta}) - 2\delta\omega \sin(\tilde{\theta})\} - \omega^2 \tilde{q} \quad (4.37)$$
and

$$ = \gamma(\gamma - 4\delta\omega\tilde{u}) + a_c\tilde{u}\{(\omega^2 - 1)\sin(\tilde{\theta}) + 2\delta\omega\cos(\tilde{\theta})\}. \quad (4.38) $$

The mapping $P_{\tilde{\theta} \to \tilde{\theta}}$ is

$$ P_{\tilde{\theta} \to \tilde{\theta}}(z; a) = P_{\tilde{\theta} \to \tilde{\theta}}(\hat{z}; a_c) + L_T(z - \hat{z}) + (I - L_T) \begin{bmatrix} \cos(\tilde{\theta}) \\ -\sin(\tilde{\theta}) \end{bmatrix} (a - a_c), \quad (4.39) $$

where

$$ L_T = e^{BT}, \quad (4.40) $$

$$ B = \begin{bmatrix} 0 & 1 \\ -\omega^2 & -2\delta\omega \end{bmatrix}, \quad (4.41) $$

$I$ is the identity matrix and $T$ is the period of interest.

Choosing $\omega = 1/4.9$, $\delta = 0.6$, $r = 1$ and $q_c = -1$, we can find a saddle-node bifurcation involving 4-periodic motion with one impact when $a = a_c = 0.91566$. The impact velocity is $\tilde{u} = -0.36815$ and the impact phase is $\tilde{\theta} = 2.72795$. To approximately find the bifurcated motion we calculate the centre manifold and the reduced dynamics. Using $\chi$ as a coordinate, we write an expansion for the centre manifold as

$$ z = \tilde{z} + \phi_1 \chi + M_{[0,1]}(a - a_c) + M_{[2,0]} \frac{1}{2} \chi^2 + \dots, \quad (4.42) $$

where $\phi_1$ is the eigenvector corresponding to the eigenvalue $\lambda_1 = 1$ of $L = L_T L_C$. There is an element of choice involved in normal form calculations. We have here taken

$$ \phi_1^* M_{[0,1]} = 0, \quad \phi_1^* M_{[2,0]} = 0, \quad (4.43) $$

where $\phi_1^*$ is the left eigenvector of $L$ with eigenvalue $\lambda_1$, to make the expansion unique. A comparison of the approximate solution with a solution found numerically can be done in two ways. This is because a numerical routine for finding periodic orbits typically assumes that $q > q_c$. Thus, either we find the bifurcated solution numerically at a convenient phase and move it to the section $\tilde{\theta}$, or else we move the approximate solution to the section used for the numerical solution.

The latter approach is taken in figure 3, where a bifurcation diagram at $\theta = 0$ is shown. The eventual fate of the branches is thoroughly investigated in Nordmark (2000). Both the stable and the unstable branch are destroyed in grazing bifurcations as $a$ is increased (the stable when $a = 0.97555$, the unstable when $a = 1$).

# 5. Results and discussion

By using the discontinuity mapping, we have shown how to obtain series expansions of mappings in impact oscillators with a unilateral displacement limiter. This is helpful in analysing bifurcations. Some extensions of this work are straightforward. We can change the assumption that the impact law can be written explicitly, by allowing an implicit relation $G(x_a, x_b)$ as long as it can be solved locally. Systems where only one impact can occur at a given time, but with many different contact possibilities, can be brought into the present framework by applying the different discontinuity mappings
<figure><img src="image_4.png" /><figcaption>Figure 3. Bifurcation diagram when $\theta = 0$. The unstable branch is dashed.</figcaption></figure>

several times, once for each contact. Other aspects are more subtle and call for future studies. Models where spatially separated parts of one or several bodies can make contact simultaneously have not received much attention. In such models two contact conditions can be fulfilled at the same time. One scenario that might be possible is that a parameter change will change the order in which the impact events occur for a periodic motion. Since linear stability is given by multiplying the Jacobians of the different mappings, which, when composed, gives the Poincaré mapping; this corresponds to switching order between the Jacobians of two discontinuity mappings. Since matrix multiplication in general does not commute, this might lead to non-differentiable dynamics exactly when the order is switched.

We thank Dr Hanno Essén for comments. This work was supported by the Swedish Research Council for Engineering Sciences (TFR).

## References

Adolfsson, J., Dankowicz, H. & Nordmark, A. B. 2000 3D passive walkers: finding periodic gaits in the presence of discontinuities. (Submitted.)

Aidanpää, J. O. & Gupta, R. B. 1993 Periodic and chaotic behaviour of a threshold-limited two-degree-of-freedom system. *J. Sound Vib.* **165**, 305–327.

Bishop, S. R. (ed.) 1994 *Impact oscillators* (Theme issue). *Phil. Trans. R. Soc. Lond. A 347*.

Budd, C. & Dux, F. 1994 Chattering and related behaviour in impact oscillators. *Phil. Trans. R. Soc. Lond. A 347*, 365–389.

Dankowicz, H. & Nordmark, A. B. 2000 On the origin and bifurcations of stick-slip oscillations. *Physica D* **136**, 280–302.
Fredriksson, M. H. 1997 Grazing bifurcations in multibody systems. *Nonlinear Analysis Theor. Meth. Applications* **30**, 4475–4483.

Fredriksson, M. H. & Nordmark, A. B. 1997 Bifurcations caused by grazing incidence in many degrees of freedom impact oscillators. *Proc. R. Soc. Lond. A* **453**, 1261–1276.

Garcia, M., Chatterjee, A., Ruina, A. & Coleman, M. 1998 The simplest walking model: stability, complexity and scaling. *ASME Jl Biomech. Engng* **120**, 281–288.

Karagiannis, K. & Pfeiffer, F. 1991 Theoretical and experimental investigations of gear-rattling. *Nonlinear Dynam.* **2**, 367–387.

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.

Lee, G. 1996 The theoretical and numerical analysis of impact oscillators. PhD thesis, University of Bristol.

McGeer, T. 1993 Dynamics and control of bipedal locomotion. *J. Theor. Biol.* **163**, 277–314.

Müller, P. C. 1995 Calculation of Lyapunov exponents for dynamic systems with discontinuities. *Chaos. Solitons Fractals* **5**, 1671–1681.

Nordmark, A. B. 1991 Non-periodic motion caused by grazing incidence in an impact oscillator. *J. Sound Vib.* **145**, 279–297.

Nordmark, A. B. 2000 Existence of periodic orbits in grazing bifurcations of impacting mechanical oscillators. (In preparation.)

Nordström Jensen, C. & True, H. 1997 On a new route to chaos in railway dynamics. *Nonlinear Dynam.* **13**, 117–129.

Pfeiffer, F. & Kunert, A. 1990 Rattling models from deterministic to stochastic processes. *Nonlinear Dynam.* **1**, 63–74.

Shaw, J. & Shaw, S. W. 1989 The onset of chaos in a two-degree-of-freedom impacting system. *J. Appl. Mech.* **56**, 168–174.

Shaw, S. W. 1985a Dynamics of harmonically excited systems having rigid amplitude constraints. I. Subharmonic motions and local bifurcations. *J. Appl. Mech.* **52**, 453–458.

Shaw, S. W. 1985b Dynamics of harmonically excited systems having rigid amplitude constraints. II. Chaotic motions and global bifurcations. *J. Appl. Mech.* **52**, 459–464.

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

Thompson, J. M. T. 1983 Complex dynamics of compliant off-shore structures. *Proc. R. Soc. Lond. A* **387**, 407–427.

Thompson, J. M. T. & Ghaffari, R. 1983 Chaotic dynamics of an impact oscillator. *Phys. Rev. A* **27**, 1741–1743.

Tung, P. C. & Shaw, S. W. 1988 The dynamics of an impact print hammer. *J. Vib. Acoust. Stress Reliability Design* **110**, 193–200.