# The use of normal forms for analysing nonlinear mechanical vibrations

Simon Neild, Alan Champneys, David Wagg, Thomas Hill, Andrea Cammarano

► To cite this version:

Simon Neild, Alan Champneys, David Wagg, Thomas Hill, Andrea Cammarano. The use of normal forms for analysing nonlinear mechanical vibrations. Philosophical Transactions of the Royal Society A: Physical and Engineering Sciences (1990–1995), 2015, 373 (2051), pp.20140404. 10.1098/rsta.2014.0404. hal-01308768

HAL Id: hal-01308768

https://hal.science/hal-01308768

Submitted on 5 May 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
PHILOSOPHICAL
TRANSACTIONS A

rsta.royalsocietypublishing.org

Research

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

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

CrossMark

click for updates

**Cite this article:** Neild SA, Champneys AR, Wagg DJ, Hill TL, Cammarano A. 2015 The use of normal forms for analysing nonlinear mechanical vibrations. *Phil. Trans. R. Soc. A* **373**: 20140404.

http://dx.doi.org/10.1098/rsta.2014.0404

Accepted: 23 June 2015

One contribution of 11 to a theme issue 'A field guide to nonlinearity in structural dynamics'.

**Subject Areas:**

structural engineering, mechanical
engineering

**Keywords:**

nonlinear dynamics, normal forms,
modal analysis, cable vibration

**Author for correspondence:**

Simon A. Neild

e-mail: simon.neild@bristol.ac.uk

The use of normal forms for
analysing nonlinear
mechanical vibrations

Simon A. Neild<fnref n="1" />, Alan R. Champneys<fnref n="1" />, David J. Wagg<fnref n="2" />,
Thomas L. Hill<fnref n="1" /> and Andrea Cammarano<fnref n="1" />,†

A historical introduction is given of the theory of normal forms for simplifying nonlinear dynamical systems close to resonances or bifurcation points. The specific focus is on mechanical vibration problems, described by finite degree-of-freedom second-order-in-time differential equations. A recent variant of the normal form method, that respects the specific structure of such models, is recalled. It is shown how this method can be placed within the context of the general theory of normal forms provided the damping and forcing terms are treated as unfolding parameters. The approach is contrasted to the alternative theory of nonlinear normal modes (NNMs) which is argued to be problematic in the presence of damping. The efficacy of the normal form method is illustrated on a model of the vibration of a taut cable, which is geometrically nonlinear. It is shown how the method is able to accurately predict NNM shapes and their bifurcations.

# 1. Introduction: nonlinear normal modes and normal forms

The problem addressed here is how to extend the well-established notion of normal modes of linear vibration systems to nonlinear systems in a mathematically consistent way that also allows for practical implementation. In recent years, there has been a lot of research related to the concept of a so-called ‘nonlinear normal mode’ (NNM). We shall not go into the complete theory here, but refer instead to Kerschen *et al.* [1] and Peeters *et al.* [2]) and Avrimov & Mikhlin [3] for surveys of the

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

<footnote n="1">¹Faculty of Engineering, University of Bristol, Bristol BS8 1TR, UK</footnote>

<footnote n="2">²Department of Mechanical Engineering, University of Sheffield, Sheffield S1 3JD, UK</footnote>

<footnote>†Present address: School of Engineering, University of Glasgow, Glasgow, G12 8QQ, UK.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Illustrating the limitation of the NNM concept in the presence of damping. Computations are for the unforced two-degree-of-freedom system whose radial parts are written in polar coordinates as $r_1 = -\delta_1 r_1 - r_1^3 + (3 + 0.5 \cos \theta_2)r_2^3$, $r_2 = -\delta_2 r_2 - r_2^3 + (0.2 \sin(\theta_1) - 1)r_1^3$, where (a) $\delta_1 = \delta_2 = 0$, (b) $\delta_1 = 0.1$, $\delta_2 = 0.05$ and (c) $\delta_1 = 0.05$, $\delta_2 = 0.1$. We show projections, ignoring the angular motion, of the invariant manifolds of the system that are tangent to the linear eigenspaces $r_1 = 0$ and $r_2 = 0$.</figcaption></figure>

state of the art. As we shall explain, this concept is of limited use when forcing and damping are present, a restriction that is not shared by the more general concept of the normal form of a nonlinear system near an equilibrium [4–6]. At present though, the application of normal form analysis to structural vibration systems [7,8], described in more detail in §3 below, is rather technical and not widely adopted. It is also not clear how these structural vibration normal forms relate to the generic theory of normal forms as treated in Murdock [6]. Thus, the aim of this paper is to explain the normal form theory for structural vibrations in a more rigorous and historical context, and to show its relevance for a variety of practical problems.

Consider mechanical structures: the dynamics can be modelled as

$$M\ddot{x} + D\dot{x} + Kx + N(x, \dot{x}, r) = Pr + c.c., \quad (1.1)$$

where $x \in \mathbb{R}^n$ represents the displacements of the $n$ degrees of freedom of the structure, $M, D$ and $K$ are, respectively, mass, damping and stiffness matrices and $N$ represents nonlinearity, which is assumed to be sufficiently smooth. The right-hand side represents the forcing on the system $r(t) \in \mathbb{C}^n$ and can be written as a sum of $e^{i\Omega_1 t}$ terms. The concept of the NNM applies well to the case where there is no forcing or damping, so that $P = D = 0$ in (1.1) and $N$ is a pure function of $x$. In this case, Shaw & Pierre [9,10] argued that an NNM is just an invariant manifold composed of periodic solutions, whose frequency in the limit of amplitude tending to zero is the same as the linear mode. Here, they can appeal to the Liapunov Centre Theorem that guarantees that such a manifold must exist. Attempts have been made to extend this definition to include damping. Here, things are more troublesome because the Liapunov centre theorem no longer applies and, in general, one should not expect there to be small periodic solutions. Nevertheless, one can still appeal to the theory of invariant manifolds, but unfortunately such manifolds will be non-unique in general and will depend crucially on the relative size of the damping in each mode; see figure 1 illustrating this point for a general two degrees of freedom system. Figure 1a illustrates NNMs the case of zero damping. These manifolds which are shown as projections onto a radial coordinate are composed entirely of periodic orbits, as given by the Liapunov centre theorem. By contrast, figure 1b,c shows damped cases where the damping coefficient is greater in either the first or the second mode. In each case, only the most strongly damped mode has a unique invariant manifold tangent to it, with solutions in the manifold now composed of decaying oscillations. For the more weakly damped mode, there are now infinitely many such invariant manifolds and there is no uniquely distinguished one. In general, a tiny perturbation in initial condition will lead to a very different manifold if one computes backwards in time.
In the presence of general forcing, there is no guarantee that these invariant manifolds persist. Instead, what one needs is a calculation tool; a way of simplifying the system that respects normal modal analysis in the limit of small amplitude response, but for moderate amplitude also allows for the inclusion of only the necessary nonlinear effects. This is precisely what the theory of normal forms does. A concise description of that theory forms the subject of §2.

The link between normal form theory and NNMs was pointed out by Montaldi et al. [11], who showed that NNMs can be constructed directly, as a consequence of the Birkhoff normal form (§2) for *undamped* systems written in Hamiltonian coordinates. More general normal analysis was carried out on various structural systems in the late 1970s and early 1980s by Holmes and co-workers, see Guckenheimer & Holmes [4] and references therein. Jezequel & Lamarque [7] showed how forcing and damping can be incorporated in the approximation of normal forms for systems equivalent to (1.1), but written in first-order, i.e. state-space, form. Their method is essentially a hybrid between normal form reduction and the harmonic balance method. Direct harmonic balance and other asymptotic methods, such as multiple scales analysis have been applied by many authors to derive approximations to systems of the form (1.1) [12].

Recently, Neild & Wagg [8] have extended the work of Jezequel & Lamarque [7] to deal directly with second-order systems of the form (1.1); see §1 for the details. One purpose of this paper is to survey their method and to demonstrate its power through applications. Primarily though, we shall show how the second-order method can be recast using the general theory of normal forms, where forcing and damping terms can be treated as specific forms of unfoldings. One does not need to make the additional approximations inherent in the harmonic balance method when deriving the normal form itself. Application of the harmonic balance method to the derived normal form, though, is shown to lead to remarkably powerful predictions for resonant responses, leading to so-called *backbone curves* in structural analysis, for mode switching in complex multimodal responses and for analytically finding NNMs.

The rest of this paper is outlined as follows. In §2, we describe the historical context of normal form transformations. Section 3 goes on to apply these concepts to engineering vibration problems of the form (1.1). In particular, we show how the Neild & Wagg method can be thought of as the choice of a particular type of normal form. We also explain how to use this method to predict resonant responses, by interpreting the normal form in the context of harmonic balance. Section 4 demonstrates the theory through application to a particular model of a taut cable. The power of the method is illustrated by showing how it can predict secondary bifurcations corresponding to switches in mode shapes. Finally, §5 draws conclusions and points to avenues of future work.

## 2. Historical background to normal form analysis

Much of our modern-day understanding of dynamical systems can be traced to the work of the French genius Henri Poincaré and the geometric theory he introduced [13] to understand systems like the three-body problem in celestial mechanics; see, for example, Barrow-Green [14], Diacu & Holmes [15], Verhulst [16] for historical reviews. His key insight was to work on approximating the system rather than producing series solutions for individual trajectories.

### (a) Birkhoff normal form

The original context for normal forms was that of conservative systems written in Hamiltonian form. Consider a Hamiltonian system in $\mathbb{R}^{2n}$ close to an elliptic equilibrium point, written in canonical form

$$ \begin{pmatrix} \dot{p} \\ \dot{q} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} \nabla H(p, q), \quad \text{where } H = \frac{1}{2} \sum_{j=1}^{n} \omega_j (q_j^2 + p_j^2) + H_1(q, p), \qquad (2.1) $$

in which $q \in \mathbb{R}^n$ represents the column vector of positions $(q_1, q_2, \dots, q_n)^T$, $p$ represents the momenta $(p_1, p_2, \dots, p_n)^T$ and $H_1$ contains all higher-power terms, starting with cubic. For the
purposes of this paper, we shall assume that all vector fields are analytic (infinitely smooth) and we shall suppress any parameter dependence for the time being. Without the nonlinear terms, the system (2.1) is easily solvable and decomposes into quasi-periodic motion with independent frequencies $\omega_j$, $j=1,2,\dots,n$. That is, in each degree-of-freedom, the motion comprises periodic orbits of period $2\pi/\omega_j$. We say that the system is completely integrable.

Would not it be nice, at least in a sufficiently small neighbourhood of the zero equilibrium, to find coordinate transformations to successively remove all the nonlinear terms from the equation (cubic and higher terms in $H$) so that the equation is once again trivially integrable up to any given power? Formally, such reductions can be performed using *canonical* (i.e. Hamiltonian-structure-preserving) *near-identity* (i.e. differing only in nonlinear terms) transformations; see, for example, Meyer *et al.* [17] for details.

But a naive approach can go wrong. Suppose, for example, $H_1$ contained a term $\gamma_j(p_j^2 + q_j^2)^2$ for some $j$ and some given non-zero coefficient $\gamma_i$, leading to terms $4\gamma_j p_j(p_j + q_j)^2$ in the $\dot{q}_j$ equation and $4\gamma_j q_j(p_j + q_j)^2$ in the $\dot{p}_j$ equation. Such a term is called *resonant* because, no matter what near-identity transformation is chosen, a nonlinear component of this term remains. The resonance occurs in terms that contain powers $p_j^{\alpha_j}$ and $q_j^{\beta_j}$, where certain integer relationships existing between the corresponding eigenvalues $\pm i\omega_j$, $j=1,\dots,n$. The specific fourth-order expression $\gamma_j(p_j^2 + q_j^2)^2$ contains resonant terms where $\alpha_j(i\omega_j)^2 + \beta_j(-i\omega_j)^2 = 0$ with $\alpha_j = \beta_j = 2$. Other resonances can arise due to internal resonances where two or more of the frequencies $\omega_j$ are rationally related (which leads to the so-called *small divisors* problem in Hamiltonian systems that had greatly exercised Poincaré and his contemporaries). A good choice of transformed variables are the so-called action-angle coordinates, and the resulting system is often called the *Birkhoff normal form* [18].

## (b) General first-order normal forms

The problem is that many vibration problems (1.1) cannot be written in Hamiltonian form owing to the presence of non-conservative forces due to excitation, damping, friction or gyroscopic effects. Even if a Hamiltonian formation does exist, the mechanical interpretation in terms of velocities and displacements is likely to be lost. Instead, there is a more general theory of normal forms that can be applied to arbitrary-dimensional dynamical systems of the form

$$ \dot{x} = Ax + N_x(x), \quad x \in \mathbb{R}^p \qquad (2.2) $$

in the neighbourhood of an equilibrium point, where $N_x$ represents nonlinear terms (with the subscript reminding us that we are in the original unscaled coordinates $x$). Here, we shall be interested in the case that the Jacobian is non-hyperbolic (that is, there are eigenvalues on the imaginary axis). Often this theory is implemented in conjunction with centre manifold theory so that only the nonlinear terms associated with the non-hyperbolic degrees of freedom are retained [4,19].

Note that the unforced version of (1.1) can be written in this form by setting the number of states $p=2n$ and

$$ A = \begin{bmatrix} 0 & I \\ -M^{-1}K & -M^{-1}D \end{bmatrix} \quad \text{and} \quad N_x = \begin{bmatrix} 0 \\ -M^{-1}N(x) \end{bmatrix}, \qquad (2.3) $$

in which $x = (x, \dot{x})^T$, as pointed out by Jezequel & Lamarque [7].

Now, let us suppose that a linear transformation $x = \mathcal{F}q$ is applied to (2.2) so that it is written in the form

$$ \dot{q} = \Lambda q + N_q(q), \quad q \in \mathbb{C}^p, \qquad (2.4) $$

where $N_q(q) = \mathcal{F}^{-1}N_x(\mathcal{F}q)$ and $\Lambda$ contains the linearization in the simplest (Jordan canonical) form.
We seek to systematically change variables by a sequence of near-identity transformations to remove non-resonant terms of successively higher powers in the Taylor series expansion of the nonlinear term $N_q$. That is, we seek a new variable $\mathbf{u} \in \mathbb{C}^p$ with

$$ \dot{\mathbf{q}} = \Lambda \mathbf{q} + N_q(\mathbf{q}) \xrightarrow{\mathbf{q} = \mathbf{u} + \mathbf{h}(\mathbf{u})} \dot{\mathbf{u}} = \Lambda \mathbf{u} + N_u(\mathbf{u}), \quad (2.5) $$

where $\mathbf{h} = o(\mathbf{u})$, and $N_u$ contains only resonant terms (in a sense to be made precise shortly). To find the unknown function $\mathbf{h}(\mathbf{u})$, we differentiate the transform, $\dot{\mathbf{q}} = (I + D\mathbf{h}(\mathbf{u}))\dot{\mathbf{u}}$, and substitute for $\dot{\mathbf{q}}$ and $\dot{\mathbf{u}}$ to give

$$ \Lambda \mathbf{u} + \Lambda \mathbf{h}(\mathbf{u}) + N_q(\mathbf{u} + \mathbf{h}(\mathbf{u})) = (I + D\mathbf{h}(\mathbf{u}))(\Lambda \mathbf{u} + N_u(\mathbf{u})). \quad (2.6) $$

On its own, (2.6) is a complicated, nonlinear functional equation for $\mathbf{h}(\mathbf{u})$. But we simplify it by expanding $\mathbf{h}(\mathbf{u})$ as a Taylor series and solve, where possible, for the unknown coefficients of $\mathbf{h}$, term by term. Thus, let $\mathbf{h}(\mathbf{u}) = \mathbf{h}^{(2)}(\mathbf{u}) + \mathbf{h}^{(3)}(\mathbf{u}) + \mathbf{h}^{(4)}(\mathbf{u}) + \dots$, where each $\mathbf{h}^{(k)}$ is a sum of homogeneous monomial terms of degree $k$. That is, the $i$th component of $\mathbf{h}^{(k)}$ takes the form

$$ h_i^{(k)} = \sum_{\mathbf{m}_k} h_{i, \mathbf{m}_k}^{(k)} \prod_{j=1}^{p} u_j^{m_j}, \quad (2.7) $$

where we use a vector notation for multi-indices such that $\mathbf{m}_k = (m_1, m_2, m_3, \dots, m_p)$, where $m_j$ is a whole number in the range $0 \le m_j \le k$ with the additional condition that $\sum_{j=1}^p m_j = k$ and the sum in (2.7) is over the set of all such indices. So, for example, for $p=2$ and $k=2$, we have $\mathbf{m}_k = (2, 0)$, (1, 1) or (0, 2) such that

$$ h_i^{(2)} = h_{i,2,0}^{(2)} u_1^2 + h_{i,1,1}^{(2)} u_1 u_2 + h_{i,0,2}^{(2)} u_2^2. $$

Now we are in a position to simplify (2.6). Suppose that coordinate transformations have taken place to remove all non-resonant terms up to $k-1$. Then, at $O(k)$, (2.6) reads

$$ \Lambda \mathbf{u} + \Lambda \mathbf{h}^{(k)}(\mathbf{u}) + N_q^{(k)}(\mathbf{u}) = \Lambda \mathbf{u} + D \mathbf{h}^{(k)}(\mathbf{u}) \Lambda \mathbf{u} + N_u^{(k)}(\mathbf{u}) + \text{h.o.t.}, $$

where again a superscript (k) on $N_q$ and $N_u$ indicates terms that are homogeneous polynomials of degree $k$ and the terms involve (k + 1)-powers and higher. This expression can be rearranged to read

$$ N_{u}^{(k)}(\mathbf{u}) - N_{q}^{(k)}(\mathbf{u}) = \Lambda \mathbf{h}^{(k)}(\mathbf{u}) - D \mathbf{h}^{(k)}(\mathbf{u}) \Lambda \mathbf{u} := [\Lambda, \mathbf{h}^{(k)}](\mathbf{u}), \quad (2.8) $$

which is known as the homological equation associated with the linear operator $\Lambda$. The operator $[\cdot, \cdot]$ on the right-hand side is known as the Lie bracket between $\Lambda$ and $\mathbf{h}$ (which is closely related to the so-called Poisson bracket in the case that the system is Hamiltonian).

It might be worthwhile to remind ourselves what we are trying to do. We want to choose the coefficients of all the terms in $\mathbf{h}^{(k)}$ to make $N_u^{(k)}$ as simple as possible. Note that the Lie bracket can be thought of as a linear operator acting on the set of homogeneous polynomials of degree $k$. If this linear operator is invertible, then the homological equation will have a unique solution for any choice of $N_u^{(k)}$. In particular, we are free to choose $N_u^{(k)} \equiv 0$, and then we will have found a unique solution for the coefficients of $\mathbf{h}^{(k)}$ to make all the $k$th power terms disappear in the simplified system. Terms that cannot be removed in $N_u^{(k)}$, the resonant terms, arise precisely when the Lie bracket is non-invertible.

Note that the Lie bracket as an operator depends on the matrix $\Lambda$, hence this matrix must contain all information necessary to define which terms, in the space of homogeneous polynomials of degree $k$, are resonant. To spell out which these terms are, it is most convenient to look at the homological equation (2.8) term by term to end up with an *indicial* form of the equation in terms of the indices $m_j$.
The indicial equation can be derived in complete generality for any Jordan canonical form of the matrix $\Lambda$ [6]. But for ease of explanation, we shall treat the simplified case where the eigenvalues of $\Lambda$ are semi-simple. That is, the geometric multiplicity of each eigenvalue is equal to its algebraic multiplicity, and hence $\Lambda$ is diagonalizable so that

$$ \Lambda = \operatorname{diag}\{\lambda_1, \lambda_2, \dots, \lambda_p\}, \qquad (2.9) $$

where $\lambda_i \in \mathbb{C}$ are the eigenvalues of $\Lambda$ allowing for multiplicity. Note how the unforced general vibration system (2.3) can be transformed using eigenvectors into (2.4) with $\Lambda$ of the form (2.9) with $(\lambda_{2i-1}, \lambda_{2i}) = (+\omega_i, -\omega_i)$, for $i = 1, \dots, n$.

Under the assumption (2.9), the indicial form of the homological equation reads

$$ \begin{aligned} N_{u,i,\mathbf{m}_k}^{(k)} \prod_{j=1}^{p} u_j^{m_j} - N_{q,i,\mathbf{m}_k}^{(k)} \prod_{j=1}^{p} u_j^{m_j} &= \lambda_i h_{i,\mathbf{m}_k}^{(k)} \prod_{j=1}^{p} u_j^{m_j} - h_{i,\mathbf{m}_k}^{(k)} \sum_{l=1}^{p} \frac{\partial}{\partial u_l} \left( \prod_{j=1}^{p} u_j^{m_j} \right) \lambda_l u_l \\ &= h_{i,\mathbf{m}_k}^{(k)} \left( \lambda_i - \sum_{l=1}^{p} m_l \lambda_l \right) \prod_{j=1}^{p} u_j^{m_j}, \end{aligned} \qquad (2.10) $$

where a subscript $(i, \mathbf{m}_k)$ on $N_q$ and $N_{it}$ represents the corresponding coefficient of the $k$th power term with vector index $\mathbf{m}_k$.

Looking at the form of equation (2.10), we find that we are free to choose $N_{u,i,\mathbf{m}_k}^{(k)} = 0$ and

$$ h_{i, \mathbf{m}_k}^{(k)} = \frac{N_{q, i, \mathbf{m}_k}^{(k)}}{(\lambda_i - \sum_{l=1}^{p} m_l \lambda_l)}, \qquad (2.11) $$

thus removing the $\prod_j u_j^{m_j}$ term from the simplified equation (the normal form), unless

$$ \lambda_i - \sum_{l=1}^{p} m_l \lambda_l = 0. \qquad (2.12) $$

The condition (2.12) is thus precisely the condition for $\prod_j u_j^{m_j}$ to be a resonant term of the $i$th equation. For resonant terms, we instead choose $h_{i,\mathbf{m}_k}^{(k)} = 0$ and then $N_{u,i,\mathbf{m}_k}^{(k)} = N_{q,i,\mathbf{m}_k}^{(k)}$ and this term remains in the simplified normal form.

The normal form method then carries out the above procedure successively for each value of $k$, starting with $k=2$. For each $k$-value, we step over all possible indices $\{i, \mathbf{m}_k\}$ and choose the coefficient $h_{i,\mathbf{m}_k}^{(k)}$ to satisfy (2.11), unless $\{i, \mathbf{m}_k\}$ satisfies the resonance condition (2.12) in which case we choose $h_{i,\mathbf{m}_k}^{(k)} = 0$. Note that at each level $k$, the near-identity transformations used to remove the non-resonant terms introduces new terms with powers greater than $k$, so before proceeding to level $k+1$, the new system must be fully calculated. The value of $k>2$ at which we stop this process is called the *degree* of the normal form. Such normal forms are also sometimes called *truncated normal forms*.

### (c) Unfolding, truncation and dynamics of the normal forms

The notion of a normal form is often applied to situations where there is a local bifurcation of the equilibrium point $\mathbf{q}=0$ (e.g. [19] or [4]). In this case, we think of the system (2.4) as depending on parameters $\mu \in \mathbb{C}^r$ (typically the parameters are real, but they need not be) and $r$ corresponds to the codimension of the bifurcation which occurs at $\mu=0$. We then define an extended set of unknowns $\tilde{\mathbf{q}} = (\mathbf{x}, \mu) \in \mathbb{C}^{p+r}$ with the additional $r$ trivial equations $\dot{\mu} = 0$. Thus, we get

$$ \dot{\tilde{\mathbf{q}}} = \tilde{\Lambda}\tilde{\mathbf{q}} + \tilde{N}_{\tilde{\mathbf{q}}}(\tilde{\mathbf{q}}), \quad \text{where } \tilde{\Lambda} \in \mathbb{C}^{(p+r) \times (p+r)} = \begin{pmatrix} \Lambda & 0 \\ 0 & 0 \end{pmatrix} \qquad (2.13) $$
and all μ-dependence appears inside the first *p* components of the new nonlinear term $\tilde{N}_{\tilde{q}}$. The final *r* components of $\tilde{N}_{\tilde{q}}$ are precisely zero. The normal form method now proceeds exactly as defined above with *p* everywhere replaced by *p* + *r*.

We can then apply the normal form method as outlined above to the system (2.13). The resulting normal form is called an *unfolded normal form*. Sometimes judicious choices can be made for certain coefficients of a normal form, or by scaling arguments certain coefficients can be set to unity or zero. Such systems are sometimes called *hyper-normal forms* [6].

The normal form provides a way of simplifying the dynamics of a system near a non-hyperbolic equilibrium, and to classify the dynamics as belonging to one of a few pre-analysed cases according to the sign of certain critical coefficients in the normal form. One of the problems with using normal forms to provide more precise details is that, in general, they represent an asymptotic expansion only. Hence, regardless to which degree *k* they are truncated, they ignore the *beyond-all-orders* terms which can lead, for example, to transverse homoclinic and heteroclinic tangles. These will never be captured by the normal form, see, for example, Champneys & Kirk [20] and references therein for the case of the normal form of the codimension-two saddle-node/Hopf bifurcation (also known as the Gavrilov-Guckenheimer bifurcation). Occasionally (roughly speaking, when the dynamics of the normal form does not include any structurally unstable homo-/heteroclinic connections), an unfolded, truncated (hyper)-normal form can be shown to contain dynamics that are topologically equivalent to all possible structurally stable dynamics that can occur in a neighbourhood of the bifurcation point. Such a normal form is called a *versal unfolding* [21] or a *topological normal form* [19], an example of which is the normal form of a Hopf bifurcation, truncated after third-degree terms.

Consider the damped Duffing equation

$$ \ddot{x} + c\dot{x} + \omega^2 x + \alpha x^3 = 0. \quad (2.14) $$

We can include this in the normal form via an additional equation $\dot{c} = 0$. Carrying out the above steps, we first diagonalize the linear part by writing

$$ q_1 = \frac{1}{2\omega}(\omega x - i\dot{x}), \quad q_2 = \frac{1}{2\omega}(\omega x + i\dot{x}) \quad \text{and} \quad q_3 = c. $$

Then, we end up with a system of the form (2.13) for $\mathbf{q} = [q_1, q_2, q_3]^T$ with

$$ \Lambda = \begin{bmatrix} i\omega & 0 & 0 \\ 0 & -i\omega & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad \text{and} \quad N_q = \begin{bmatrix} \frac{i\alpha}{2\omega}(q_1+q_2)^3 - \frac{q_3}{2}(q_1-q_2) \\ -\frac{i\alpha}{2\omega}(q_1+q_2)^3 + \frac{q_3}{2}(q_1-q_2) \\ 0 \end{bmatrix}. $$

Now, from the form of the linearization of the extended system up to $O(3)$, we see from (2.12) that quadratic terms of the form $q_3q_i$ are resonant, as are cubic terms of the form $\alpha q_1^2 q_2$ in the $\dot{q}_1$ equation or of the form $\alpha q_1 q_2^2$ in the $\dot{q}_2$ one. We can remove all non-resonant cubic terms with a transformation of the form

$$ q_1 = \frac{\alpha}{8\omega}(-2u_1^3 + 3u_1u_2^2 + u_2^3) \quad \text{and} \quad q_2 = \frac{\alpha}{8\omega}(u_1^3 + 3u_2u_1^2 - 2u_2^3). $$

Hence, the third-degree normal form of the damped, unforced Duffing equation can be written as

$$ \dot{u}_1 = i\omega u_1 + \frac{3\alpha i}{2\omega} u_1^2 u_2 - \frac{c}{2}(u_1 - u_2) \quad \text{and} \quad \dot{u}_2 = -i\omega u_2 - \frac{3\alpha i}{2\omega} u_1 u_2^2 + \frac{c}{2}(u_1 - u_2), \quad \dot{c} = 0, \quad (2.15) $$

cf. Jezequel & Lamarque [7].

Jezequel & Lamarque [7] extended this normal form method in combination with the harmonic balance method to consider forced vibration problems of the form (1.1). This enabled them to predict resonant responses in the case that the forcing frequency $\Omega$ is close to one of the natural frequencies of the system. There remains a philosophical weakness that their method does not produce an approximation to the system itself, as Poincaré's envisaged, it simply uses the normal
form method to predict resonant (periodic) responses within a harmonic balance framework. This does not allow the complete dynamics close to certain multiple resonances to be explored, which in general will include quasi-periodic and possibly chaotic motion.

### 3. Normal forms for second-order mechanical systems

One of the weaknesses of using the general normal form method for vibration problems, as in Jezequel & Lamarque [7], is that this does not use the specific structure of the matrix A in (2.3). In view of this, Neild & Wagg [8] produced an extension to the method that can be applied directly to second-order differential equations of the form (1.1), preserving the decomposition into velocity and position variables. We term this the second-order normal forms, referring to the order of the differential equations rather than the level of accuracy achieved. As shown in §3c, this can lead to more accurate predictions of resonant amplitudes. However, their method still relies on the harmonic balance method and so one cannot appeal directly to the mathematical theory of normal forms. The key contribution of this paper then is to show this second-order normal form method can be re-derived in a similar spirit to the analysis of the previous section.

#### (a) Theoretical development

Consider a system of the form (1.1). We shall specifically treat the case where $r(t) = e^{i\Omega t}$ for some fixed forcing frequency $\Omega$, although the method is easily generalizable to quasi-periodic forcing. The first step is to diagonalize the system as much as possible. To do this, we set $x = \Phi q$, where $\Phi$ is a matrix of eigenvectors of $M^{-1}K$, giving

$$ \ddot{q} + \Lambda q + N_q(q, \dot{q}) = -\delta_1 \hat{D} \dot{q} + \delta_2 (\hat{P} r(t) + \text{c.c.}), \quad \text{with: } \Lambda = \text{diag}\{\omega_1^2, \omega_2^2, \dots, \omega_n^2\}, $$

and where $\hat{D}$, $\hat{P}$ and $N_q$ are the transformed damping matrix, forcing vector and nonlinearity, respectively, and $\omega_i$ are the linear frequencies of the system. Here, we have introduced two small parameters $\delta_1$ and $\delta_2$ which will play the role of unfolding parameters. Furthermore, we shall be interested in the critical situation where there is a resonance between the forcing frequency and one of the natural, undamped frequencies of the system, without loss of generality, $\omega_1$. That is, we assume that $\Omega = \omega_1 + \delta_3$, where $\delta_3$ is a third unfolding parameter. Finally, we suppose that there are no other linear resonances when $\delta = (\delta_1, \delta_2, \delta_3) = 0$; that is $\Omega \neq \omega_i$ for $i = 2, \dots, n$. Then, we can perform a further linear change of variables to remove the non-resonant forcing terms by writing

$$ \ddot{q} + \Lambda q + \sum_k N_q^{(k)}(q, \dot{q}) = P_q r(t) \xrightarrow{q = v + er} \ddot{v} + \Lambda v + \sum_k N_v^{(k)}(v, \dot{v}, r) = P_v r(t), \quad (3.1) $$

where $e_1 = 0$ and $e_k = P_k/(\omega_k^2 - \Omega^2)$ so that the forcing term is in the first, resonant equation only. Here, nonlinear terms have been expressed in monomial form.

We shall now use a tilde to represent an extended variable, so that $\tilde{v} \in \mathbb{C}^{n+4}$ is the vector $(v, r, \delta)$ and write the complete system in the form

$$ \ddot{v} + \tilde{\Lambda}\tilde{u} + \sum_k \tilde{N}_v^{(k)}(\tilde{v}, \dot{\tilde{v}}) = \tilde{P}_v r(t), \quad (3.2) $$

where the first $n$ components of $\ddot{\tilde{v}}$ are given by the transformed equation in (3.1) and the last four components by $\ddot{r} = -\Omega^2 r$ and $\ddot{\delta} = 0$. Here, the first $n$ components of $\tilde{N}_v \in \mathbb{C}^{n+4}$ are equal to the nonlinearity $N_v$ and the last four components are zero. Also,

$$ \tilde{\Lambda} \in \mathbb{C}^{n+4} \times \mathbb{C}^{n+4} = \text{diag}\{\omega_1^2, \omega_2^2, \dots, \omega_n^2, \Omega^2, 0, 0, 0\}. \quad (3.3) $$

Note that forcing, damping and frequency mismatch are all now thought of as nonlinear terms, because they contain terms that involve an unfolding parameter $\delta_i$ times a state variable $u$ or $r$ or their derivatives, and we think of the $\delta$s as state variables too. Observe how the system (3.2) looks remarkably similar to the unfolded first-order normal form (2.13) apart from the second
derivatives on the left-hand side. This difference will be important when we try to compute the normal form for such a system.

Now we focus on transforming the equations of motion to the form

$$ \ddot{\mathbf{v}} + \mathcal{A}\tilde{\mathbf{u}} + \sum_k \tilde{N}_v^{(k)}(\tilde{\mathbf{v}}, \dot{\tilde{\mathbf{v}}}) = \tilde{P}_v r \xrightarrow{k} \tilde{\mathbf{h}}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}, r) \rightarrow \ddot{\mathbf{u}} + \mathcal{A}\tilde{\mathbf{u}} + \sum_k \tilde{N}_{u}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}, r) = \tilde{P}_u r, \quad (3.4) $$

where, as with $\tilde{\mathbf{v}}, \tilde{\mathbf{u}} = (\mathbf{u}, r, \delta)$. As before, we suppose that all non-resonant nonlinear terms have been removed up to degree $k-1$ using the transformation to remove terms at $O(k)$. Eliminating $\tilde{\mathbf{v}}$ leads to the $O(k)$ equation

$$ \tilde{N}_{u}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}) - \tilde{N}_{v}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}) = \mathcal{A}\tilde{\mathbf{h}}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}) + \frac{d^2}{dt^2}\tilde{\mathbf{h}}^{(k)}(\tilde{\mathbf{u}}, \dot{\tilde{\mathbf{u}}}). \quad (3.5) $$

To proceed, we next complexify the system by writing

$$ u_i = u_i + \bar{u}_i, \quad \dot{u}_i = i\omega_i(u_i - \bar{u}_i) \quad \text{and} \quad r = \tau + \bar{\tau}, \quad (3.6) $$

where $u_i$ and $\bar{u}_i$ are the $i$th elements in $\mathbf{u}$ and $\mathbf{u}$, respectively, and the bar indicates the complex conjugate. This results in $2n+5$ arguments for $\tilde{\mathbf{h}}^{(k)}$, $\tilde{N}_u$ and $\tilde{N}_v$ which we temporarily call $\mathbf{z}$ such that $\mathbf{z} = (\mathbf{u}, \tilde{\mathbf{u}}, \tau, \bar{\tau}, \delta_1, \delta_2, \delta_3)$. As with the first-order case, we wish to write (3.5) in indicial form and solve the corresponding homological equation. To do this, we again assume a particular component for the pre- and post-transformed nonlinear terms and the transform term of

$$ \tilde{N}_{v,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j}, \quad \tilde{N}_{u,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} \quad \text{and} \quad \tilde{h}_{i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j}, \quad (3.7) $$

respectively. Here, we define the augmented vector $\tilde{\mathbf{m}}_k$ of length $2n+5$, corresponding to the $2n+5$ arguments, $\mathbf{z}$, and may be written using the notation $\tilde{\mathbf{m}} = (m_1, \dots, m_n, m_{-1}, \dots, m_{-n}, m_r, m_{-r}, m_{\delta_1}, m_{\delta_2}, m_{\delta_3})$, with the $j$th element referred to as $\tilde{m}_j$. Using this indicial representation of the $n+4$ element vectors $\tilde{N}_v^{(k)}$, $\tilde{N}_u^{(k)}$ and $\tilde{h}^{(k)}$, (3.5) may be written in homological form as

$$ \tilde{N}_{u,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} - \tilde{N}_{v,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} = \tilde{\lambda}_i \tilde{h}_{i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} + \tilde{h}_{i,\tilde{\mathbf{m}}_k}^{(k)} \sum_{l=1}^{2n+5} \frac{\partial}{\partial z_l} \left( \sum_{l=1}^{2n+5} \frac{\partial}{\partial z_l} \left( \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} \right) \right) \gamma_l z_l \\ = \tilde{h}_{i,\tilde{\mathbf{m}}_k}^{(k)} \left( \lambda_i + \left[ \sum_{j=1}^{2n+5} \tilde{m}_j \gamma_j \right]^2 \right) \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j}. \quad (3.8) $$

Here, $\tilde{\lambda}_i$ is the $i$th diagonal element in $\tilde{\mathcal{A}}$, see (3.3), and $\tilde{z}_i = \gamma_i z_i$ such that $\gamma_i$ is given by

$$ \gamma_i = \begin{cases} i\omega_i & \text{for: } 1 \le i \le n, \\ -i\omega_{i-n} & \text{for: } n+1 \le i \le 2n, \\ i\Omega & \text{for: } i=2n+1, \\ -i\Omega & \text{for: } i=2n+2, \\ 0 & \text{for: } 2n+3 \le i \le 2n+5. \end{cases} \quad (3.9) $$

By considering the form of $\tilde{\mathcal{A}}$ and $\gamma_i$, equation (3.8) can be simplified to give for $1 \le i \le n$

$$ N_{u,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} - N_{v,i,\tilde{\mathbf{m}}_k}^{(k)} \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j} = h_{i,\tilde{\mathbf{m}}_k}^{(k)} \left( \omega_i^2 - \left[ (m_r - m_{-r})\Omega + \sum_{j=1}^{n} (m_j - m_{-j})\omega_j \right]^2 \right) \prod_{j=1}^{2n+5} z_j^{\tilde{m}_j}, \quad (3.10) $$
where, for these values of $i$, $N_{v,i,\bar{m}_k}^{(k)} = \tilde{N}_{v,i,\bar{m}_k}^{(k)}$ and we define $N_{u,i,\bar{m}_k}^{(k)} = \tilde{N}_{u,i,\bar{m}_k}^{(k)}$ and $h_{u,i,\bar{m}_k}^{(k)} = \tilde{h}_{u,i,\bar{m}_k}^{(k)}$. For $n+1 \le i \le n+4$, $\tilde{N}_{v,i,\bar{m}_k}^{(k)} = 0$, so we set $\tilde{N}_{u,i,\bar{m}_k}^{(k)} = \tilde{h}_{i,\bar{m}_k}^{(k)} = 0$.

The homological equation can be inverted and solved for $h_{i,\bar{m}_k}^{(k)}$ provided we avoid the resonant terms for which

$$ \omega_i^2 = \left[ (m_r - m_{-r})\Omega + \sum_{j=1}^{n} (m_j - m_{-j})\omega_j \right]^2. \quad (3.11) $$

Note that $m_r$ or $m_{-r} = 1$ with $m_{\delta,1} = 1$ and $i=1$ automatically leads to resonance because of the fundamental resonance we have assumed between the forcing and natural frequency $\omega_1$; that is, $\omega_1 = \Omega$. But there will in general be other resonances; for example, whenever $m_i = m_{-i} + 1$ as with the normal form of the Duffing equation. As we shall see in later examples, there can be further internal resonances if $\omega_i = l\omega_j$ for some integers $i, j, k$ and $l$.

## (b) Practical implementation

Having shown in principle that a second-order version of the normal form can be derived for structural vibration systems, we now turn to a practical implementation of it, as discussed in detail in Neild & Wagg [8] and Neild [22]. This implementation avoids the need to explicitly form the extended $2n+5$ dimensional state vectors. The key difference in this implementation is that, rather than seeking a specific normal form, we seek to project the system onto a nonlinear equivalent of normal modes, in a harmonic-balance-like approach. As we shall see this approach relies on computation of exactly the same normal-form coefficients as in the previous subsection. So, we specifically suppose that in the complexification step (3.6), we seek a solution of the form

$$ u_i = u_i + \bar{u}_i, \quad \text{with } u_i = (U_i e^{-i\phi_i/2}) e^{i\omega_i r}, $$

where $\omega_{ri}$ is the response frequency in the $i$th modal coordinate direction. Here, $\omega_{ri}$ is, in general, a nonlinear function of the variables $u$ that must be determined as part of the calculation. It corresponds, for an unforced system to the nonlinear natural frequency for moderate amplitude periodic orbits. This concept can be made precise using the Liapunov centre theorem $\omega_{ri} = 2\pi/T$, where $T$ is the period for a given amplitude periodic orbit in the NNM invariant manifold (in the sense of Shaw & Pierre [9]) tangent to the $\omega_i$ eigenspaces. In the presence of forcing, the response frequency of a resonant mode is taken to be $\omega_{ri} = \Omega$.

Now the analysis is the same as before using (3.2) and then the equivalent to (3.1), but without using the extended states. The equivalent of the expression (3.5) for the $k=1$ case reads

$$ N_v^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) + \ddot{\mathbf{h}}^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) + \Gamma \mathbf{h}^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) = N_u^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r). \quad (3.12) $$

Here, we have introduced diagonal matrix $\Gamma$ with the $i$th diagonal element being $\omega_{ri}^2$. To first approximation, $k=1$, this equals $\Lambda$ as $\omega_{ri}^2 \approx \omega_i^2$. Making this substitution has no algebraic effect on the lowest degree approximation to the normal form, but results in more accurate prediction of the harmonics of the responses [23]. As a consequence, the $k>1$ equations contain $\Gamma - \Lambda$ correction terms (although an acceptable approximation is normally achieved without evaluating these equations).

As before, after complexification of the equations, we now seek to solve (3.12). This can be for each of the terms in the $O(k)$ nonlinearity $N_{u}^{(k)}$. Now, instead of introducing multi-indices $\{i, m_k\}$, we simply place all the combinations of states that exist in $N^{(k)}$ into vector $\mathbf{u}^*$ (of length $\ell$) and introduce coefficient matrices $\mathbf{h}^*, \mathbf{n}^*$ and $\mathbf{n}_{u}^*$ in $\mathbb{C}^{n \times \ell}$, so that

$$ \left. \begin{array}{l} N_v^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) = \mathbf{n}^* \mathbf{u}^*(\mathbf{u}, \bar{\mathbf{u}}, \tau, \bar{\tau}), \quad \mathbf{h}^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) = \mathbf{h}^* \mathbf{u}^*(\mathbf{u}, \bar{\mathbf{u}}, \tau, \bar{\tau}) \\ N_u^{(1)}(\mathbf{u}, \dot{\mathbf{u}}, r) = \mathbf{n}_u^* \mathbf{u}^*(\mathbf{u}, \bar{\mathbf{u}}, \tau, \bar{\tau}). \end{array} \right\} \quad (3.13) $$

and
This is an alternative to the arbitrary product notation used earlier, see (3.7) for the indicial form of this. Using the new representation, (3.12) may be rewritten as

$$ \mathbf{n}_{u}^{*} - \mathbf{n}^{*} = \tilde{\mathbf{h}}^{*} \quad \text{where, in indicial form: } \tilde{h}_{i,\ell}^{*} = -\beta_{i,\ell}^{*} \mathbf{h}_{i,\ell}^{*}. \qquad (3.14) $$

Considering the indicial form, it can be shown that

$$ \beta_{i,\ell}^{*} = \left( [m_r - m_{-r}]\Omega + \sum_{j=1}^{n} [(m_j - m_{-j})\omega_{rj}] \right)^2 - \omega_{r\ell}^2 \quad (3.15) $$

where the form of the ℓth term in **u**\* is given by

$$ u_{\ell}^{*} = \mathbf{v}^{m_r} \bar{\mathbf{r}}^{m_{-r}} \prod_{j=1}^{n} (\mathbf{u}_{k}^{m_{k}} \bar{\mathbf{u}}_{k}^{m_{-k}}). \qquad (3.16) $$

Noting the minus sign in (3.14), introduced to maintain consistency of the definition of $\beta^*$ with previous publications, (3.15) is exactly equivalent to (3.10). The appearance of $\omega_{r\ell}$ rather than $\omega_i$ is due to the introduction of $\Gamma$.

Considering the indicial version of (3.14), $n_{u,i,\ell}^* + \beta_{i,\ell}^* \mathbf{h}_{i,\ell}^* = n_{i,\ell}^*$, if $\beta_{i,\ell}^*$ is zero, it is resonant and must be retained in the post-transformed nonlinear vector, otherwise the nonlinear term can be removed using the transform such that

$$ \left. \begin{array}{l} \text{if } \beta_{i,\ell}^* = 0: \quad n_{u,i,\ell}^* = n_{i,\ell}^*, \quad h_{i,\ell}^* = 0 \quad (\text{a resonant term}) \\ \text{otherwise: } \quad n_{u,i,\ell}^* = 0, \quad h_{i,\ell}^* = \frac{n_{i,\ell}^*}{\beta_{i,\ell}^*} \quad (\text{a non-resonant term}). \end{array} \right\} \qquad (3.17) $$

Having identified $\mathbf{h}^*$ and $\mathbf{n}_{u,\ell}^*$, and hence $\mathbf{h}^{(1)}$ and $N_u^{(1)}$ using (3.13), we can approximate the transform and transformed dynamic equations to

$$ \mathbf{v} = \mathbf{u} + \mathbf{h}^{(1)}(\mathbf{r}, \bar{\mathbf{r}}, \mathbf{r}), \quad \ddot{\mathbf{u}} + \Lambda \dot{\mathbf{u}} + N_u^{(1)}(\mathbf{r}, \bar{\mathbf{r}}, \mathbf{r}) = \tilde{P}_u \mathbf{r}. \qquad (3.18) $$

The approximation here is that we have used only the $k=1$ terms, $\mathbf{h}^{(1)}$ and $N_u^{(1)}$ rather than a summation over all $k$ as was used in (3.4). To refine these expressions, the higher $k$ terms can be included using a similar approach (e.g. [23]), however this is not normally necessary.

In the case of forced systems, when the forcing is close to resonance for all the modes, the forcing transform (3.1) reduces to $\mathbf{q} = \mathbf{v}$ and the nonlinear transform $\mathbf{v} \to \mathbf{u}$ remains the same as the unforced case, this is discussed further in the example that follows. When non-resonant forcing occurs in any mode then additional terms appear in $\mathbf{u}^*$ due to the forcing transform— for examples which include non-resonant forcing, see Wagg & Neild [24] and Neild & Wagg [8], in which a two-mode forced system is assessed. Wagg & Neild [24] also contains a discussion and comparison of the second-order normal forms alongside the harmonic balance, averaging and multiple scales techniques for the vibration of a forced single mode system. The relationship between the forced response and the NNMs of the unforced system is discussed in Hill et al. [25].

Two examples are now considered, the first is the Duffing oscillator which will be used to compare results from the first-order, state-space, variant of the normal forms with the second-order, oscillator equation, variant. Both responses will be considered using just the $k=1$ terms. This example will also be used to demonstrate the forced and damped case. The second example considers the NNMs for a two-mode model of a cable under free vibration. Both the resonant and the harmonic responses are calculated to yield the response frequencies and NNM mode-shapes.
(c) Example: the Duffing equation

Let us consider the unforced, undamped Duffing equation

$$\ddot{x} + \omega_n x + \alpha x^3 = 0. \quad (3.19)$$

As this system has one degree-of-freedom $x = \mathbf{x} = \mathbf{q} = \mathbf{v}$. Using the matrix formulation, we keep the cubic term $\alpha x^3$, so that $N_v^{(3)}(\mathbf{u}) = \alpha u_1^3$. Using (3.13), $\mathbf{n}^*$ and $\mathbf{u}^*$ may be defined, then $\beta^*$ can be calculated using (3.15) and from this the transformed nonlinear terms and the transform can be found using (3.17). This gives

$$\mathbf{n}^* = \begin{bmatrix} \alpha \\ 3\alpha \\ 3\alpha \\ \alpha \end{bmatrix}^T, \quad \mathbf{u}^* = \begin{bmatrix} u_1^3 \\ u_1^2 \bar{u}_1 \\ u_1 \bar{u}_1^2 \\ \bar{u}_1^3 \end{bmatrix} \rightarrow \beta^* = \omega_{r1}^2 \begin{bmatrix} 8 \\ 0 \\ 0 \\ 8 \end{bmatrix}^T \rightarrow \mathbf{n}_u^* = \begin{bmatrix} 0 \\ 3\alpha \\ 3\alpha \\ 0 \end{bmatrix}^T, \quad \mathbf{h}^* = \frac{\alpha}{8\omega_{r1}^2} \begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \end{bmatrix}^T,$$

Using (3.1), (3.4) and (3.13) allows the transformed equation of motion along with the transform to be written as

$$\ddot{u}_1 + \omega_{n1}^2 u_1 + 3\alpha(u_1^2 \bar{u}_1 + u_1 \bar{u}_1^2) = 0, \quad x = q_1 = u_1 + \frac{\alpha}{8\omega_{r1}^2}(u_1^3 + \bar{u}_1^3),$$

respectively. As the response for the ith mode will just be at the resonant frequency $\omega_{ri}$, we can write the steady-state solution

$$u_i = u_i + \bar{u}_i = \left(\frac{U_i}{2} e^{-i\phi_i}\right) e^{i\omega_{ri}t} + \left(\frac{U_i}{2} e^{i\phi_i}\right) e^{-i\omega_{ri}t}. \quad (3.20)$$

Hence,

$$x = U_1 \cos(\omega_{r1} t) + \frac{\alpha}{32\omega_{r1}^2} U_1^3 \cos(3\omega_{r1} t), \quad \text{with: } \omega_{r1}^2 = \omega_1^2 + \frac{3\alpha}{4} U_1^2, \quad (3.21)$$

where we have defined t such that $\phi_1 = 0$. Note that as the harmonic terms were removed using the near-identify transform, $U_1$ conveniently represents the amplitude of the resonant response.

Considering the case where the system is damped and forced at a frequency close to resonance, such that $\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n x + \alpha x^3 = R \cos(\Omega t)$, the forcing transform is $\mathbf{q} = \mathbf{u}$. As the forcing is near-resonant, the response frequency is set as $\omega_{r1} = \Omega$ and the near-identity transform is the same as that derived above, it is unaffected by the introduction of the forcing or the damping (which can be thought of as a resonant term). The resulting resonant dynamics and transform are

$$\ddot{u}_1 + 2\zeta\omega_{n1}u_1 + \omega_{n1}^2 u_1 + 3\alpha(u_1^2 \bar{u}_1 + u_1 \bar{u}_1^2) = R \cos(\Omega t) \quad \text{and} \quad x = q_1 = u_1 + \frac{\alpha}{8\omega_{r1}^2}(u_1^3 + \bar{u}_1^3),$$

respectively. Using the steady-state solution (3.20) and balancing the complex exponential terms allows the amplitude relationship $[(\omega_n^2 - \Omega^2)U_1 + (3\alpha/4)U_1^3]^2 + [2\zeta\omega_n\Omega U_1]^2 = R^2$ and phase relationship $\tan(\phi_1) = 2\zeta\omega_n\Omega / (\omega_n^2 - \Omega^2 + (3\alpha/4)U_1^2)$ to be found. Once $U_1$ and $\phi_1$ have been found for a given forcing, then the response, including harmonics, can be calculated using the transform equation and (3.20).

The unforced solution can be compared with the solution derived using the first-order version of the normal form technique in which the state-space, rather than oscillator, form of the equations are used. Here, the first step is to rewrite the equation of motion in state-space form using $x = (\mathbf{x}, \dot{\mathbf{x}})$ and then to apply a linear transform giving

$$\dot{\mathbf{x}} = \begin{bmatrix} 0 & 1 \\ -\omega^2 & 0 \end{bmatrix} \mathbf{x} + \begin{bmatrix} 0 \\ -\alpha x_1^3 \end{bmatrix} \rightarrow \dot{\mathbf{q}} = \begin{bmatrix} i\omega & 0 \\ 0 & -i\omega \end{bmatrix} \mathbf{q} + \frac{i\alpha}{2\omega} \begin{bmatrix} (q_1 + q_2)^3 \\ -(q_1 + q_2)^3 \end{bmatrix},$$

where, as the system is written in first-order form, we have

$$x = \begin{bmatrix} 1 & 1 \\ i\omega_n & -i\omega_n \end{bmatrix} \mathbf{q}.$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Duffing oscillator example with $\omega = 1$ and $\alpha = 0.5$: (a) backbone curves and (b) forced response curves when $\zeta = 0.01$ and $R = 0.1$. The solid line shows the exact computation using AUTO; the dashed line shows the result of the second-order normal form calculation; and the dashed-dotted line the result of the first-order normal form.</figcaption></figure>

Applying the normal form transformation in a similar way to the second-order matrix formulation, now using $\mathbf{u} = [\mathbf{u} \ \bar{\mathbf{u}}]^T$ gives

$$x = U_1 \left( 1 - \frac{3\alpha U_1^2}{16\omega_{n1}^2} \right) \cos(\omega_{r1} t) + \frac{\alpha}{32\omega_{n1}^2} U_1^3 \cos(3\omega_{r1} t), \quad \text{with: } \omega_{r1} = \omega_{n1} + \frac{3\alpha}{8\omega_{n1}} U_1^2 \quad (3.22)$$

[7,12]. Note that here $U_1$ does not fully represent the amplitude of the resonant response. This is because, when identifying the transform, only the $e^{i\omega_{r1}t}$ terms are kept in the first equation of motion in $\mathbf{U}$ with $e^{-i\omega_{r1}t}$ terms being represented in the transform and vice versa for the second equation of motion [8].

The response (or backbone) curves for both the second- and the first-order variants of the normal forms, (3.21) and (3.22), respectively, are shown in figure 2a. For comparison, the solution derived using numerical continuation, using AUTO [26], in which no assumptions are made regarding the smallness of terms, is also shown. It can be seen that at low amplitudes, corresponding to weaker nonlinear terms, there is good agreement in the prediction of the response frequency, however as the amplitude increases the higher power terms become significant, firstly for the first-order normal form approximation at about 0.6 and then also for the second-order version at about 1.8. Figure 2b compares the second-order normal form and the AUTO predictions of a forced response, again it can be seen that good agreement is achieved.

## 4. Application to nonlinear normal modes

We now consider how the application of normal forms to unforced, undamped second-order differential equations can be used to calculate the NNMs of a system. Second-order normal forms can be used to derive the backbone curves of a system via the resonant responses (e.g. Wagg & Neild [24] or Hill et al. [27]). Here, we extend this to finding the NNM mode shapes, which requires both the resonant response and also the harmonics captured in $\mathbf{h}$. While the technique is general, to facilitate the discussion we apply it to a well-known nonlinear system—the dynamics of a taut cable. Specifically, we will consider the interactions between the first out-of- and in-plane modes (we use mode to indicate the modes of the linearized system). In doing this, to calculate the associated NNMs, we will derive the harmonics excited in all modes due to resonant responses in the first out-of- and in-plane modes.

To examine the dynamics of a taut cable, we use the modal equations of motion derived by Warnitchai et al. [28] and discussed in detail in Wagg & Neild [24], which accounts for
gravitational sag and tension effects. The motion of the cable out-of- and in-plane (where in-plane relates to the plane in which gravitational sag occurs) are represented as

$$v(x, t) = \sum_{n} \phi_{n}(x)y_{n}(t) \quad \text{and} \quad w(x, t) = w_{s}(x) + \sum_{n} \psi_{n}(x)z_{n}(t), \qquad (4.1)$$

where $x$ is the distance from one of the supports and the other support is positioned at $x = \ell$. Here, $\phi_n(x)$ and $\psi_n(x)$ are the nth out-of- and in-plane linear mode shapes, respectively, $y_n(t)$ and $z_n(t)$ represent the modal contributions for the nth modes and $w_s(x)$ captures the static sag in the cable (noting that $v$ and $w$ are defined as zero on the chord line between the supports). The dynamics for the nth out-of- and in-plane modes may be expressed as

$$\left.
\begin{aligned}
& \ddot{y}_n + \omega_{yn}^2 y_n + \sum_k \frac{\nu_{nk}}{m} y_n(y_k^2 + z_k^2) + \sum_k 2 \frac{\beta_{nk}}{m} y_n z_k = 0 \\
& \ddot{z}_n + \omega_{zn}^2 z_n + \sum_k \frac{\nu_{nk}}{m} z_n(y_k^2 + z_k^2) + \sum_k 2 \frac{\beta_{nk}}{m} z_n z_k + \sum_k \frac{\beta_{kn}}{m} (y_k^2 + z_k^2) = 0,
\end{aligned}
\right\} \quad (4.2)$$

and respectively. Here, modal damping and external forcing terms have been removed from the original derivation in line with investigating the backbone curves of the system. Parameters $m$, $\beta_{ij}$ and $v_{ij}$ are given in Warnitchai et al. [28] and Gonzalez-Buelga et al. [29], but importantly $\beta_{ij}$ is zero for even $j$. The natural frequencies of the out-of-plane modes are proportional to the mode number $n$, and for even $n$, natural frequencies for the in-plane modes match those of the out-of-plane ones. For odd $n$, the in-plane natural frequencies are slightly higher than the out-of-plane ones, again see [28,29]. The in-plane axis is defined as positive down, and so gravitational sag is positive.

Gonzalez-Buelga et al. [29] considered these equations in terms of internal resonance showing that 1:1 resonance occurs between the second out-of- and in-plane modes and Macdonald et al. [30] generalized this for the nth out-of- and in-plane modes. In addition, they both show that 2:1 resonance can occur, but only for the case where the cable is inclined. Hill et al. [27] identified the backbone curves when the system is reduced to the first out-of- and in-plane modes. Here, we build on this work by deriving algebraic expressions for the NNMs associated with these backbone curves, which requires not only the resonant responses of the two modes but also the harmonic response in all modes. To do this, we first calculate the backbone curves along with expressions that capture the harmonics contained in the response. These are then used to find the NNMs in terms of just the two first modes and finally the additional modal contribution due to harmonics in other modes is added to give the full NNM expressions.

### (a) Resonant equation of motion and harmonic response

Letting $\mathbf{q} = (y_1 \ z_1)^T$ the modal equation of motion for the reduced two-mode model may be written in the form $\ddot{\mathbf{q}} + \Lambda\mathbf{q} + N_q(\mathbf{q}) = 0$, an unforced version of the left-hand expression in (3.1), where

$$\mathbf{q} = \begin{pmatrix} q_1 \\ q_2 \end{pmatrix}, \quad \Lambda = \begin{bmatrix} c \omega_{y1}^2 & 0 \\ 0 & \omega_{z1}^2 \end{bmatrix} \quad \text{and} \quad \mathbf{N}_q^{(1)} = \begin{pmatrix} c \left( \frac{\nu_{11}}{m} \right) (q_1^3 + q_1 q_2^2) + 2 \left( \frac{\beta_{11}}{m} \right) (q_1 q_2) \\ \left( \frac{\nu_{11}}{m} \right) (q_2^3 + q_1^2 q_2) + \left( \frac{\beta_{11}}{m} \right) (q_1^2 + 3q_2^2) \end{pmatrix}.$$

Here, taking $N_q(\mathbf{q}) = \sum_k N_q^{(k)}(\mathbf{q})$, all the nonlinear terms have been placed in $N_q^{(1)}$. As there is no forcing, $\mathbf{v} = \mathbf{q}$ and so $N_v^{(1)}(\mathbf{v}) = N_q^{(1)}(\mathbf{q})$. To calculate the nonlinear transform, see (3.4), $N_v^{(1)}(\mathbf{u})$ is considered. Rewriting this nonlinear vector in terms of $\mathbf{u}$, where, for the kth coordinate $u_k = u_k + \bar{u}_k$, gives

$$N_v^{(1)}(\mathbf{u}) = \frac{\nu_{11}}{m} \begin{pmatrix} (u_1 + \bar{u}_1)^3 + (u_1 + \bar{u}_1)(u_2 + \bar{u}_2)^2 \\ (u_2 + \bar{u}_2)^3 + (u_1 + \bar{u}_1)^2(u_2 + \bar{u}_2) \end{pmatrix} + \frac{\beta_{11}}{m} \begin{pmatrix} 2(u_1 + \bar{u}_1)(u_2 + \bar{u}_2) \\ (u_1 + \bar{u}_1)^2 + 3(u_2 + \bar{u}_2)^2 \end{pmatrix}.$$
This can now be expressed in matrix form $\mathbf{n}_v^* u^*(\mathbf{u}, \bar{\mathbf{u}}, \bar{\mathbf{r}}, \bar{\mathbf{v}})$, see (3.13) along with (3.6), and the matrix $\boldsymbol{\beta}^*$ can be derived using (3.15), where it can be assumed that, due to the closeness of the natural frequencies, we can write the response frequencies as $\omega_{r1} = \omega_{r2}$. Using (3.17), the nonlinear vector in the transformed equation of motion and the transformation vector can be written as

$$N_u^{(1)} = \frac{\nu_{11}}{m} \left( \begin{array}{c} 3u_1 \bar{u}_1 u_1 + 2u_2 \bar{u}_2 u_1 + u_1 \bar{u}_2^2 + \bar{u}_1 u_2^2 \\ 3u_2 \bar{u}_2 u_2 + 2u_1 \bar{u}_1 u_2 + u_1^2 \bar{u}_2 + \bar{u}_1^2 u_2 \end{array} \right) \quad (4.3)$$

and

$$\mathbf{h}^{(1)} = \frac{\nu_{11}}{8m\omega_{r1}^2} \begin{pmatrix} u_1^3 + \bar{u}_1^3 + u_1 u_2^2 + \bar{u}_1 \bar{u}_2^2 \\ u_2^3 + \bar{u}_2^3 + u_1^2 u_2 + \bar{u}_1^2 \bar{u}_2 \end{pmatrix} + \frac{\beta_{11}}{3m\omega_{r1}^2} \begin{pmatrix} 2(u_1 u_2 + \bar{u}_1 \bar{u}_2) - 6(u_1 \bar{u}_2 + \bar{u}_1 u_2) \\ u_1^2 + \bar{u}_1^2 + 3(u_2^2 + \bar{u}_2^2) - 6u_1 \bar{u}_1 - 18u_2 \bar{u}_2 \end{pmatrix}, \quad (4.4)$$

respectively. The modal response of the system can now be written as $\mathbf{q} = \mathbf{v} = \mathbf{u} + \mathbf{h}^{(1)}$, where $\mathbf{h}^{(1)}$ captures the harmonic content of the response.

Using the steady-state solution (3.20), the transformed equations of motion, $\ddot{\mathbf{u}} + \Lambda \mathbf{u} + N_u^{(1)} = 0$, can be written as

$$\{\omega_{y1}^2 - \omega_{r1}^2 + \frac{\nu_{11}}{4m}[3U_1^2 + (2+p)U_2^2]\} U_1 = 0 \quad (4.5)$$

and

$$\{\omega_{z1}^2 - \omega_{r1}^2 + \frac{\nu_{11}}{4m}[(2+p)U_1^2 + 3U_2^2]\} U_2 = 0, \quad (4.6)$$

where $p = \cos(2(\phi_1 - \phi_2))$ with the condition $\sin(2(\phi_1 - \phi_2)) = 0$ to ensure (4.3) is real. This allows possible solutions $p = \pm 1$. It can be shown that only $p = -1$ gives physically meaningful solutions (as $\omega_{z1} \neq \omega_{y1}$) which means that when both linear modes are present, they are $\pm 90^{\circ}$ out-of-phase, see Hill et al. [27] for a more detailed discussion of this.

## (b) Backbone curves and modal response

Taking the physically meaningful $p = -1$ case, there are two semi-trivial solutions for (4.5) and (4.6) in which only one of the two linear modes is resonant

$$S1: \quad \omega_{r1}^2 = \omega_{y1}^2 + \frac{3\nu_{11}}{4m} U_1^2, \quad U_2 = 0 \quad (4.7)$$

and

$$S2: \quad \omega_{r1}^2 = \omega_{z1}^2 + \frac{3\nu_{11}}{4m} U_2^2, \quad U_1 = 0, \quad (4.8)$$

and two further solutions exist in which both linear modes are resonant

$$S3^{\pm}: \quad \omega_{r1}^2 = \omega_{y1}^2 + \frac{\nu_{11}}{4m}(3U_1^2 + U_2^2), \quad U_2^2 = U_1^2 - \frac{2m}{\nu_{11}}(\omega_{z1}^2 - \omega_{y1}^2). \quad (4.9)$$

As $\omega_{z1} > \omega_{y1}$, it can be seen that valid solutions for $S3$ only exist if $U_1 \ge 2m(\omega_{z1}^2 - \omega_{y1}^2)/\nu_{11}$. The point at which $U_1 = 2m(\omega_{z1}^2 - \omega_{y1}^2)/\nu_{11}$ results in $U_2 = 0$ and lies on backbone curve $S1$, see (4.8). Hence these $S3$ solutions (there are two solutions relating to the relative phase of the two linear modes being $\pm 90^{\circ}$, solutions $S3^{\pm}$) are branches from $S1$ following a bifurcation at $U_2 = 0$. This bifurcation from the $S1$ solution is shown in figure 3 using the parameters defined in Gonzalez-Buelga et al. [29]. Note that the $S3^{\pm}$ backbone curves lie on top of each other in this projection.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Backbone curves S1 and S3± for a cable with inset panels showing the motion in projection $q_1 = y_1$ versus $q_2 = z_1$. The inset panel for the response near the bifurcation contains two responses; the response line (a) corresponding to a point on S1s just below the bifurcation and the response loop (b) corresponding to a point on S3+ just above the bifurcation. Subscripts s and u indicate that backbone curve S1 is stable below and unstable above the bifurcation.</figcaption></figure>

Using the transform equation, $q = v = u + h^{(1)}$, and the steady-state solution (3.20) gives the response of the modal coordinates

$$ y_1 = U_1 \cos(\omega_{r1} t) + \frac{\nu_{11}}{32m\omega_{r1}^2} U_1(U_1^2 - U_2^2) \cos(3\omega_{r1} t) \\ + \frac{\beta_{11}}{3m\omega_{r1}^2} U_1 U_2 \sin(2\omega_{r1} t) \quad (4.10) $$

and

$$ z_1 = U_2 \sin(\omega_{r1} t) + \frac{\nu_{11}}{32m\omega_{r1}^2} U_2(U_1^2 - U_2^2) \sin(3\omega_{r1} t) \\ + \frac{\beta_{11}}{6m\omega_{r1}^2} (U_1^2 - 3U_2^2) \cos(2\omega_{r1} t) - \frac{\beta_{11}}{2m\omega_{r1}^2} (U_1^2 + 3U_2^2), $$

when the phase difference is $\phi_2 - \phi_1 = 90^\circ$ and $\phi_1$ has been set to zero.

The S2 solution consists of purely in-plane resonant motion as $U_1 = 0$. In this case, the response given by (4.10) reduces to purely in-plane motion. Harmonics exist in this motion including a static displacement which, as it is negative, reduces the apparent sag of the cable during resonance. For the S1 solution, which is resonant purely in the out-of-plane direction as $U_2 = 0$, there are also harmonic components in the out-of-plane response. In addition, there is an amplitude-dependent in-plane non-resonant response. For the S1 solution example of the motion of the mid-span of the cable (where the mode shapes are unity) in the $y_1$ versus $z_1$ plane are shown as in figure 3 as inset plots. Note that in contrast to Hill et al. [27], these inset panels include the harmonic components as well as the resonant ones. Examples of the more complex motion of the S3 solution, in which there is 90° out-of-phase resonant motion in both planes, are also shown in figure 3. Here, it can be seen that away from the bifurcation (indicated with a dot), the cable is
whirling in a near elliptical shape—it would be exactly circular if $\omega_{y1} = \omega_{z1}$, albeit offset from the origin in the vertical direction.

These inset panels show the motion at the mid-span assuming that the harmonics of the cable response are limited to the first out-of- and in-plane modes. Once additional harmonics from other modes are included in the response, the full NNM response can be derived. This is done in §4c for the NNMs that contain resonant responses solely in either or in both the first out-of- or in-plane modes.

### (c) Additional modal contributions

In the previous analysis, only the first out-of- and in-plane modes have been considered. It has been shown that in-plane motion is present throughout (figure 3), even for the S1 solution, where there is no resonant in-plane response, but there are twice- and zero-frequency components. By inspection of (4.2), it can be seen that this is due to the term $\beta_{11}y_1^2/m$ that arises from the $\sum_k(\beta_{kn}/m)(y_k^2 + z_k^2)$ term which is due to the variation in tension in the cable during oscillations. For the S1 solution, where the first in-plane mode is non-resonant, this term may loosely be thought of as a non-resonant 'forcing' of the in-plane mode giving rise to twice- and zero-frequency response components—being non-resonant components they appear in $h^{(1)}$, (4.4), as $\beta_{11}(u_1^2 + \bar{u}_1^2 - 6u_1\bar{u}_1)/(3m\omega_{r1}^2)$.

We now assess whether non-resonant response terms can occur in higher modes of the cable due to the resonant response of the first out-of- and in-plane modes. For the nth out-of-plane mode, all the terms in the equation of motion (4.2) are multiples of $y_n$, hence there can only be a response if the mode is resonant. For the in-plane modes, a response is possible due to the $\beta_{n1}(y_1^2 + z_1^2)/m$ terms from the summation $\sum_k(\beta_{kn}/m)(y_k^2 + z_k^2)$. These terms can result in a non-zero $h^{(1)}$, even with zero resonant response in the nth in-plane mode, that is captured by

$$ \left. \begin{array}{l} \text{non-resonant nth in-plane mode} \\ \text{nth} \end{array} \right\} \mathbf{h}_{\text{non-res}}^{(1)} = \frac{\beta_{n1}}{3m\omega_{r1}^2} [u_1^2 + \bar{u}_1^2 + u_2^2 + \bar{u}_2^2 - 6(u_1\bar{u}_1 + u_2\bar{u}_2)], \quad (4.11) $$

following the same approach as was used to calculate (4.4) and then setting $u_n = 0$. Applying the resonant response solution for the first in- and out-of-plane modes using (3.20), noting that $\beta_{1n} = \beta_{11}/n$ for odd n and zero for even n, and applying the phase condition that $\phi_2 - \phi_1 = 90^\circ$, and $\phi_1 = 0$, gives the response for the nth in-plane mode

$$ z_n = \begin{cases} -\frac{\beta_{11}}{2n(4-n^2)m\omega_{r1}^2} U_1^2 \cos(2\omega_{r1}t) - \frac{\beta_{11}}{2n^3m\omega_{r1}^2} U_1^2 & \text{odd } n \text{ with } n > 1 \\ 0 & \text{even } n \end{cases}. \quad (4.12) $$

Here, we have noted that, when $n > 1$, $z_n = h_{\text{non-res}}^{(1)}$ if only the resonant response of the first-in- and out-of-plane modes are considered.

### (d) Nonlinear normal modes

These modal contributions can be added to those from the first-in- and out-of-plane modes, (4.10). Using (4.1), along with the mode shapes, the NNMs can be identified in terms of the lateral displacement $v(x, y)$ and the vertical displacement $w(x, y)$. By approximating the in-plane modes as sinusoidal, an acceptable approximation if the cable is taut [24], the response may be written in terms of amplitude of the resonant responses $U_1$ and $U_2$ as

$$ v(x,t) = \sin\left(\frac{\pi x}{\ell}\right) \left[ U_1 \cos(\omega_{r1}t) + \frac{\nu_{11}}{32m\omega_{r1}^2} U_1(U_1^2 - U_2^2) \cos(3\omega_{r1}t) + \frac{\beta_{11}}{3m\omega_{r1}^2} U_1 U_2 \sin(2\omega_{r1}t) \right], \quad (4.13) $$
in the out-of-plane direction and

$$ w(x, t) = w_s(x) + \sin\left(\frac{\pi x}{\ell}\right) \times \left[ \begin{aligned} & U_2 \sin(\omega_{r1} t) + \frac{\nu_{11}}{32m\omega_{r1}^2} U_2(U_1^2 - U_2^2) \sin(3\omega_{r1} t) - \frac{\beta_{11}}{3m\omega_{r1}^2} U_2^2 \cos(2\omega_{r1} t) - \frac{\beta_{11}}{m\omega_{r1}^2} U_2^2 \\ & - \frac{\beta_{11}}{2m\omega_{r1}^2} (U_1^2 - U_2^2) \cos(2\omega_{r1} t) \sum_{\tilde{n}=1}^{\infty} \frac{1}{(2\tilde{n}-1)(4-(2\tilde{n}-1)^2)} \sin\left((2\tilde{n}-1)\frac{\pi x}{\ell}\right) \\ & - \frac{\beta_{11}}{2m\omega_{r1}^2} (U_1^2 + U_2^2) \sum_{\tilde{n}=1}^{\infty} \frac{1}{(2\tilde{n}-1)^3} \sin\left((2\tilde{n}-1)\frac{\pi x}{\ell}\right) \end{aligned} \right], \quad (4.14) $$

in the in-plane direction. These expressions along with the response frequency and resonant amplitude relationships, (4.7), (4.8) and (4.9), define the NNMs.

In the out-of-plane direction, the NNM mode shape is straightforward, a half sine wave matching the linear mode shape; however, the response contains not only the response frequency but also harmonics. In the in-plane direction, the NNM is more complex with modal components from all the odd in-plane modes. This expression can be significantly simplified using the Fourier series expansions

$$ \begin{aligned} & \sum_{\tilde{n}=1}^{\infty} \frac{1}{(2\tilde{n}-1)(4-(2\tilde{n}-1)^2)} \sin\left((2\tilde{n}-1)\frac{\pi x}{\ell}\right) = \frac{\pi}{8} \sin^2\left(\frac{\pi x}{\ell}\right), \\ & \sum_{\tilde{n}=1}^{\infty} \frac{1}{(2\tilde{n}-1)^3} \sin\left((2\tilde{n}-1)\frac{\pi x}{\ell}\right) = \frac{\pi^3}{8\ell^2}(x\ell - x^2), \end{aligned} $$

such that

$$ w(x, t) = w_s(x) + \sin\left(\frac{\pi x}{\ell}\right) \times \left[ \begin{aligned} & U_2 \sin(\omega_{r1} t) + \frac{\nu_{11}}{32m\omega_{r1}^2} U_2(U_1^2 - U_2^2) \sin(3\omega_{r1} t) - \frac{\beta_{11}}{3m\omega_{r1}^2} U_2^2 \cos(2\omega_{r1} t) - \frac{\beta_{11}}{m\omega_{r1}^2} U_2^2 \\ & - \frac{\beta_{11}\pi}{16m\omega_{r1}^2} (U_1^2 - U_2^2) \cos(2\omega_{r1} t) \sin^2\left(\frac{\pi x}{\ell}\right) - \frac{\beta_{11}\pi^3}{16m\omega_{r1}^2} (U_1^2 + U_2^2) \left[\frac{x}{\ell} - \left(\frac{x}{\ell}\right)^2\right] \end{aligned} \right]. \quad (4.15) $$

It is, perhaps, interesting to note that the last term, which is not a function of time, takes the same shape along the cable as the static sag ($w_s(x) \propto (\ell x - x^2)$ from [24]). It can therefore be seen as reducing the sag effect due to the increased average dynamic tension in the cable as it oscillates. The in-plane NNM response, at the level of accuracy of the normal form transformation, therefore consists of the linear mode shape oscillating at three frequencies along with a dynamically corrected static sag and a further term in the shape of the square of the linear mode shape. Due to this final term, the shape of the NNM in the in-plane direction is dependent on the amplitudes of the resonant responses in the first-in- and out-of-plane modes and varies over each oscillation.

To examine these NNMs, first consider the S1 solution, in which $U_2 = 0$. Figure 4 shows the trajectory of the mid-span deflections $v(\ell/2, t)$ and $w(\ell/2, t)$ for a range of values of $U_1$ up to the bifurcation amplitude in the $y$ versus $-z$ (to account for the fact that $z$ is defined positive downwards) plane using (4.13) and (4.15). For each trajectory, the resonant frequency is different, as shown in figure 3. The dashed line represents the maxima of the trajectories and corresponds to the initial deflection points from which, provided the initial shape is imposed, the system can be released into the pure S1 NNM. This initial shape is a half sine wave out-of-plane regardless of the amplitude of $U_1$ but is a more complex weakly amplitude-dependent shape in-plane. Note that the in-plane deflection plotted here includes static sag—the value of which corresponds to the minimum of the dashed line where the initial out-of-plane deflection is zero. Any deflection along the NNM from this point results in the cable oscillating above the static sag position. For the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Trajectory of mid-span deflections for a range of $U_1$ values on the S1 solution up to the bifurcation amplitude (solid lines) and the locus of points from which a stationary cable can be released onto a pure S1 solution provided the deflected shape along the length of the cable is appropriate (dashed line).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Trajectory of the cable for the S3 solution (a) at the bifurcation point and for a $U_1$ amplitude (b) 0.1% and (c) 1% greater than that at which the bifurcation occurs. Vertical deflections are plotted relative to the static sag shape. Green lines show cable locations at quarter time period points, magenta lines show the eighth-span traces of the cable for all time, the red line show the sag shape and the blue line the chord between the two supports.</figcaption></figure>

S2 solution, as $U_1 = 0$, the motion is confined to the vertical plane such that both the trajectories and the initial deflection points would lie on the vertical axis in the figure.

Considering the S3± solution, figure 5 shows the full cable response in space when $U_1$ is (a) at the bifurcation point and (b) at 0.1% greater than and (c) 1% greater than the amplitude at which the bifurcation occurs. Here, the vertical displacements are plotted relative to the static sag
and hence the sag equilibrium is the straight (red) line between the supports. The cable location is shown at quarter time-period spacing as thick (green) curves and the trace of the cable over all time at eighth-span locations are shown as thin (magenta) loops. The initial transition into whirling, panel (b), is a solution that closely envelops the pre-whirling response, panel (a), as is suggested by the responses (a) and (b) in figure 3. Figure 5a,b also shows that the motion over all time occurs above the sag location. This is perhaps contrary to the expectation that whirling is initiated when the vertical momentum becomes sufficient for the cable to wrap around the zero deflection point rather like a pendulum going 'over the top' as it undergoes a transition from swinging to rotating motion. Instead, the initial whirling motion, panel (b), envelops the pre-whirling solution, and then rapidly becomes a near-elliptical orbit that encircles the sag equilibria as can be seen in panel (c).

## 5. Concluding remarks

The purpose of this paper has been to review the method of normal forms in the context of second-order mechanical vibration problems, to demonstrate its power, and to show how recent developments fit with the historical account of the subject. To this end, the main novel contribution is the re-derivation of the second-order normal form technique using the general theory (§3). This shows that, provided forcing and damping are considered to be weak, then forced-damped systems can be treated in a mathematically consistent way. Moreover, this treatment provides insight into the region of validity of the method. Frequency detuning, forcing amplitude and damping coefficients can be treated as *unfolding* terms of the same order of magnitude as the nonlinear terms. While this assumption of weak forcing and damping may appear restrictive, for smooth nonlinear systems close to resonance, one is usually dealing with internal forces generated from structural stiffnesses that are much larger than excitation and damping forces.

Other novelties of this work include a re-derivation of the homological equation (3.11) without reference to harmonic balance or an artificial book-keeping parameter $\varepsilon$ as was used in Neild & Wagg [22]. Also, we have applied the method to the computation of the NNMs of a multiple-degree-of-freedom model of a cable close to the resonance between in-phase and out-of-phase fundamental modes. This has enabled a general expression for the responses, (4.13) and (4.15), including contributions from arbitrarily many higher order modes to be derived. These fully parametrized expressions have enabled us to draw novel insight into the mechanism of transition from planar to whirling solutions in the problem.

As discussed elsewhere, the second-order normal form method keeps the transform physically relevant. It could be argued that the concept of the normal form is the natural generalization of linear modal analysis for finite degree of freedom systems. Linear modal analysis essentially identifies eigenvectors and their associated natural frequencies and damping constants. In addition, the normal form identifies the key nonlinear terms that must be included for a complete description of the dynamics. Normal form analysis can be compared to other methods for providing approximate solutions to weakly nonlinear vibration problems. There are a wide class of perturbation methods; averaging, harmonic balance, multiple scales, etc. (e.g. [31] or [32]). These have been applied successfully by many authors to produce approximate solutions close to many different types of resonance. While the methods often give the same algebraic solutions, for example for the prediction of backbone curves in the absence of forcing and damping, it is often hard to predict the region of validity of the analysis, and they tend not to explain the full dynamics including harmonics. The normal form is philosophically different; the method works not by looking for specific kinds of solution, but by providing sequences of approximations to the equations of motion themselves. Thus, one still has access to the complete dynamics. Moreover, the method is algorithmic and can readily be carried out computationally.

It is also useful to compare normal form analysis with direct calculation of NNMs. Our application to cable dynamics in §4 has shown that the normal form provides a natural way to compute local (i.e. weak amplitude) approximations to NNMs in undamped systems.
However, as argued in the caption to figure 1, NNMs are not uniquely defined in the presence of damping, which provides both conceptual and computational problems. For example, in a recent paper by Renson *et al.* [33], a finite-element method was proposed for computing NNMs directly. In example 5.2 of that paper, a large-amplitude approximation to a NNM is computed under two different assumptions on the magnitude of damping for a two-degree-of-freedom nonlinear oscillator. Not only is it found that the NNM shape depends strongly on the damping coefficient, but the NNM is only computed for the most strongly damped mode; it is not uniquely defined for the other mode. Unfortunately, in reality, the least damped mode is the one that is most likely to be excited in practice.

By contrast, the normal form method is not limited by restriction to the most heavily damped mode. Indeed, normal form methods are routinely used in conjunction with Melnikov's method where the forcing and damping terms play the role of unfolding parameters. This enables the complete dynamics (not just that restricted to putative NNMs) to be studied in the neighbourhood of a resonance (e.g. [34]). Furthermore, in this paper, it has been shown that damping and forcing can naturally be included as unfolding parameters in the second-order normal form. Examples of the second-order normal form method applied to forced and damped vibration problems were presented in Wagg & Neild [24].

Finally, we note that the second-order normal form method is naturally generalizable to higher degrees of freedom, and to cases of multi-frequency excitation, parametric resonance and multiple or internal resonances. These will simply result in different detuning parameters $\delta$ in the notation of §3a. An exploration of such cases is left for future work, as is application of the method to systems of higher degrees of freedom using computational methods.

**Authors' contributions.** A.R.C., S.A.N. and D.J.W. working on recasting the second-order technique in the general theory of normal forms; S.A.N., T.L.H. and A.C. developed the cable application example. All authors helped to draft the manuscript and gave final approval for publication.

**Competing interests.** We declare we have no competing interests.

**Funding.** The authors acknowledge the support of the Engineering and Physical Sciences Research Council. S.A.N. and T.L.H. are supported by EP/K005375/1, and A.R.C., D.J.W. and A.C. are supported by EP/K003836/1.

## References

1. Kerschen G, Peeters M, Golinval JC, Vakakis AF. 2009 Nonlinear normal modes, Part I: a useful framework for the structural dynamicist. *Mech. Syst. Signal. Process.* **23**, 170–194. (doi:10.1016/j.ymssp.2008.04.002)

2. Peeters M, Viguié R, Sérandour G, Kerschen G, Golinval JC. 2009 Nonlinear normal modes, Part II: toward a practical computation using numerical continuation techniques. *Mech. Syst. Signal. Process.* **23**, 195–216. (doi:10.1016/j.ymssp.2008.04.003)

3. Avramov KV, Mikhlin YV. 2013 Review of applications of nonlinear normal modes for vibrating mechanical systems. *Appl. Mech. Rev.* **65**, 020801. (doi:10.1115/1.4023533)

4. Guckenheimer J, Holmes PJ. 1983 *Nonlinear oscillations, dynamical systems, and bifurcations of vector fields*. New York, NY: Springer.

5. Chow SN, Li C, Wang D. 1994 *Normal forms and bifurcations of planar vector fields*. Cambridge, UK: Cambridge University Press.

6. Murdock J. 2002 *Normal forms and unfoldings for local dynamical systems*. New York, NY: Springer.

7. Jezequel L, Lamarque CH. 1991 Analysis of non-linear dynamic systems by the normal form theory. *J. Sound Vib.* **149**, 429–459. (doi:10.1016/0022-460X(91)90446-Q)

8. Neild SA, Wagg DJ. 2011 Applying the method of normal forms to second order nonlinear vibration problems. *Proc. R. Soc. A* **467**, 1141–1163. (doi:10.1098/rspa.2010.0270)

9. Shaw S, Pierre C. 1991 Non-linear normal modes and invariant manifolds. *J. Sound Vib.* **150**, 170–173. (doi:10.1016/0022-460X(91)90412-D)

10. Shaw SW, Pierre C. 1994 Normal modes of vibration for non-linear continuous systems. *J. Sound Vib.* **169**, 319–347. (doi:10.1006/jsvi.1994.1021)
11. Montaldi J, Roberts M, Stewart I. 1990 Existence of nonlinear normal modes of symmetric Hamiltonian systems. *Nonlinearity* **3**, 695–730. (doi:10.1088/0951-7715/3/3/009)

12. Nayfeh AH. 1993 *Method of normal forms*. New York, NY: Wiley.

13. Poincaré H. 1892–9 New methods of celestial mechanics. New York, NY: American Institute of Physics. [English transl. edited by Goroff DL 1993.]

14. Barrow-Green J. 1996 *Poincaré and the three body problem*. History of Mathematics, vol. 11. Providence, RI: American Mathematical Society.

15. Diacu F, Holmes PJ. 1999 *Celestial encounters*. Princeton, NJ: Princeton University Press.

16. Verhulst F. 2012 *Henri Poincaré: impatient genius*. New York, NY: Springer.

17. Meyer KR, Hall GR, Offin D. 2009 *Introduction to Hamiltonian dynamical systems and the N-body problem*, 2nd edn. New York, NY: Springer.

18. Birkhoff GD. 1927 *Dynamical systems*, vol. IX. Providence, RI: American Mathematical Society Colloquium Publication.

19. Kuznetsov YuA. 2004 *Elements of applied bifurcation theory*. New York, NY: Springer.

20. Champneys AR, Kirk V. 2004 The entwined wiggling of homoclinic curves emerging from saddle-node/Hopf instabilities. *Phys. D* **195**, 77–105. (doi:10.1016/j.physd.2004.03.004)

21. Arnol'd VI. 1988 *Geometrical methods in the theory of ordinary differential equations*. New York, NY: Springer.

22. Neild SA. 2012 Approximate methods for analysing nonlinear structures. In *Exploiting nonlinear behaviour in structural dynamics* (eds LN Virgin, DJ Wagg), pp. 77–84. New York, NY: Springer.

23. Neild SA, Wagg DJ. 2013 A generalized frequency detuning method for multidegree-of-freedom oscillators with nonlinear stiffness. *Nonlinear Dyn.* **73**, 649–663. (doi:10.1007/s11071-013-0818-7)

24. Wagg DJ, Neild SA. 2014 *Nonlinear vibration with control: for flexible and adaptive structures*, 2nd edn. New York, NY: Springer.

25. Hill TL, Cammarano A, Neild SA, Wagg DJ. 2015 Interpreting the forced responses of a two-degree-of-freedom nonlinear oscillator using backbone curves. *J. Sound Vib.* **349**, 276–288. (doi:10.1016/j.jsv.2015.03.030)

26. Doedel E, Oldeman BE, Champneys AR, Dercole F, Fairgrieve T, Kuznetsov Yu, Wang X, Zhang C 2007 AUTO-07P: Continuation and bifurcation software for ordinary differential equations. See http://sourceforge.net/projects/auto-07p/.

27. Hill TL, Cammarano A, Neild SA, Wagg DJ. 2015 Out-of-unison resonance in weakly nonlinear coupled oscillators. *Proc. R. Soc. A* **471**, 20140659. (doi:10.1098/rspa.2014.0659)

28. Warnitchai Y, Fujino T, Susumpov A. 1995 A nonlinear dynamic model for cables and its application to a cable structure-system. *J. Sound Vib.* **187**, 695–712. (doi:10.1006/jsvi.1995.0553)

29. Gonzalez-Buelga A, Neild SA, Wagg DJ, Macdonald JHG. 2008 Modal stability of inclined cables subjected to vertical support excitation. *J. Sound Vib.* **318**, 565–579. (doi:10.1016/j.jsv.2008.04.031)

30. Macdonald JHG, Dietz MS, Neild SA, Gonzalez-Buelga A, Crewe AJ, Wagg DJ. 2010 Generalised modal stability of inclined cables subjected to support excitations. *J. Sound Vib.* **329**, 4515–4533. (doi:10.1016/j.jsv.2010.05.002)

31. Nayfeh AH. 2000 *Perturbation methods*. New York, NY: Wiley.

32. Holmes HM 2012 *Introduction to perturbation methods*, 2nd edn. New York, NY: Springer.

33. Renson L, Deliege G, Kerschen G. 2014 An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems. *Meccanica* **49**, 1901–1916. (doi:10.1007/s11012-014-9875-3)

34. Haller G 1999 *Chaos near resonance*. New York, NY: Springer.