How do Conservative Backbone Curves Perturb into Forced Responses? A Melnikov Function Analysis

TL;DR

This paper provides a mathematical framework using Melnikov functions to predict how conservative backbone curves influence forced responses in weakly damped mechanical systems, clarifying when certain periodic orbits dominate the response.

Contribution

It introduces a systematic analysis predicting which conservative periodic orbits form backbone curves in forced responses and establishes conditions for approximate methods and isolas.

Findings

Predicts which conservative orbits serve as backbone curves

Derives conditions for the validity of energy balance and force appropriation methods

Provides criteria for the emergence of isolated response branches (isolas)

Abstract

Weakly damped mechanical systems under small periodic forcing tend to exhibit periodic response in a close vicinity of certain periodic orbits of their conservative limit. Specifically, amplitude frequency plots for the conservative limit have often been noted, both numerically and experimentally, to serve as backbone curves for the near resonance peaks of the forced response. In other cases, such a relationship between the unforced and forced response was not observed. Here we provide a systematic mathematical analysis that predicts which members of conservative periodic orbit families will serve as backbone curves for the forced-damped response. We also obtain mathematical conditions under which approximate numerical and experimental approaches, such as energy balance and force appropriation, are justifiable. Finally, we derive analytic criteria for the birth of isolated response branches (isolas) whose identification is otherwise challenging from numerical continuation.

Equations

$$H(q,\dot{q})=E(q,\dot{q})+V(q)=\frac{1}{2}⟨\dot{q},M(q)\dot{q}⟩+V(q).$$

$$M(q)\ddot{q}+G(q,\dot{q})+DV(q)=\epsilon Q(q,\dot{q},t;T,\epsilon ),$$

$$\dot{x}=f(x)+\epsilon g(x,t;T,\epsilon ),$$

$$\begin{array}[]{lcr}\displaystyle f(x)=\begin{pmatrix}\dot{q}\\ -M^{-1}(q)(DV(q)+G(q,\dot{q}))\end{pmatrix},&&\displaystyle g(x,t;T,\varepsilon)=\begin{pmatrix}0\\ M^{-1}(q)Q(q,\dot{q},t;T,\varepsilon)\end{pmatrix}.\end{array}$$

$$\dot{Y}=Df(x_0(t;p))Y,Y(0)=I.$$

$${\Delta }_l:{\mathbb{R}}^{n+2}\to {\mathbb{R}}^n,{\Delta }_l(\xi ,T,\epsilon )=x(lT;\xi ,T,\epsilon )-\xi ,{\Delta }_l\in C^r,$$

$$\xi =x_0(s;p)+O(\epsilon ),p\in \mathcal{Z},T=\tau m/l+O(\epsilon ),$$

$$L(\xi ,lT)=L(p,m\tau ).$$

$${\Delta }_{l,L}:{\mathbb{R}}^{n+2}\to {\mathbb{R}}^{n+1},{\Delta }_{l,L}(\xi ,T,\epsilon )=\{\begin{array}{l}{\Delta }_l(\xi ,T,\epsilon )\\ L(\xi ,lT)-\lambda \end{array},{\Delta }_{l,L}(\xi ,T,\epsilon )=0.$$

$$\displaystyle M^{m:l}(s)=\int_{0}^{m\tau}\big{\langle}\,DH(x_{0}(t+s;p))\,,\,g(x_{0}(t+s;p),t;\tau m/l,0)\,\big{\rangle}\,dt,$$

$$\begin{array}[]{lcr}M^{m:l}(s_{0})=0,&&DM^{m:l}(s_{0})\neq 0,\end{array}$$

$$\begin{array}[]{lcr}\displaystyle M^{m:l}(s_{0},\kappa_{0})=D_{s}M^{m:l}(s_{0},\kappa_{0})=0,&&\displaystyle D^{2}_{ss}M^{m:l}(s_{0},\kappa_{0})\neq 0.\end{array}$$

$$E_{b}(\xi,\varepsilon)_{[0,lT]}=\varepsilon\int_{0}^{lT}\big{\langle}\,\dot{q}(t;\xi,T,\varepsilon)\,,\,Q(q(t;\xi,T,\varepsilon),\dot{q}(t;\xi,T,\varepsilon),t;T,\varepsilon)\big{\rangle}dt,$$

$$E_b(\xi ,\epsilon {)}_{[0,lT]}=0,$$

$$E_b(\xi ,\epsilon {)}_{[0,lT]}=\epsilon M^{m:l}(s)+O({\epsilon }^2),$$

$$M^{m:l}(s)=\int_{0}^{m\tau}\big{\langle}\,\dot{q}_{0}(t+s;p)\,,\,Q(q_{0}(t+s;p),\dot{q}_{0}(t+s;p),t;\tau m/l,0)\big{\rangle}dt,$$

$$DH(x)=DH(q,\dot{q})=\left(\begin{array}{cc}DV(q)+D_{q}E(q,\dot{q}),& M(q)\dot{q}\end{array}\right).$$

$$\begin{array}[]{lcr}Q(q,\dot{q},t;T,e,0)=ef_{e}\cos(\Omega t)-C(q,\dot{q}),&&\Omega=2\pi/T,\end{array}$$

$$\begin{array}[]{rl}M^{m:l}(s,e)=&\displaystyle e\int_{0}^{m\tau}\langle\dot{q}_{0}(t+s;p),f_{e}\rangle\cos\left(\frac{2l\pi}{m\tau}t\right)dt+\\ \\ &\displaystyle-\int_{0}^{m\tau}\langle\dot{q}_{0}(t+s;p),C(q_{0}(t+s;p),\dot{q}_{0}(t+s;p))\rangle dt\\ \\ =&\displaystyle w^{m:l}(s,e)-mR,\end{array}$$

$$R={\int }_0^{\tau }⟨{\dot{q}}_0(t;p),C(q_0(t;p),{\dot{q}}_0(t;p))⟩dt,$$

$$w^{m:l}(s,e)=e{\int }_0^{m\tau }⟨{\dot{q}}_0(t+s;p),f_e⟩\cos \left(\frac{2l\pi }{m\tau }t\right)dt$$

$$\begin{array}[]{lcr}\displaystyle q_{0}(t;p)=\frac{a_{0}}{2}+\sum_{k=1}^{\infty}a_{k}\cos\left(k\omega t\right)+b_{k}\sin\left(k\omega t\right),&&\displaystyle\omega=\frac{2\pi}{\tau},\end{array}$$

$$w^{m:l}(s,e)=\{\begin{array}{ll}0& \mathrm{if}m\ne 1\\ W^{1:l}(e)\cos \left(l\omega s-{\alpha }_{l,e}\right)& \mathrm{if}m=1\end{array},$$

$$\begin{array}[]{lcccr}W^{1:l}(e)=eA_{l,e},&&A_{l,e}=l\pi\sqrt{\displaystyle\left\langle a_{l},f_{e}\right\rangle^{2}+\left\langle b_{l},f_{e}\right\rangle^{2}},&&\displaystyle\alpha_{l,e}=\arctan\left(\frac{\langle a_{l},f_{e}\rangle}{\langle b_{l},f_{e}\rangle}\right).\end{array}$$

$$M^{1:l}(s,e)=W^{1:l}(e)\cos \left(l\omega s-{\alpha }_{l,e}\right)-R.$$

$$|D_e{M}^{1:l}(s_{qz},e_0)|=|D_e{W}^{1:l}(e_0)|=A_{l,e}>0$$

$${\Gamma }_l(\lambda )=R(\lambda )/A_{l,e}(\lambda ).$$

$$M^{1:l}(s,e,\lambda)=A_{l,e}(\lambda)\Big{(}e\cos\big{(}l\omega(\lambda)s-\alpha_{l,e}(\lambda)\big{)}-\Gamma_{l}(\lambda)\Big{)},$$

$$D_{\lambda}M^{1:l}\big{(}s_{qz}(\lambda_{0}),\Gamma_{l}(\lambda_{0}),\lambda_{0}\big{)}=-A_{l,e}(\lambda_{0})D\Gamma_{l}(\lambda_{0})\neq 0,$$

$$W^{1:l}(e_0,{\lambda }_1)=e_0{A}_{l,e}({\lambda }_1)={\Gamma }_l({\lambda }_0)A_{l,e}({\lambda }_1)>{\Gamma }_l({\lambda }_1)A_{l,e}({\lambda }_1)=R({\lambda }_1),$$

Full text

1 Introduction

Conservative families of periodic orbits, broadly known as nonlinear normal modes (NNMs) in the field of structural engineering, often appear to shape the behaviour of mechanical systems even in the presence of additional damping and time-periodic forcing. Not only has this influence been noted for low-amplitude vibrations of a small number of coupled oscillators, but it also appears to hold for large amplitude motion in arbitrary degrees of freedom. Various descriptions of this phenomenon are available in the literature, ranging from the original introduction of NNMs by Rosenberg [1] to the more recent reviews lead by Vakakis [2, 3, 4], Avramov and Mikhlin [5, 6] and Kerschen [7].

These studies suggest that forced-damped frequency responses might bifurcate from conservative NNMs. To summarise features of such bifurcations, Fig. 1 qualitatively represents possible nonlinear phenomena in the frequency response. Dark and light red curves show steady-state solutions for a periodically forced-damped mechanical system for low and high forcing amplitudes, respectively. Blue curves correspond to amplitude-frequency relations of periodic orbit families of the underlying conservative system, which we refer to as conservative backbone curves. Grey curves depict backbone curves of the forced-damped response, defined as the frequency locations of amplitude maxima under variation of the forcing amplitude [8]. The first and last peaks in Fig. 1 show the classic hardening and softening resonance trends respectively. As most frequently observed behaviours, these two phenomena have been broadly studied: see, e.g. [8, 7] for analytical and numerical treatments and [9, 10] for experimental results. In these settings, as the response amplitude increases, the backbone curves of the conservative limit and those of the forced-damped system have been noted to pull apart [11].

The backbone curves of the second peak from the left in Fig. 1 feature a non-monotonic trend in frequency, i.e., softening for lower amplitudes and hardening for higher ones. This type of behaviour is relevant for shallow-arch systems [12], MEMS devices [13] and structural elements [14], while a simple mechanical example is available at [15, 16]. The third peak of Fig. 1 shows an isolated response curve (isola) that may occur due to the influence of nonlinear damping [17] or symmetry breaking mechanisms [18]. In the former case, the isola can join with the main branch when the forcing amplitude exceeds a certain threshold [19], as indicated by the red light curve. Subharmonic responses, displayed along the second blue line from the right, also show up as isolated branches, as they cannot originate from the linear limit [8]. Finally, as highlighted in [20] with the analysis of a simple model of nonlinear beam, not all NNMs contribute to shaping the forced-response, as illustrated by the rightmost backbone curve in Fig. 1.

Analytic relations between conservative families of periodic orbits and frequency responses are only available for specific, low-dimensional oscillators from perturbation expansions that assume the conservative periodic limit to have small amplitudes. These expansions may arise from the method of multiple scales [8], averaging [21], the first-order normal form technique [22] or the second-order normal form technique [23]. Based on this latter method, Hill et al. [24] developed an energy-transfer formulation for locating maxima of the frequency response along conservative backbone curves. However, the authors a priori postulate a relation between conservative oscillations and frequency responses, and also discuss potential limitations arising from this assumption in [24]. Vakakis and Blanchard in [25] show the exact steady states of a strongly nonlinear periodically forced and damped Duffing oscillator. They also clarify how these forced-damped periodic motions are related the conservative backbone curve, but they restrict the discussion to specific types of periodic forcing.

Relying on exact mathematical results, spectral submanifold (SSM) theory [26, 27] has been developed for the local analysis of damped, nonlinear oscillators. This approach is insensitive to the number of degrees of freedom thanks to the automated procedure developed in [28] and can be hence used for exact nonlinear model reduction. SSMs, however, do not exist in the limit of zero forcing and damping, in which case they are replaced by Lyapunov subcenter manifolds (LSMs) [29]. The relationship between the dynamics on damped, unforced SSMs and LSMs is now established in a small enough neighbourhood of the unforced equilibrium point [27, 19].

Even though diverse numerical options are available to explore forced responses, analytical tools applicable to arbitrary degrees of freedom and motion amplitudes are still particularly valuable. Not only can such tools help with the analysis of several perturbation types by relying only on the knowledge of conservative orbits, but they can also overcome limitations of numerical routines. For a thorough review of these limitations, we refer the reader to [30]. For example, direct numerical integration needs long computational time for high degrees of freedom systems with small damping and it is limited to stable periodic orbits. Despite being very accurate, numerical continuation (shooting methods [31], harmonic balance [32] or collocation [33]) suffers from the curse of dimensionality. Furthermore, it can efficiently compute the main branch of the forced response for frequencies away from resonance, but fails to find isolated branches when their existence and location is not a priori known.

In a parallel development, SSMs, backbone curves, main and isolated branches can now be computed efficiently for general, forced-damped, multi-degree-of-freedom mechanical systems up to any required order of accuracy [19]. This approach also yields analytic approximations for backbone curves and isolas, as long as these stay within the domain of convergence of Taylor expansions constructed for the forced-damped SSMs [19].

Beyond these numerical approaches, experimental methods would also be aided by a rigorous mathematical relation between conservative backbone curves and forced-damped responses. One of such methods has been developed by Peeters et al. [34, 9], who propose an extension of the phase-lag quadrature criterion, well-known for linear systems, to nonlinear systems in order to isolate conservative NNMs experimentally. Their method assumes that the nonlinear periodic motion is synchronous [1], the damping is linear (or at least odd in the velocity) and the excitation is multi-harmonic and distributed. The idea is to subject the system to forcing that exactly balances damping along a periodic orbit of the conservative limit. Peeters et al. find that such a balance holds approximately when each harmonic of the conservative periodic orbit has a phase lag of $90^{\circ}$ relative to the corresponding forcing one. Force appropriation [35, 36] and control-based continuation [37] exploit the phase-lag quadrature criterion for systematic tracking of backbone curve.

Another related experimental method in need of a mathematical justification is resonance decay [34], which uses a force appropriation routine to isolate a periodic orbit, then turns off the forcing and assumes the response to converge to the equilibrium position along a two-dimensional SSM, sometimes called a damped-NNM [38]. This technique is expected to provide an approximation of conservative backbone curves, but it remains partially unjustified for two reasons. On the analytical side, it assumes a yet unproven purely, parasitic effect of damping on the response. On the experimental side, decoupling the shaker from the system remains a challenge that affects the accuracy of the results.

Despite available experimental and numerical observations, it is unclear if and when conservative NNMs perturb into forced-damped periodic responses. This is because periodic orbits in conservative systems are never structurally stable under generic perturbations, which tend to destroy them [39]. Indeed, any conservative periodic orbit has at least two Floquet multipliers equal to $+1$ due to the conservation of energy [40], rendering the orbit structurally unstable. Classic analytic approaches [41, 42, 43] to generic, non-autonomous perturbations of normally hyperbolic periodic orbits are therefore inapplicable in this setting. Conservative NNMs exist in families that are only guaranteed to persist under small, smooth conservative perturbations [40, 44].

In its simplest form, the study of dissipative perturbations on a conservative family of periodic orbits dates back to Poincaré [45], developed further by Arnold [46]. An important contribution was made by Melnikov [47], who focused on dissipative perturbations of planar Hamiltonian systems. His approach reduces the persistence problem of periodic orbits to the analysis of the zeros of the subharmonic Melnikov function [39, 48]. As extensions of Melnikov’s approach, studies on two-degree-of-freedom Hamiltonian systems are available: [49, 50] consider the fate of periodic orbit families in an integrable system subject to Hamiltonian perturbations, while [51] analyses two fully decoupled oscillators under generic dissipative perturbations. Subharmonic Melnikov-type analysis for non-smooth systems is also available; see [52] for examples of planar oscillators and [53, 54] for a system with two degrees of freedom. All these results, therefore, require low-dimensionality and integrability before perturbation, neither of which is the case for the conservative limits of nonlinear structural vibrations problems arising in practice.

As an alternative to these analytic methods, Chicone [55, 56, 57] established a perturbation method for manifolds of isochronous periodic orbits without any restriction on their Floquet multipliers or assumptions on integrability/coupling before perturbation. This elegant approach exploits the Lyapunov-Schmidt reduction to obtain a generic multi-dimensional bifurcation function for the persistence of single orbits. Furthermore, this method has also been extended for non-smooth (but Lipschitz) dynamical systems [58]. However, these results are not directly applicable to the typical setting of nonlinear structural vibrations. Moreover, an exact resonance condition is required in Chicone’s method, even though perturbed periodic orbits are often observed when the forcing frequency clocks in near-resonance with the frequency of a periodic orbit of the conservative limit.

In this paper, we develop an exact analytical criterion for the perturbation of conservative NNMs into forced-damped periodic responses, thereby predicting the variety of behaviours depicted in Fig. 1. We assume that the conservative limit of the system has a one-parameter family of periodic orbits, satisfying generic nondegeneracy conditions. We then study the persistence or bifurcation of these periodic orbits under small damping and time-periodic forcing. Utilising ideas from Rhouma and Chicone [57], we reduce this perturbation problem to the study of the zeros of a Melnikov-type function, generalising therefore the original Melnikov method to multi-degree-of-freedom systems. Our approach relies on the smallness of dissipative and forcing terms which is generally satisfied in structural dynamics, but we will not assume that the unperturbed periodic orbit has small amplitude. This distinguishes our approach from various classic perturbation expansions that assume closeness from the unforced equilibrium. Our analysis also differ from classic Melnikov-type approaches in that it does not require the conservative limit of the system to be integrable. Indeed, for our unperturbed conservative limit, we only require the existence of a generic family of periodic orbits that may only be known from numerical continuation.

When our Melnikov-type method is applied to mechanical systems, it provides a rigorous justification for the classic energy principle, a broadly used but heuristic necessary condition for the existence of periodic response in forced-damped nonlinear oscillations [18, 24, 34]. Our analysis shows that under further conditions, the energy principle becomes a rigorous sufficient condition for nonlinear periodic response and extends to arbitrary number of degrees of freedom, multi-harmonic forcing, large-amplitude periodic orbits and higher-order external resonances.

We first give a mathematical formulation for general dynamical systems, then apply it to the classic setting of nonlinear structural vibrations. In that context, our results reveal how the near-resonance part of the main and isolated branches of the periodic response diagram are born out of the conservative backbone curve. We also discuss how our results justify the phase-lag quadrature criterion under more general conditions than prior studies assumed. Finally, we illustrate the power of our analytic predictions on a six-degree-of-freedom mechanical system.

2 Setup

We consider a mechanical system with $N$ degrees of freedom and denote its generalised coordinates by $q\in U\subset\mathbb{R}^{N}$, $N\geq 1$. We assume that this system is a small perturbation of a conservative limit that conserves the total energy $H:U\times\mathbb{R}^{N}\rightarrow\mathbb{R}$, defined as

[TABLE]

Here, $M(q)$ is the positive definite, symmetric mass matrix, $E(q,\dot{q})$ is the kinetic energy and $V(q)$ the potential. The equations of motion for the system take the form

[TABLE]

where $\varepsilon\geq 0$ is the perturbation parameter, $G(q,\dot{q})=D_{t}(M(q))\dot{q}-D_{q}E(q,\dot{q})$ contains inertial forces and $\varepsilon Q(q,\dot{q},t;T,\varepsilon)=\varepsilon Q(q,\dot{q},t+T;T,\varepsilon)$ denotes a small perturbation of time-period $T$. System in Eq. 2 is then a weakly non-conservative mechanical system.

Introducing the notation $x=(q,\dot{q})\in\mathbb{R}^{n}$ with $n=2N$, the equivalent first-order form reads

[TABLE]

where we assume that $f\in C^{r}$ with $r\geq 2$, while $g$ is $C^{r-1}$ in $t$ and $C^{r}$ with respect to the other arguments. These vector fields are defined as

[TABLE]

We assume any further parameter dependence in our upcoming derivations to be of class $C^{r}$. Trajectories of (3) that start from $\xi\in\mathbb{R}^{n}$ at $t=0$ will be denoted with $x(t;\xi,T,\varepsilon)=(q(t;\xi,T,\varepsilon),\dot{q}(t;\xi,T,\varepsilon))$. We will also use the shorthand notation $x_{0}(t;\xi)=(q_{0}(t;\xi,),\dot{q}_{0}(t;\xi,))=x(t;\xi,T,0)$ for trajectories of the unperturbed (conservative) limit of system (3). We recall that, for $\varepsilon=0$, energy is conserved, i.e., $H(x_{0}(t;\xi))=H(\xi)$ holds as long as $x_{0}(t;\xi)\in U$.

3 Non-autonomous resonant perturbation of normal families of conservative periodic orbits

In this section, we first state our main mathematical results for single conservative orbits then for orbit families. We also discuss the physical meaning of these results in the context of mechanical systems.

For the $\varepsilon=0$ limit of system (3), we assume that there exists a periodic orbit $\mathcal{Z}\subset U$ of minimal period $\tau>0$ and we denote by $\Pi(p)$ the monodromy matrix111The monodromy matrix, or linearised period-$\tau$ mapping, $\Pi(p):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is defined as $\Pi(p)=Y(\tau;p)$ where $Y$ solves the equation of variations [59]

$\dot{Y}=Df(x_{0}(t;p))Y,\,\,\,\,\,\,\,\,\,\,\,Y(0)=I.$

based at any point $p\in\mathcal{Z}$. We consider $m\in\mathbb{N}^{+}$ multiples of the period and let $\mu_{a,m}$ denote the algebraic multiplicity of the $+1$ eigenvalue of $\Pi^{m}(p)$ and $\mu_{g,m}$ denote its geometric multiplicity. Note that these two multiplicities are invariant under translations along the orbit, while they may change for different values of $m$. We will need the following definition from [60].

Definition 3.1**.**

A conservative periodic orbit $\mathcal{Z}$ is $m$-normal if one of the following holds:

1. (a)

$\mu_{g,m}=1$* ;* 2. (b)

$\mu_{g,m}=2$* and $f(p)\notin$ range*$(\Pi^{m}(p)-I)$* .*

This normality is a nondegeneracy condition under which a one-parameter family, $\mathcal{P}$, of $m$-normal periodic solutions of the vector field $f$ emanates from $\mathcal{Z}$ (see Theorem 4 of [60] or Theorem 7 of [61]). We denote with $\lambda\in\mathbb{R}$ the parameter identifying each individual orbit in the family $\mathcal{P}$.

Figure 2 describes the types of $m$-normal periodic orbits covered by Definition 3.1, with their associated backbone-curve geometry, as given in Theorem 5 of [60]. The backbone curve can be parametrised as $(\tau(\lambda),h(\lambda))$, with $\tau$ denoting the orbit period and $h$ the value of the first integral. The value of the parameter $\lambda$ is given by a scalar mapping $\lambda=L(\xi,\tau)$ depending on the initial condition $\xi\in\mathbb{R}^{n}$ and the period $\tau\in\mathbb{R}^{+}$. We also require $L$ to be invariant under translations of $\xi$ along the orbit. For an $m$-normal periodic orbit belonging to case (a) in Definition 3.1, one can simply choose $L(\xi,\tau)=\tau$. Instead, when $\mu_{a,m}=2$ (see (i) and (ii) in Fig. 2), the orbit family can be locally parametrised with the value of the first integral $h$, i.e., $L(\xi,\tau)=H(\xi)$. Other possible parametrisations include the value of a coordinate determined by a Poincaré section, the $L^{2}$ norm of the trajectory or the maximum value of a coordinate along the trajectory. For continuation through cusp points of backbone curves, i.e., (iv) in Fig. 2, $L$ may be chosen to provide a pre-defined relation between energy and period (see [60]), but that is outside the scope of this paper.

3.1 Perturbation of a single orbit

Our starting point in the analysis of the fate of perturbed periodic orbits is the displacement map

[TABLE]

whose zeros correspond to $lT$-periodic orbits for system (3) for $l\in\mathbb{N}^{+}$. We aim to smoothly continue normal periodic orbits in the family $\mathcal{P}$ that exists at $\varepsilon=0$. We consider an $m$-normal periodic orbit $\mathcal{Z}\subset\mathcal{P}$ and assume that $l$ and $m$ are relatively prime integers, i.e., $1$ is their only common divisor. We then look for zeros of (5) that can be expressed as

[TABLE]

under the additional constraint

[TABLE]

Equation (7) represents a resonance condition as it relates, either explicitly or implicitly, the periods of the perturbation with that of the periodic orbit $\mathcal{Z}$. With the notation $L(p,m\tau)=\lambda$, the zero problem to be solved reads

[TABLE]

Defining the smooth, $L$-independent, $m\tau$-periodic function $M^{m:l}:\mathbb{R}\rightarrow\mathbb{R}$ as

[TABLE]

we obtain the following main result.

Theorem 3.1**.**

If the Melnikov function $M^{m:l}(s)$ has a simple zero at $s_{0}\in\mathbb{R}$, i.e.,

[TABLE]

then the $m$-normal periodic orbit $\mathcal{Z}$ of the $\varepsilon=0$ limit smoothly persists in system (3) for small $\varepsilon>0$. Moreover, in this case, there exists at least another topologically transverse zero in the interval $(s_{0},\,s_{0}+m\tau)$. If $M^{m:l}(s)$ remains bounded away from zero, then $\mathcal{Z}$ does not smoothly persists for small $\varepsilon>0$.

We prove Theorem 3.1 in Appendix A. The proof reduces the $(n+1)$-dimensional persistence problem to the analysis of the zeros of the scalar function (9). This function formally agrees with the one derived originally by Melnikov for a planar oscillator, but the proof for $n>2$ is more involved compared to the simple geometric treatment in [39] for $n=2$.

When the Melnikov function has a simple or transverse zero, the perturbed orbit emanating from the $m$-normal periodic orbit $\mathcal{Z}$ and its period are $O(\varepsilon)$-close to $\mathcal{Z}$ and to $\tau m$, respectively. Since topologically transverse zeros of functions are generically simple, we expect from Theorem 3.1 that an even number of perturbed periodic orbits bifurcates from the $m$-normal periodic orbit at the $\varepsilon=0$ limit, as indeed typically observed in literature. Moreover, since the Melnikov function (9) does not depend on the parametrisation function $L$ used in Eq. (7), Theorem 3.1 and its consequences hold for any possible parametrising direction used for the unperturbed periodic orbit family.

Theorem 3.1 can be interpreted directly in terms of the backbone curve of $\mathcal{P}$ and the frequency response of system (3). Suppose that the backbone curve of $\mathcal{P}$ shows only regular points near the $m$-normal periodic orbit $\mathcal{Z}$ so that we can select $L(p,m\tau)=m\tau=\lambda$. In this case, Eq. (7) imposes the exact resonance condition $lT=m\tau$. For a pair of simple zeros of $M^{m:l}$, Theorem 3.1 guarantees that the point in the backbone curve corresponding to $\mathcal{Z}$ bifurcates in two frequency responses along the isochronous direction, as depicted in Fig. 3(a). If, instead, $\mathcal{Z}$ corresponds to a fold point in $\tau$, then Theorem 3.1 does not hold for this choice of $L$, but we can still use the energy $h$ as parametrisation variable. In that case, our perturbation method constrains the perturbed initial condition $\xi$ to lie in the same energy level as that of $\mathcal{Z}$. At the same time, the time period $T$ for the perturbed orbit is $O(\varepsilon)$-close to $\tau m/l$, i.e., a near-resonance condition is satisfied. For two simple zeros of $M^{m:l}$, $\mathcal{Z}$ can be smoothly continued in two frequency responses along the isoenergetic direction, as shown in Fig. 3(b) and 3(c).

Remark 3.1**.**

While for the classic planar oscillator case the Melnikov function is also able to predict the stability of perturbed orbits [48], the stability analysis of persisting periodic orbits is more involved for $n>2$. Indeed, their stability depends on all the Floquet multipliers of the conservative limit, as well as on the nature of the perturbation.

3.2 Perturbation of a family and parameter continuation

Here we consider an additional parameter $\kappa\in\mathbb{R}$ in Eq. (8), where $\kappa$ is either a feature of the vector fields in system (3) or the family parameter $\lambda$. The Melnikov function $M^{m:l}$ in (9) clearly inherits this smooth parameter dependence.

Next we investigate the fate of the $m$-normal periodic orbit $\mathcal{Z}$ in the family $\mathcal{P}$ for which the Melnikov function features a quadratic zero at $(s_{0},\kappa_{0})$ defined as:

[TABLE]

The following theorem describes what generic bifurcations may arise in this setting.

Theorem 3.2**.**

Assume that $M^{m:l}(s,\kappa)$ has a quadratic zero at $(s_{0},\kappa_{0})$, as defined in Eq. (11). If $D_{\kappa}M^{m:l}(s_{0},\kappa_{0})\neq 0$, then $\kappa_{sn}=\kappa_{0}+O(\varepsilon)$ is a bifurcation value at which a saddle-node bifurcation of periodic orbits occurs. If $D_{\kappa}M^{m:l}(s_{0},\kappa_{0})=0$ and $\mathrm{det}(D^{2}M^{m:l}(s_{0},\kappa_{0}))>0$ (resp. $<0$), then isola births (resp. simple bifurcations) arise for small $\varepsilon>0$.

We prove Theorem 3.2 in Appendix A. Note that the bifurcations described in the last sentence of Theorem 3.2 are singular ones. Under these, the local, qualitative behaviour of the solutions of Eq. (8) may change for different small $\varepsilon>0$ as we describe below in an example. On the other hand, periodic orbits arising from either simple zeros or quadratic and $\kappa$-nondegenerate zeros persist for small $\varepsilon>0$. We refer the reader to [62, 63, 64] for detailed analyses of such singular bifurcations.

Figure 4 illustrates the bifurcations described in Theorem 3.2 when the family $\mathcal{P}$ can be locally parametrised with the period, which is also the selected continuation parameter, $\kappa=\tau$. This means sweeping along orbits of the family, indicated with a blue line, and analysing when these orbits give rise to perturbed ones in the frequency response, denoted in red.

Plot (a) shows a saddle-node bifurcation in $\tau$, also known as limit or fold bifurcation. For this type of quadratic zero, the conservative orbit at $\tau_{0}$ and the period $\tau_{0}$ itself are $O(\varepsilon)$-close to a locally unique saddle-node periodic orbit in $\tau$ of the frequency response. This unique orbit originates as two solutions branches of Eq. (8) join together. After this point, conservative orbits of $\mathcal{P}$ do not smoothly persist, at least locally.

The singular case of isola birth, illustrated in Fig. 4(b) and Fig. 4(c), has three possible bifurcation outcomes, depending on the value of the parameters. It is typically observed that either no solution persists from the ones in $\mathcal{P}$ (not show in Fig. 4) or a closed branch of solutions appears, i.e., an isola as shown in Fig. 4(c). Instead, the single solution case of Fig. 4(b) may occur, but it is non-generic.

Similarly, simple bifurcations may manifest themselves in three scenarios. The bottleneck, in Fig. 4(d), and the isola detachment, in (f), are generic, while the node singularity of Fig. 4(e) is an extreme case.

Remark 3.2**.**

The results we have presented in Theorems 3.1-3.2 apply to general, non-autonomous perturbations of conservative systems with a normal family of periodic orbits, not just to mechanical systems. Moreover, the perturbation may be also of type $g(x,\dot{x},t;T,\varepsilon)$.

3.3 The Melnikov function for mechanical systems

The underlying physics of the full system in Eq. (2) implies that any periodic solution with minimal period $lT$ must necessarily experience zero energy balance in one oscillation cycle. Defining the energy function

[TABLE]

we deduce that, along such a periodic orbit, we must have

[TABLE]

given that the work done by the non-conservative forces must vanish over one cycle of that orbit. Equation (13) is commonly called the energy principle in literature [24, 18, 36].

By imposing $\xi=(q_{0}(s;p),\dot{q}_{0}(s;p))+O(\varepsilon)$ and $lT=m\tau+O(\varepsilon)$ in the energy balance equation, one can easily verify that the Melnikov function of Eq. (9) is the leading-order term of the Taylor expansion of Eq. (13), i.e.,

[TABLE]

[TABLE]

where we have used Eq. (4) and the relation

[TABLE]

Before exploring the implications of this peculiar form of the Melnikov function in specific cases, it is useful to recall the definitions of subharmonic and superharmonic resonances [8] in terms of $l$ and $m$. These integers define the relation between the minimal period of the orbit and that of the perturbation. A subharmonic resonance occurs when the forcing frequency is a multiple of the orbit frequency, i.e., $l\neq 1$ and $m=1$. The converse holds for a superharmonic resonance, for which we have $l=1$ and $m\neq 1$. The attribute ultrasubharmonic [39] indicates higher-order resonances, when both $m$ and $l$ are different from 1.

4 Monoharmonic forcing with arbitrary dissipation

Due to their importance in the structural vibrations context, we now consider perturbations $Q$ in Eq. (2) whose leading-order term in $\varepsilon$ is of the form

[TABLE]

where $e\in\mathbb{R}$ is a forcing amplitude parameter, $f_{e}\in\mathbb{R}^{N}$ is a constant vector of unit norm and $C(q,\dot{q})$ is a smooth, dissipative vector field. The actual forcing amplitude and the dissipative vector field are both rescaled by the value of the perturbation parameter $\varepsilon$.

First, we discuss the possible bifurcations that single orbits can experience in such systems when perturbed into forced-damped frequency responses, then we discuss the fate of periodic orbit families. Finally, we also illustrate the implications of the Melnikov method for the phase-lag quadrature criterion used in experimental vibration analysis.

4.1 Bifurcations from single orbits

We consider bifurcations from the conservative orbit $\mathcal{Z}$ with $p\in\mathcal{Z}$ and seek to perform continuation with the parameter $e$. In this case, the Melnikov function takes the form

[TABLE]

where we have introduced the resistance

[TABLE]

measuring the dissipated energy along one period $\tau$ of $\mathcal{Z}$. This function is independent of $s$ since $C$ does not explicitly depend on time and the factor $m$ arises in (17) because (18) is $\tau$-periodic. In contrast,

[TABLE]

is the work done by the force along $m$ periods of the conservative solution.

To simplify Eq. (17) further, we express the conservative periodic solution $\mathcal{Z}$ using the Fourier series

[TABLE]

where $a_{k},\,b_{k}\in\mathbb{R}^{N}$ are the Fourier coefficients of the displacement coordinates. By inserting this expansion in Eq. (17), we obtain for $w^{m:l}(s,e)$ the expression

[TABLE]

where

[TABLE]

We provide the details of these derivations in Appendix B. The quantity $W^{1:l}(e)$ measures the maximum work done by the forcing along one cycle of the conservative periodic orbit. This work depends linearly on the forcing amplitude parameter $e$. Equation (21) implies that superharmonics or ultrasubharmonics cannot occur for the considered perturbation, which is consistent with literature observations. As a consequence, we have the following proposition characterising primary and subharmonic resonances, where the relation between the forcing frequency $\Omega$ and the conservative orbit frequency $\omega$ reads $\Omega=l\omega+O(\varepsilon)$.

Proposition 4.1**.**

The Melnikov function for the perturbation in Eq. (16) takes the specific form

[TABLE]

Assuming $M^{1:l}(s,e_{0})\not\equiv 0$ for some $e_{0}\neq 0$, the following bifurcations of the conservative periodic orbit $\mathcal{Z}$ are possible for small $\varepsilon>0$:

1. (i)

if $|W^{1:l}(e_{0})|<|R|$, the conservative solution $\mathcal{Z}$ does not smoothly persist; 2. (ii)

if $|W^{1:l}(e_{0})|>|R|$, two periodic orbits bifurcate from $\mathcal{Z}$; 3. (iii)

if $|W^{1:l}(e_{0})|=|R|>0$, there exist a forcing amplitude parameter $e_{sn}=e_{0}+O(\varepsilon)$ for which a unique periodic orbit emanates from $\mathcal{Z}$.

Proof.

Since $M^{1:l}(s,e_{0})$ remains bounded away from zero for $|W^{1:l}(e_{0})|<|R|$, statement (i) follows from Theorem 3.1. When $|W^{1:l}(e_{0})|>|R|$, $M^{1:l}(s,e_{0})$ features $2l$ simple zeros for $s\in[0,\tau)$ for which Theorem 3.1 applies again. Considering that the forcing signal passes $l$ times the zero phase in $[0,\tau)$, $\,l$ of these zeros correspond to a single perturbed orbit so that two periodic solutions bifurcate from $\mathcal{Z}$, proving statement (ii). Finally, we will argue that statement (iii) is a direct consequence of Theorem 3.2. First, note that the Melnikov function (23) features $l$ quadratic zeros in $s$ as defined in Eq. (11), corresponding to the $l$ maxima or minima of $\cos\left(l\omega s-\alpha_{l,e}\right)$ for $s\in[0,\tau)$, depending on the signs of $W^{1:l}(e)$ and $R$. Considering a location $s_{qz}$ among these quadratic zeros, we obtain

[TABLE]

by the assumption $|W^{1:l}(e_{0})|=|R|>0$. Thus, these quadratic zeros are nondegenerate in $\kappa$ and, since $l$ of them correspond again to a single orbit, we conclude that a saddle-node bifurcation occurs from Theorem 3.2. More precisely, there exists a value $e_{sn}=e_{0}+O(\varepsilon)$ for which a periodic orbit $O(\varepsilon)$-close to $\mathcal{Z}$ corresponds to a fold point for continuations in $e$.

Due to the specific form of the function in Eq. (23), no further degeneracies in $s$ are possible (e.g. cubic zeros) for $M^{1:l}(s,e)$ so that the cases (i-iii) are the only possible bifurcations. ∎

From Proposition 4.1, we can derive necessary conditions for the persistence of a periodic orbit under forcing and damping. Either for case (ii) and (iii), $W^{1:l}(e)$ must be nonzero, i.e., $e\neq 0$ and $A_{l,e}>0$. The latter quantity is zero if the $l$-th harmonic is not present in Eq. (20) or if $f_{e}$ is orthogonal to both its Fourier vectors. In the non-generic case of $M^{1:l}(s,e_{0})\equiv 0$, the Melnikov function does not give any information on the persistence problem.

4.2 Bifurcations from normal families

We now investigate possible bifurcations that a conservative, 1-normal family $\mathcal{P}$ of periodic orbits may exhibit when perturbed with Eq. (16) into frequency responses at fixed $e$. Specifically, we study phenomena that occur with respect to the forcing frequency $\Omega$ and an amplitude measure $a$ of interest.

We assume that either $\omega$ or $a$ can be locally used as the family parameter $\lambda$ for $\mathcal{P}$ and we denote $\mathcal{B}$ the conservative backbone curve in the plane $(l\omega,a)$. We then introduce the following definition.

Definition 4.1**.**

A ridge $\,\mathcal{R}_{l}$ is the curve in the plane $(e,\lambda)$ identifying the forcing amplitudes and the orbits of $\mathcal{P}$ at which quadratic zeros of $M^{1:l}$ in $s$ occur.

The significance of ridges for frequency responses is clarified by the following proposition.

Proposition 4.2**.**

Assume that $eR(\lambda_{0})>0$ and $A_{l,e}(\lambda_{0})>0$ hold for the periodic orbit $\mathcal{Z}$ identified by $\lambda_{0}$. Then, the explicit local definition of $\,\mathcal{R}_{l}$ becomes $e=\Gamma_{l}(\lambda)$, where

[TABLE]

If $D\Gamma_{l}(\lambda_{0})>0$ (resp. $<0$), then the forced-damped response for $e_{0}=\Gamma_{l}(\lambda_{0})$ shows a maximal (resp. minimal) response with respect to $\lambda$ $O(\varepsilon)$-close to $\mathcal{B}$ at $\mathcal{Z}$. If $D\Gamma_{l}(\lambda_{0})=0$ and $D^{2}\Gamma_{l}(\lambda_{0})>0$ (resp. $<0$), then the forced-damped response for $e_{0}=\Gamma_{l}(\lambda_{0})$ includes an isola birth (resp. simple bifurcation) in $\lambda$ which is $O(\varepsilon)$-close to $\mathcal{B}$ at $\mathcal{Z}$.

Proof.

We rewrite the Melnikov function as

[TABLE]

which features $l$ quadratic zeros in $s$ for $e=\Gamma_{l}(\lambda)$. When $D\Gamma_{l}(\lambda_{0})\neq 0$, Theorem 3.2 identifies a saddle-node bifurcation because

[TABLE]

at any of the $l$ locations $s_{qz}(\lambda_{0})$ of quadratic zeros of $M^{1:l}$ in $s$. As already discussed in Proposition 4.1, there exists a unique periodic orbit $O(\varepsilon)$-close to $\mathcal{Z}$, corresponding to a fold for continuations in $\lambda$. If $D\Gamma_{l}(\lambda_{0})>0$, we can choose a small positive $\epsilon$ defining a $\lambda_{1}=\lambda_{0}-\epsilon$ for which

[TABLE]

so that, according to Proposition 4.1, two periodic orbits bifurcate at $e_{0}=\Gamma_{l}(\lambda_{0})$ from the orbit of $\mathcal{P}$ described by $\lambda_{1}$. For $\lambda_{2}=\lambda_{0}+\epsilon$, we can similarly conclude that no orbit persists smoothly. Thus, the periodic orbit at the fold in $\lambda$ represents a maximal response. An analogous reasoning holds for the minimal response case arising for $D\Gamma_{l}(\lambda_{0})<0$.

The last statement of Proposition 4.2 holds again, based on Theorem 3.2, since we have that D_{\lambda}M^{1:l}\big{(}s_{qz}(\lambda_{0}),\Gamma_{l}(\lambda_{0}),\lambda_{0}\big{)}=0 from Eq. (27) and

[TABLE]

∎

Ridges, as introduced in Definition 4.1, are effective tools for the analysis of forced-damped responses in the vicinity of backbone curves as they can track fold bifurcations and generations of isolated responses. These phenomena are the most generic bifurcations for the perturbation type in Eq. (16). Ridge points may be used to detect further singular bifurcation behaviours under additional degeneracy conditions on $\lambda$ [63].

4.3 The phase-lag quadrature criterion

We now discuss the relevance of the phase of the Melnikov function and the next proposition illustrates an important result in this regard.

Proposition 4.3**.**

Consider a perturbed periodic orbit $q_{qz}(t;\xi,T,\varepsilon)$ corresponding to a quadratic zero of the Melnikov function (23) related to the conservative limit $\mathcal{Z}$. Then, the $l$-th harmonic of the function $\langle q_{qz}(t;\xi,T,\varepsilon),f_{e}\rangle$ has a phase lag (resp. lead) of $\,90^{\circ}+O(\varepsilon)$ with respect to the forcing signal if $eR>0$ (resp. $eR<0$).

Proof.

To determine the phase lag, we consider, without loss of generality, the phase condition for $\mathcal{Z}$

[TABLE]

under which the $l$-th term in the Fourier series of the function $\langle q_{0}(t;p),f_{e}\rangle$ is equal to $\langle a_{l},f_{e}\rangle\cos(l\omega t)$, having the same phase of the forcing. In that case, the Melnikov function becomes

[TABLE]

Eq. (31) shows that the $l$ quadratic zeros of the Melnikov function occur for $|W^{1:l}(e)|=|R|$ and $l\omega s_{qz}=-\mathrm{sign}(eR)\pi/2+2k\pi$ with $k=0,1,...\,l-1$. Thus, we obtain

[TABLE]

whose $l$-th harmonic is equal to $\langle a_{l},f_{e}\rangle\cos(l\omega t-\mathrm{sign}(eR)\pi/2)+O(\varepsilon)$, independent of $k$. ∎

In numerical or experimental continuation, one can track the relation between the forcing amplitude parameter $e$ and either the amplitude $a$ or the forcing frequency $\Omega$ under the phase criterion of Proposition 4.3. The resulting curve of points is an $O(\varepsilon)$-approximation of the ridge curve $\,\mathcal{R}_{l}$ whose interpretation is available in Proposition 4.2.

Proposition 4.3 relaxes some restrictions of the phase-lag quadrature criterion derived in [34]. Indeed, Eq. (20) allows for arbitrary periodic motion, not just synchronous ones along which all displacement coordinates reach their maxima at the same time. Moreover, our criterion is not limited to velocity-dependent, odd damping, but it admits arbitrary, smooth dissipations. For this general case, we proved that the phase-lag must be measured in co-location, i.e., when the output (displacement response) is observed at the same location where the input (force) excites the system.

5 Examples

In this section, we study a conservative multi-degree of freedom system subject to non-conservative perturbations in the form of Eq. (16). First, we consider frequency responses with monoharmonic forcing and linear damping. Then, we introduce nonlinear damping to investigate the presence of isolas. In both cases, we show how the Melnikov analysis can predict forced-damped response bifurcations under a $1:1$ resonance between the forcing and periodic orbits of the conservative limit.

We analyse a system composed of six masses $m_{i}$ with $i=1,\,2,\,...\,6$ that are connected by seven nonlinear massless elements, as shown in Fig. 5. All masses are assumed unitary and the external forcing acts on the first degree of freedom only. The seven nonlinear elements exert a force depending on the elongation $\Delta l$ and its speed $\dot{\Delta l}$, modelled as

[TABLE]

for $i=1,\,2,\,...\,7$. The coefficients $k_{i,1}$, $k_{i,3}$ and $k_{i,5}$ are reported in the table in Fig. 5, while the values of $\alpha$, $\beta$, $\gamma$ and $\varepsilon$ will vary from case by case below. The equations of motion read

[TABLE]

To compute conservative periodic orbits and frequency responses for system (34), we use the MATLAB-based numerical continuation package coco [33]. We specifically exploit its periodic orbit toolbox that solves the continuation problem via collocation. In this method, solutions to the governing ordinary differential equations are approximated by piecewise polynomial functions and continuation is performed using a refined pseudo-arclength algorithm.

First, we focus on the study of the conservative limit ($\varepsilon=0$), in which the nonlinear elements are springs with convex potentials and the origin is an equilibrium whose eigenfrequencies $\omega_{i}$ are reported in the table of Fig. 5. Since no resonance arises among these frequencies, the system features six families of periodic orbits emanating from the origin by the Lyapunov subcenter manifold theorem [29]. Using numerical continuation starting from small-amplitude linearised periodic motions, we compute the conservative backbone curve for each mode, shown in Figure 6(a). We plot these curves using the normalised frequency $\bar{\omega}=\omega/\omega_{1}$ and the $L^{2}$ norm $||x||_{L^{2},[0,T]}$ of the conservative periodic orbits. We consider the latter norm as the amplitude measure $a$. With the exception of the first periodic orbit family, the monodromy matrix of the periodic orbits has two Floquet multipliers equal to $+1$, whose geometric multiplicity is 1 in the selected frequency-amplitude range. Therefore, these five orbit families are $1$-normal, precisely belonging to case (a) of Definition 3.1, and showing a hardening trend ($Da$, $D\omega>0$). The first family also shows normality with hardening behaviour up to the magenta point, where branching phenomena takes place and a further family originates from the continuation of the first linearised mode. As $1$-normality does not hold in the vicinity of the branch point, depicted in magenta in Fig. 6(a), we restrict our analysis of the first family to amplitudes below the branch point amplitude.

5.1 Resonant external forcing with linear damping

In this first example, we take the damping to be linear with $\alpha=0.04$ and $\beta=\gamma=0$ in Eq. (33). We focus on the orbits surviving from primary resonance conditions, where $m=l=1$, and we perform the analysis of $M^{1:1}$ for each mode of the system as described in section 4.

Figure 6(b) shows the work done by non-conservative contributions along the first modal family of conservative orbits parametrised with the non dimensional frequency $\bar{\omega}$. The black solid line is the resistance $R(\bar{\omega})$, while coloured lines represent $W^{1:1}(e,\bar{\omega})$ for three forcing amplitudes. According to Proposition 4.1, we find that two orbits bifurcate from the conservative one when the lines illustrating $W^{1:1}(e,\bar{\omega})$ lay in the grey zone of this plot, i.e., when $W^{1:1}(e,\bar{\omega})>R(\bar{\omega})$. No orbit bifurcates in the white area and unique solutions appear at intersection points between coloured lines and the black one. We also note that $A_{1,e}$ is never zero, except when $\bar{\omega}=1$. Similar trends and considerations hold for the other modes, except for the last one. For that mode, the active work contribution of the forcing is very small compared to the dissipative terms: the forcing is nearly orthogonal to the mode shape. Thus, no orbits arise from the conservative limit for the forcing amplitude ranges under investigation.

Figure 7(c) shows the curves $\Gamma_{1}(\bar{\omega})$ for the first five modes of the system using different colours; all of them show a strictly increasing monotonic trend. Thus, according to Proposition 4.2, ridge orbits are $O(\varepsilon)$ approximations for maximal responses in $\omega$ and $a$, since all the conservative backbone curves can be parametrised with both quantities. By selecting a forcing value in Fig. 7(c), we can predict the frequencies and the amplitudes of the maximal responses in the forced-damped setting. Moreover, since the damping is linear and proportional, ridges are defined as

[TABLE]

where we denoted with $K$ the stiffness matrix of system (34) and expressed initial conditions $p$, periods $\tau$ and coefficients $A_{1,e}$ as functions of the non-dimensional frequency $\bar{\omega}$. From Eq. (35), we obtain that the location of maximal frequency responses close to backbone curves is determined by the ratio between the forcing amplitude parameter $e$ and the damping term $\alpha$, with $O(\varepsilon)$ accuracy.

These theoretical findings are confirmed by the direct numerical computation of frequency responses presented in Fig. 7. To obtain them, we continue in frequency an initial guess acquired through numerical integration for a forcing frequency away from resonance with any of the linearised natural frequencies. The existence of this orbit is guaranteed by the asymptotic stability of the origin when $\varepsilon>0$ and $e=0$. Plot (a) in Fig. 7 shows two frequency sweeps for $e=1$ and for $\varepsilon=0.05$ (grey line), $0.1$ (black line), while this plot is zoomed in (b) around the first and fifth peaks. The sixth mode shows some tiny responses with the rightmost peaks in these two frequency sweeps, more evident for the case $\varepsilon=0.05$ where the physical damping is lower. Figures 7(a) and 7(b) are completed with our analytic predictions in green for the maxima. By imposing the forcing parameter in the ridges of Fig. 7(c), we obtain the frequencies of each mode around which maximal response occur that are validated when carried over with green dotted lines in Fig. 7(a). Moreover, Fig. 7(d) shows the frequency response surface keeping $e=1$ and varying the damping value222For purposes of better illustration, we decided to sweep with the damping parameter $\alpha$ instead of the forcing amplitude one. for two orders of magnitude. Green curves show analytic predictions for the maxima that closely approximate the peaks of this surface.

5.2 Resonant external forcing with nonlinear damping

We now repeat the analysis of the previous section including also the nonlinear damping characteristic of the connecting elements. In order to break the monotonic trend of the resistance in the linear damping case, cf. Fig. 6(b), we select $\alpha=0.2481$, $\beta=-1.085$ and $\gamma=0.8314$. We also restrict our attention solely to the first mode of the system.

The Melnikov analysis is reported in Figure 8(a), which outlines a behaviour change for increasing forcing amplitudes. Indeed, an isola birth occurs at $e\approx 0.4$ as was also displayed in Fig. 4. The branch persists up to connecting with the main branch for $e\approx 1$ through a simple bifurcation.

These predictions are confirmed by the numerical simulations shown in Figures 8(b) and 8(c). The former illustrates several frequency responses for different physical forcing amplitudes, with $\varepsilon=0.1$. The green line is the ridge $\mathcal{R}_{1}$, also plotted in Fig. 8(c), and the two singular bifurcations show up at its folds, as explained in Proposition 4.2. From a computational perspective, main branches of the frequency response are computed with the same strategy of the previous section. For isolated branches, we obtain initial guesses from a numerical continuation in $e$ of saddle-node periodic orbits333This functionality is directly available in the periodic orbit toolbox of coco [33] through the constructor ode_SN2SN. that started near the maximal response of the frequency sweep at $e=1.3$. We also plot the relation between frequency and forcing amplitude in this latter numerical continuation with the black line of Fig. 8(c). This curve is $O(\varepsilon)$-close to the ridge (in green), which was obtained solely from the knowledge of the conservative limit. We remark that a 18-core workstation with 2.3 GHz processors required 18 minutes and 15 seconds to compute the black curve, while the green curve took 1 minute and 45 seconds to compute.

6 Conclusion

We have developed an analytic criterion that relates conservative backbone curves to forced-damped frequency responses in multi-degree-of-freedom mechanical system with small external forcing and damping. Our procedure uses a perturbation approach starting from the conservative limit to evaluate the persistence or bifurcation of periodic orbits in the forced-damped setting. We have shown that this problem can be reduced to the analysis of the zeros of a Melnikov-type function. In a general setting, we proved that, if a simple zero of the Melnikov function exists, generically two periodic orbits bifurcate the conservative limit. We also characterised quadratic zeros and eventual singular bifurcations that may arise. Our results assume the forcing to be periodic and small, but otherwise allow for arbitrary types of damping and forcing. In addition, our analysis yields analytic criteria for the creation of subharmonics, superharmonics and ultrasubharmonics arising from small forcing and damping.

When applied specifically to mechanical systems, the Melnikov function turns out to be the leading-order term in the equation expressing energy balance over one oscillation period. In this context, we have worked out the Melnikov function in detail for the typical case of purely sinusoidal forcing combined with an arbitrary dissipation. Our method shows that either two, one or no orbits can arise from an orbit of the conservative limit. Moreover, ridge curves allow to identify forcing amplitudes and orbits of conservative backbone curves that are close to bifurcations phenomena of the frequency response. Thus, saddle-node bifurcations of frequency continuations, maximal responses and isolas can be efficiently predicted directly from the analysis of the conservative limit of the system. Our analysis also justifies the phase-lag quadrature criterion of [34] in a general setting, without the assumptions of synchronous motion and linear damping.

We have confirmed these theoretical findings by numerical simulations. Specifically, we have considered a nonlinear mechanical system with six degrees of freedom, and implemented our Melnikov analysis on the six families of periodic orbits emanating from an equilibrium. We have verified our results both for linear and nonlinear damping. In the latter case, we successfully predicted the generation of isolated branches in the frequency response. Our six-degree-of-freedom example illustrates that the analytic tools developed here do not require the conservative limit of the mechanical system to be integrable. Indeed, one can apply the present Melnikov function approach directly to periodic orbit families obtained from numerical continuation in the conservative limit of the system.

Our present analysis is limited to well-defined conservative one-parameter families of periodic orbits subject to small damping and periodic external forcing. Perturbations of degenerate cases or resonance interactions, which can be identified by analysing the monodromy matrix of a conservative orbit, would require the analysis of a more general, multi-dimensional bifurcation function. A rigorous approach for tackling such problems can be found in [57].

Acknowledgements. We are grateful to Harry Dankowicz and Jan Sieber for organizing the Advanced Summer School on Continuation Methods for Nonlinear Problems 2018 with the financial support from the ASME Design Engineering Division and the ASME Technical Committee on Multibody Systems and Nonlinear Dynamics. We are also thankful to Thomas Breunung and Shobhit Jain for their careful proofreading of this paper.

Appendix A Proofs of the theorems in section 3

A.1 Preparatory results

We first need some technical results to set the stage for the proofs of the theorems stated in section 3. For a one-parameter family $\mathcal{P}$ of periodic orbits emanating from a $m$-normal periodic orbit $\mathcal{Z}$, the smooth map $\mathcal{T}:\mathbb{R}\rightarrow\mathbb{R}$ describes the minimal period $\mathcal{T}(\lambda)$ of each orbit. Introducing a Poincaré section $\mathcal{S}$ passing through the point of $z\in\mathcal{P}$, we can find a smooth curve $\vartheta_{\mathcal{S}}:\mathbb{R}\rightarrow V\cup\mathcal{P}\pitchfork\mathcal{S}$ parametrising initial conditions under $\lambda$ where $V\subset\mathbb{R}^{n}$ is a open neighbourhood of $z$. For more detail on these mappings, we refer the reader to [65]. We denote the tangent space of the 2-dimensional manifold $\mathcal{P}$ at the point $z$ by $\mathbb{T}_{z}\mathcal{P}$, to which $f(z)$ belongs due to invariance. We consider vectors as column ones and we use the superscript ∗ to denote transposition. We refer to the column and row spaces of a matrix $A$ with the notations col$(A)$ and row$(A)$, respectively.

Next, we discuss a useful result on the properties of the monodromy matrices for normal periodic orbits. Specifically, we restate Proposition 2.1 of [57] for the setting of $m$-normal period orbits.

Proposition A.1**.**

Consider an $m$-normal periodic orbit $\mathcal{Z}$ of period $m\tau$ in the periodic orbit family $\mathcal{P}$. The smooth invertible matrix families $K,\,R:\mathcal{Z}\rightarrow\mathbb{R}^{n\times n}$, defined as

[TABLE]

satisfy the identity

[TABLE]

where $A_{r}(z)\in\mathbb{R}^{(n-2)\times(n-2)}$ is always invertible and the value $\tau_{v}\in\mathbb{R}$ describes the shear effect within $\mathcal{P}$ being zero if $\mathcal{Z}$ is a normal periodic orbit of case (b) and nonzero for case (a) of Definition 3.1.

Proof.

For proving this factorisation result, we need to characterise kernel and range of the monodromy operator for $\mathcal{Z}$ based at $z$. First, by [61], we have that

[TABLE]

Without loss of generality, we introduce a Poincaré section $\mathcal{S}$ orthogonal to $f(z)$ and the value $\lambda_{z}$ identifies $\mathcal{Z}$ leading to $z=\vartheta_{\mathcal{S}}(\lambda_{z})$, $\tau=\mathcal{T}(\lambda_{z})$. Consider the identity

[TABLE]

whose differentiation in $\lambda$ and evaluation at $\lambda=\lambda_{z}$ yields

[TABLE]

We then have $v(z)=D\vartheta_{\mathcal{S}}(\lambda_{z})/||D\vartheta_{\mathcal{S}}(\lambda_{z})||$ leading to a parametrisation-independent relation and to the definition of $\tau_{v}$ in the form

[TABLE]

If the orbit $\mathcal{Z}$ belongs to the case (a) of Definition 3.1, $\tau_{v}$ cannot be zero, otherwise the kernel of $\Pi^{m}(z)-I$ is two dimensional. Instead, for case (b), $\tau_{v}$ must be zero, otherwise there exists a nonzero vector $v(z)$ whose image is parallel to $f(z)$. In both cases, the column space of $K_{r}(z)$ always lays in the complement of the kernel of $\Pi^{m}(z)-I$ by construction, so it maps through $\Pi^{m}(z)-I$ a $n-2$ dimensional linear subspace $V_{z}$ such that $f(z),\,DH(z)\notin V_{z}$. Since the row space of the matrix $R_{r}(z)$ does not contain the latter vectors, the matrix $A_{r}(z)=R_{r}(z)(\Pi^{m}(z)-I)K_{r}(z)$ is invertible $\forall z\in\mathcal{Z}$. ∎

We can now state and prove the following reduction theorem.

Theorem A.2**.**

Perturbed solutions of Eq. (8) in the form of Eq. (6) are (locally) in one-to-one correspondence with the zeros of the bifurcation function

[TABLE]

where the leading-order term, defined in Eq. (9), is independent from the choice of the mapping $L$ used in the last equation of system (8).

Proof.

With the shorthand notation $z=x_{0}(s;p)$, we consider the following change of coordinates

[TABLE]

where $K(z)$ is the matrix defined in Proposition A.1. By construction, $D\hat{\varphi}(0,\hat{\sigma},s)$ is invertible for any $s\,,\,\hat{\sigma}\in\mathbb{R}$. Then, we rescale $\hat{\delta}=\varepsilon\delta$, $\hat{\sigma}=\varepsilon\sigma$ and, by calling $\eta=(\delta,\sigma,s)$, we denote $\varphi(\eta,\varepsilon)=\hat{\varphi}(\varepsilon\delta,\varepsilon\sigma,s)$. Note that $\varphi(\cdot,\varepsilon)$ is a family of diffeomorphisms for $\varepsilon$ nonzero small enough. Note also that col$(K_{\mathbb{T}\mathcal{Z}}(z))=\mathbb{T}^{\perp}_{z}\mathcal{Z}$.

By imposing this coordinate change and Taylor expanding in $\varepsilon$ Eq. (8), we obtain $\Delta_{l,L}(\varphi(\eta,\varepsilon),\varepsilon)=\varepsilon\Delta_{1}(\eta,\varepsilon)$. The latter mapping is of $C^{r-1}$ class and reads

[TABLE]

in which

[TABLE]

where we have denoted by $x_{\kappa}(m\tau;z,m\tau l,0)$ the solution of the first variational problem in the parameter $\kappa$ at time $m\tau$. The solution of the first variation in the period is zero since the period dependence of the vector field only appears at $O(\varepsilon)$. Exploiting Proposition A.1, we project $\Delta_{1}$ using the invertible matrix $R_{ext}(z)$ defined as

[TABLE]

in order to obtain

[TABLE]

We now aim to show that $A(z)$ is an invertible matrix for any $s$. Due to its block matrix structure, its determinant reads

[TABLE]

where the first factor is nonzero due to Proposition A.1. For the second factor, we use the identity

[TABLE]

whose differentiation in $\lambda$, evaluation $\lambda=\lambda_{z}$ and division by $||D\vartheta_{\mathcal{S}}(\lambda_{z})||$ yields

[TABLE]

proving that $A(z)$ is then invertible. Hence, we can solve for $\delta$ and $\sigma$ in the leading order term of $\Delta_{r}$ for any $s$, so that the implicit function theorem ($\Delta_{1}^{\prime}\in C^{r-1}$ with $r\geq 2$) guarantees that we can locally express $\delta=\delta_{r}(s,\varepsilon)$ and $\sigma=\sigma_{r}(s,\varepsilon)$ such that $\Delta_{r}(\delta_{r}(s,\varepsilon),\sigma_{r}(s,\varepsilon),s,\varepsilon)=0$. Thus, we have shown that the perturbed solutions of $\Delta_{l,L}(\xi,T,\varepsilon)=0$ have a one-to-one correspondence with the zeros of the bifurcation function $B^{m:l}_{L}(s,\varepsilon)$ defined as

[TABLE]

where $M^{m:l}(s)=\langle DH(z),\chi(z)\rangle$. Moreover, this function does not depend on the mapping $L$ used as a constraint in Eq. (8).

We now aim to simplify the Melnikov-type function $M^{m:l}(s)$ to the form in Eq. (9). Denoting with $Y(t;x_{0}(s;p))$ the solution of the first variational problem for the vector field $f(x)$, we write explicitly the solution of the first variational problem in $\varepsilon$ (see [59]) leading to

[TABLE]

and we recall that the dynamics on an energy surface $H(x)=H(p)$ (that acts as a codim. 1 invariant manifold), is characterised by (see Proposition 3.2 in [66] for a proof)

[TABLE]

Equation (53) leads to

[TABLE]

and this function is clearly smooth and $m\tau$-periodic. ∎

A.2 Proof of Theorem 3.1

Proof.

Thanks to Theorem A.2, we are able to reduce the persistence problem of Eq. (8) to the study of $B^{m:l}_{L}(s,\varepsilon)$. The zeros of this function mark the existence of periodic orbits for $\varepsilon$ small enough which smoothly connect to $\mathcal{Z}$ at $\varepsilon=0$. Note that, if the $M^{m:l}(s)\equiv 0$, then no conclusions for persistence can be drawn solely from $M^{m:l}(s)$. Indeed, we need to analyse the $O(\varepsilon)$ term in $B^{m:l}_{L}$.

If the Melnikov function remains bounded away from zero, then we conclude the last statement thanks to the fact that no zeros of the bifurcation function exists for $\varepsilon$ small enough.

We now analyse the case of simple zeros. Assuming that the conditions in Eq. (10) hold for $s_{0}$, the implicit function theorem guarantees that we can express $s=s(\varepsilon)$ from the bifurcation function $B^{m:l}_{L}(s,\varepsilon)$. According to the proof of Theorem A.2, we can define

[TABLE]

such that $\Delta_{l,L}(\varphi(\eta(\varepsilon),\varepsilon),\varepsilon)=0$ for a sufficiently small neighbourhood $C_{0}\subset\mathbb{R}$. Hence, we can express the initial conditions and the periods

[TABLE]

of periodic orbits solving system (3) and satisfying $L(\xi(\varepsilon),lT(\varepsilon))=\lambda$ for small enough $\varepsilon>0$.

Finally, the second statement of Theorem 3.1 is a direct consequence of the intermediate value theorem. Namely, the existence of a simple zero for the Melnikov function implies that $M^{m:l}(s)$ is not constant and there exist points $s_{1}=s_{0}-\epsilon$ and $s_{2}=s_{0}+\epsilon$ such that $M^{m:l}(s_{1})M^{m:l}(s_{2})<0$ for $\epsilon>0$ small enough. Due to periodicity, we also have $M^{m:l}(s_{2})M^{m:l}(s_{1}+m\tau)<0$. Thus, there exists at least another $\hat{s}_{0}\in(s_{2},\,s_{1}+m\tau)$ such that $M^{m:l}(\hat{s}_{0})=0$ due Bolzano’s theorem and it must be a zero at which the function changes sign, i.e., a topologically transverse zero. ∎

Remark A.1**.**

Theorem 3.1 guarantees smooth persistence only. There may be degenerate cases where there exist periodic orbits of system (3) that are still $O(\varepsilon)$-close to $\mathcal{Z}$, but they can only be continuously connected to the latter or not connected at all. The Melnikov function in (9) cannot prove existence of such orbits.

Remark A.2**.**

To analyse the type of zeros of the Melnikov function, it is convenient to have closed formulae for its derivatives. The first derivative can be computed as

[TABLE]

given that

[TABLE]

which is again a smooth periodic function. Thus, a transverse zero $s_{0}$ of $M^{m:l}(s)$ must satisfy:

[TABLE]

Assuming enough smoothness, the second derivative of $M^{m:l}(s)$ is likewise

[TABLE]

A similar formula follows for high-order derivatives.

Remark A.3**.**

Note that if the orbit family can be parametrised with the period, one can directly insert the exact resonance condition into the displacement map. In this case, the method developed in [57] applies in a straightforward way in what the authors call a non-degenerate case. Compared with the discussion in that reference, we simplified the final Melnikov function.

A.3 Proof of Theorem 3.2

Once the reduction to a scalar bifurcation function has been performed as in Theorem A.2, the statements in Theorem 3.2 follow from results of the bifurcation analysis outlined in [63]. Specifically, one can look at Theorem 2.1 and Table 2.3 in Chapter IV to recognise the bifurcation problem. In that reference, the singular bifurcation isola birth is called isola centre.

We further remark that a saddle-node bifurcation persists in the perturbed setting. Indeed, defining

[TABLE]

we find that

[TABLE]

Therefore, the implicit function theorem applies, guaranteeing that a locally unique orbit persists at $s_{sn}=s_{0}+O(\varepsilon)$, $\kappa_{sn}=\kappa_{0}+O(\varepsilon)$.

Appendix B Melnikov function with monoharmonic, space-independent forcing

In this Appendix, we show the derivations that lead to Eq. (21). By substituting the Fourier series of Eq. (20) into (21), we find that

[TABLE]

Expanding using trigonometric addition formulae, we obtain

[TABLE]

We recall the following trigonometric integral identities with $k\neq j$ :

[TABLE]

Thus, the integrals in the first and last summations in Eq. (64) are always zero. We first discuss the case $m\neq 1$. We call $m\tau=\tau_{o}$ so that Eq. (64) becomes

[TABLE]

Therefore, to obtain nonzero integrals in Eq. (66), we need that $km=l$ according to Eq. (65). However, since we choose $l$ and $m$ to be positive integers and relatively prime, that condition will never hold. We then conclude that $w^{m:l}(s,e)\equiv 0$ for $m\neq 1$.

For $m=1$, only the terms for $k=l$ can be nonzero in Eq. (66), resulting in

[TABLE]

Thus, we recover Eq. (21) with the proper definitions of $A_{l,e}$ and $\alpha_{l,e}$.