# Snaking bifurcations in a self-excited oscillator chain with cyclic symmetry

Antonio Papangelo, Aurelien Grolet, Loic Salles, Norbert Hoffmann, Michele Ciavarella

► To cite this version:

Antonio Papangelo, Aurelien Grolet, Loic Salles, Norbert Hoffmann, Michele Ciavarella. Snaking bifurcations in a self-excited oscillator chain with cyclic symmetry. Communications in Nonlinear Science and Numerical Simulation, 2017, 44, pp.108-119. 10.1016/j.cnsns.2016.08.004. hal-02121450

HAL Id: hal-02121450

https://hal.science/hal-02121450

Submitted on 14 Dec 2019

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
# Snaking bifurcations in a self-excited oscillator chain with cyclic symmetry

A. Papangelo*(1,2), A. Grolet(2), L. Salles(2), N. Hoffmann(2,3), M. Ciavarella(1)

(1) Polytechnic of Bari, 70126 Bari, Italy

(2) Imperial College London, Exhibition Road, London SW7 2AZ, UK

(3) Hamburg University of Technology, 21073 Hamburg, Germany

(*) corresponding author. e-mail: antonio.papangelo@poliba.it

## Abstract

Snaking bifurcations in a chain of mechanical oscillators are studied. The individual oscillators are weakly nonlinear and subject to self-excitation and subcritical Hopf-bifurcations with some parameter ranges yielding bistability. When the oscillators are coupled to their neighbours, snaking bifurcations result, corresponding to localised vibration states. The snaking patterns do seem to be more complex than in previously studied continuous systems, comprising a plethora of isolated branches and also a large number of similar but not identical states, originating from the weak coupling of the phases of the individual oscillators.

Key words:

snaking bifurcation, subcritical Hopf bifurcation, localised vibration, self-excitation, bistability, nonlinear dynamics

## 1. Introduction

Spatially localised states of dynamical systems have been studied in a large number of different fields in the sciences and in engineering. While for linear systems Anderson localisation was the key to quite a satisfactory understanding, in nonlinear dynamical systems the quest to understand localisation seems far from settled. For a long time progress seems to have been largely confined to conservative nonlinear systems, where solitons and breathers made their appearance. Only later, dissipative systems have come
into focus, with first work based on tracing solitons into the driven and dissipative regime, introducing dissipative solitons. In parallel to the study of solitary states in conservative and dissipative systems, another breakthrough to the understanding of spatial localisation in dissipative localisation was accomplished in the study of subcritical bifurcations in pattern-forming systems, where the concept of branching has emerged and is a well-established field of study today.

Branching is today well known in a number of disciplines, amongst others in optics [1], granular matter [2], structural mechanics [3][4][5][6][7][8], and mostly in fluid dynamics [1][9][10][11][12][13], and magnetohydrodynamics [14][15]. The first studies into the topic have probably emerged in the field of binary-fluid convection, where spatially localised convection rolls have been observed in water-ethanol mixtures [9] or helium [10]. There localised convection domains of arbitrary length are found to be stable, being surrounded by the conductive state.

In terms of bifurcation diagrams, the localised states have shown to be arranged in a unique and fascinating way, giving birth to what is called a snaking structure [11][12][15][6]. The snaking structure has e.g. been studied a lot in the one dimensional [16] and the two dimensional [17] Swift-Hohenberg equation, which is a convenient and generic model system to study fundamental properties of the arising dynamics.

A typical snaking bifurcation diagram involves two snaking solution branches, intertwined into each other. Figure 1 gives an example from convection [12] to illustrate the phenomenon. For the bifurcation diagram (left panel), the average kinetic energy “E” of the flow is plotted versus the Rayleigh number “Ra” and two intertwined branches appear. In the right panel, nineteen solutions are shown, which correspond to the numbers positioned close to the snaking structure, in which spatially localised convective rolls can be identified. Notice that the higher the energy of the solution, the larger the number of convection cells. Often the two snaking branches are also interconnected through a number of unstable branches, and a ladder like pattern emerges [17].

Although snaking bifurcations are now generally known and studied in many fields of dynamical systems, it seems that there is hardly any study into the phenomenon in the context of structural vibrations in engineering. In many respects this is quite surprising, since non-linear oscillators with subcritical Hopf bifurcations, often coupled to neighbouring oscillators of the same type into chains or arrays, are actually very common models for
a number of systems from engineering vibrations. And also the appearance
of bi- or multi-stability, which is obviously at the core of the phenomenon
[13][16][17], is well established in many of these engineering systems. More-
over, the emergence of spatially localised vibration states in structural dy-
namics is also a well known observational fact: e.g. in turbo-machinery,
there is the so-called effect of 'mis-tuning in rotors'[18][19]. Traditionally,
the origin of the localisation is thought to have its root in slight system in-
homogeneities, leading to linear localisation in the sense of Anderson. From
testing, strong localisation is confirmed, but proper validation of the theory
has up to now not been accomplished in the linear framework. In a sense
it is tempting to hypothesise that one of the key reasons behind might be
the non-linearity involved, which definitely becomes substantial for the large
local vibration amplitudes observed. To the best of our knowledge, in model
systems for turbo-machinery dynamics, snaking behaviour has never been
investigated. Also systems from fluid-structure-interaction, may show weak
non-linearity, Hopf bifurcation, and bi-stability, like models for aerofoil flap
dynamics [20][21][22][23]. Similarly in friction induced vibrations the emer-
gence of snaking could be well expected, with all the necessary ingredients
like flutter instability and bi-stability already known to exist, cf. e.g. [24],
[25].

We will thus consider a model system as simple as we can think of, but
derived from models actually in use in the turbo-machinery community and
the field of fluid-structure interaction and friction-induced vibration. We
choose a chain of (weakly non-linear) oscillators coupled into a linear oscil-
lator chain. For simplicity we close the chain into a cyclically symmetric
ring, which moreover has the advantage of bringing it even closer to models
used widely in turbo-machinery for rotors with a small but finite number of
blades attached. To obtain or model instability and spatially local bifurca-
tion, i.e. the bi-stability of the individual oscillator, we introduce non-linear
damping terms, i.e. non-linear terms depending on velocities. This approach
is heuristic and rather for simplicity at the present stage of understanding,
but can be thought of bringing into our purely structural model the corre-
sponding non-linear forcing and dissipation terms from surrounding flow, or
an involved friction interface.

With a velocity dependent force arising from a fifth order polynomial
representation, our system results in individual uncoupled oscillators show-
ing subcritical Hopf bifurcations and bistability, while the oscillator chain,
i.e. the coupled oscillators, yield snaking bifurcations that we determine
by solving algebraic equations or time-integration. Interestingly, the results
do turn out different to many of the hitherto reported snaking bifurcation
patterns. The snaking observed in our study seems to show more than a
single snake-and-ladders pattern, and much of the pattern as a whole seems
to have disintegrated into isolated branches, now usually called isolas. Re-
viewing where these differences might come from, shows that the key features
where our system is different, can easily be identified. First of all, our system
is discrete. Individual self-excited oscillators are coupled, as e.g. in the work
of Yulin and Champneys in [26], where a one-dimensional periodic array of
optical cavities pumped by coherent light were studied. As in our results,
the effect of discreteness was studied, showing that the pinning region (the
parameter interval where the snaking occurs) gets progressively narrower as
the continuum limit is approached. The second characteristic of our system
that is slightly unusual is its finite size due to the cyclic symmetry, which
does not allow arbitrary wave-numbers or wavelengths to appear, and so
does put a constraint onto the system. Formally similar, Taylor and Dawes
in [27], studied snaking and localised states in spatially discrete problems for
modified periodic or phase shifted boundary conditions. Also in their case,
isola rather than continuous snakes have been observed.
<figure><img src="image_2.png" /><figcaption>Figure 1. An example for snaking in a convection system. Adapted from [12]. Left: typical snaking pattern with two intertwined solution branches in the bifurcation diagram with kinetic energy “E” of the fluid plotted versus the Rayleigh number “Ra”. Right: The stream-functions for solutions as marked in the snaking pattern. The graphs show the spatially localised convection patches.</figcaption></figure>

## 2. The mechanical system

We consider a cyclic system of $N_{dof}$ non-linear oscillators, see Figure 2, which are coupled via a weak linear spring of stiffness $k_Δ$. Each oscillator has mass $m$ and is linked to the ground via a linear spring $k$ and a non-linear damper which introduces a velocity proportional force of the form

$$F_v = -c_1 \dot{x} + c_3 x - c_5 x^3. \quad (1)$$

Here $x$ denotes the displacement of the individual oscillator, $\dot{x}$ the velocity, and we introduced the coefficients $c_1, c_3, c_5$ to parametrise the velocity-dependent force.
<figure><img src="image_3.png" /><figcaption>Figure 2. The model system under study.</figcaption></figure>

The evolution equations for the individual oscillators read

$$m\ddot{x}_n + c_1\dot{x}_n - c_3\dot{x}_n^3 + c_5\dot{x}_n^5 + kx_n - k_\Delta(x_{n+1} + x_{n-1} - 2x_n) = 0, \quad (2)$$

where the stiffness $k_\Delta$ couples the n-th mass with the neighbouring ones. We introduce the quantities $\omega_0 = \sqrt{k/m}$, $\eta_\Delta = k_\Delta/k$, $\xi_i = \frac{c_i}{2\sqrt{km}}$, $\tau = \omega_0 t$, divide the equation (2) using the group $\omega_0^2 x_0$, with $x_0$ a reference displacement, obtaining

$$q_1\ddot{x}_n + q_2\dot{x}_n - q_3\dot{x}_n^3 + q_4\dot{x}_n^5 + q_5\dot{x}_n - q_6(\tilde{x}_{n+1} + \tilde{x}_{n-1} - 2\tilde{x}_n) = 0, \quad (3)$$

where

$$q_1 = 1, \quad q_2 = 2\xi_1, \quad q_3 = 2\xi_3\omega_0^2 x_0^2, \quad q_4 = 2\xi_5\omega_0^4 x_0^4, \quad q_5 = 1, \quad q_6 = \eta_\Delta, \quad (4)$$

and the $\Box$ superposed indicates that the new displacements are dimensionless, $\tilde{x}(\tau) = x(\tau)/x_0$. In (3) we defined a dimensionless time $\tau = \omega_0 t$, which allows to replace $\frac{d}{dt}$ with $\omega_0 \frac{d}{d\tau}$. Notice that we choose $\xi_3, \xi_5 > 0$, thus the third degree term of the velocity-dependent force introduces a destabilizing force into the system, while the fifth degree term tends to stabilize it. Figure 3(a) lists the arbitrary but characteristic parameters that will be used in the next sections. We will choose $\xi_1$, i.e. the linear damping coefficient, as our primary control parameter in a range from -0.4 to +0.6. In this range the velocity dependent force changes its shape in a way such that for low $\xi_1$ values a negative damping is introduced, which is often used in the literature to model self-excited vibrations, such as in fluid- or friction-induced flutter [28] or squeal [29].
### 3. Numerical algorithm

#### 3.1. Harmonic Balance Method (HBM)

The problem to solve is composed of $N_{dof}$ second order differential equations. In this study periodic vibrations will be considered only. We therefore apply the Harmonic Balance Method (HBM) as an efficient numerical technique to obtain an approximation to the steady-state solution of the system. In the following we briefly recall the main steps of the HBM, further details can e.g. be found in [30]. Consider a time dependent signal $x(t)$ and express it in a Fourier series

$$x(t) = \frac{a_0}{2} + \sum_{h=1}^{N_h} (a_h \cos(\omega h t) + b_h \sin(\omega h t)), \quad (5)$$

where $\omega$ is the fundamental frequency, $N_h$ is the number of harmonics considered and

$$a_h = \frac{2}{T} \int_{-T/2}^{T/2} x(t) \cos(\omega h t) dt, \quad (6)$$

$$b_h = \frac{2}{T} \int_{-T/2}^{T/2} x(t) \sin(\omega h t) dt. \quad (7)$$

A set of non-linear algebraic equations is obtained by substituting the dimensionless form of eq. (5) into (3) and projecting the equilibrium equations on the Fourier basis $1, \sin(\rho h \tau), \cos(\rho h \tau)$ for $h = 1, ..., N_h$. The projection gives back a system of algebraic equations in which the unknowns are the Fourier coefficients $a_h, b_h$ of the $N_{dof}$ degrees of freedom (dof) considered. In our particular case the dynamical system is autonomous (there is no forcing term) and therefore the angular frequency of the solution $\omega$ is also treated as an unknown.

#### 3.2. Numerical solution and continuation

The number of algebraic equations to be solved is $(2N_{dof}N_h + 1)$. If only one harmonic is considered ($N_h = 1$) the system of equations for the n-th mass is
$$ \left\{ \begin{array}{l} \frac{5}{8} q_4 \rho^5 (a_n^2 + b_n^2)^2 b_n - \frac{3}{4} q_3 \rho^3 (a_n^2 + b_n^2) b_n + \rho (q_2 b_n - \rho q_1 a_n) + \\ \qquad + a_n (q_5 + 2q_6) - q_6 (a_{n-1} + a_{n+1}) = 0, \\[2ex] \frac{5}{8} q_4 \rho^5 (a_n^2 + b_n^2)^2 a_n - \frac{3}{4} q_3 \rho^3 (a_n^2 + b_n^2) a_n + \rho (q_2 a_n + \rho q_1 b_n) + \\ \qquad - b_n (q_5 + 2q_6) + q_6 (b_{n-1} + b_{n+1}) = 0, \end{array} \right. \qquad (8) $$

where $a_n (b_n)$ is the first harmonic cosine (sine) coefficient. In subsection 4.1 it will be shown that in most parameter ranges even with just one harmonic the approximate solution is very close to the exact one, therefore we will mostly use one harmonic only, which leads to a system of $2N_{dof}$ polynomial equations of the fifth order. The other equation comes from the projection of the equilibrium equation on 1, which would allow for evaluation of the mean displacement over a period. As the non-linearities are odd and there is no constant term in the equations, the mean displacement always vanishes. The system of $2N_{dof}$ polynomial equations has been solved using a Newton-Raphson scheme implemented in the MATLAB® function `fsolve`. For a fixed value of $\xi_1$, a set of initial conditions has been provided using steady state solutions obtained from a time integration algorithm. Localised as well as non localised solutions were used as a starting point for the Newton-Raphson algorithm. To solve the system we set for one mass $b_1 = 0$, which allows us to solve for the fundamental frequency of the solution too. A continuation algorithm, also implemented in MATLAB®, continued the solution using a pseudo arc-length scheme [31] which allowed to follow the solution trajectory even when turning points were encountered [31].

# 4. Simulation Results

## 4.1. Single oscillator dynamics

In this subsection we concentrate our attention on the behaviour of a single oscillator when it is isolated from the rest of the chain. The following parameters will be used:

$$ \omega_0 = 2\pi, \quad x_0 = 1, \quad \xi_3 = 0.3, \quad \xi_5 = 0.1, \quad \eta_{\Delta} = 0. \qquad (9) $$

In Figure 3(b) we plot the maximum potential energy of the mass in dimensionless form: $\tilde{U}_{\max} = \frac{1}{2}\tilde{x}_{\max}^2$, where $\tilde{x}_{\max}$ is the vibration amplitude at the steady-state condition. The results have been obtained using three different approaches: circles represent the result of the Time Integration (TI), the solid
line shows the result obtained using the HBM where $x(t)$ is approximated with only the first harmonic, while the dashed line is the result obtained using the HBM with two harmonics, the first and the third. The choice of the first and third harmonics takes into account the fact that the non-linearities introduced are odd and symmetry breaking bifurcations are not considered in this work. We start focusing on the TI results (circles). When $\xi_1 \gtrsim 0.2$ only one solution exists in which for every initial condition the vibration is damped and vanishes. Decreasing $\xi_1$, in the range $0 < \xi_1 \lesssim 0.2$ another solution appears with $\tilde{U}_{\max} > 0$ which corresponds to a stable limit cycle of finite amplitude. The single oscillator experiences a subcritical Hopf bifurcation, and in the interval $0 < \xi_1 \lesssim 0.2$ two different stable solutions exist. Notice that using TI, only the stable solutions can be found (Figure 3(b) circles), while HBM also allows to obtain the unstable ones and to follow them. If $\xi_1 < 0$ the steady state is not stable and the vibration amplitude grows up to the upper branch with $\tilde{U}_{\max} > 0$. In this range the only stable solution is a limit cycle for which the amplitude grows when $\xi_1$ is further reduced. Figure 3(b) also shows the results obtained using HBM with one harmonic (solid line) and with two harmonics (dashed line). First we note that in the interval $0 < \xi_1 \lesssim 0.2$ another solution appears which represents an unstable limit cycle. Differences in the two curves obtained with HBM appear for the upper branch when $\xi_1$ decreases, as higher vibration amplitude implies higher contribution of the non-linearities. The solution which includes two odd harmonics approximates very well the exact solution obtained with TI scheme. Nevertheless if we just focus on the interval $0 < \xi_1 \lesssim 0.2$ we can see that even the single harmonic approximation is never too far from the exact solution<fnref n="1" />. In the following sections the oscillator chain dynamics is studied in the bistability zone ($0 < \xi_1 \lesssim 0.2$), and all the results presented will be obtained using HBM with one harmonic.

<footnote n="1">¹Also in Figure 8 (left panel) we will plot in the bifurcation diagram TI solutions superposed to the HB solutions to demonstrate that the first order truncation is able to capture the basic features of the exact solution. We are aware that this approximation can smooth out some fine details of the true solution, but the aim here is to focus on the overall result.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3. (a) Dimensionless velocity dependent force plotted versus the velocity for $\xi_1 = [-0.4, -0.25, -0.1, ..., 0.6]$, $\xi_3 = 0.3$, $\xi_5 = 0.1$, $\omega_0 = 2\pi$, $x_0 = 1$. (b) Maximum potential energy of one single oscillator plotted against the bifurcation parameter $\xi_1$. Blue dots: results of the time integration at the steady state. Solid line: HB solutions with a single harmonic. Dashed line: HB solutions with two odd harmonics.</figcaption></figure>

## 4.2. Oscillator chain dynamics: linear system

Before studying the dynamical behaviour of a cyclically symmetric chain of $N_{dof} = 12$ non-linear oscillators we analyse the response of the underlying undamped linear system with a dimensionless coupling stiffness $\eta_\Delta = 0.01$. Due to the cyclic symmetry this system presents mostly pairwise degenerate eigenfrequencies that can be computed as [32]

$$ \omega_p = 1 + 2\eta_\Delta \left(1 - \cos\left(\frac{2p\pi}{N_{dof}}\right)\right), \quad (10) $$

where $p \in \mathbb{N}$ and $0 \le p \le N_{dof}/2$ for even $N_{dof}$ or $0 \le p \le (N_{dof} - 1)/2$ for odd $N_{dof}$. Due to the weak coupling among the oscillators the natural frequencies will lie on a narrow band. A possible set of orthogonal normal
modes is [32]

$$ \phi_0 = [1, 1, \dots, 1]^T, \qquad (11) $$

$$ \phi_p^c = [\cos(\theta_p), \cos(2\theta_p), \dots, \cos(N_{dof}\theta_p)]^T, \qquad (12) $$

$$ \phi_p^s = [\sin(\theta_p), \sin(2\theta_p), \dots, \sin(N_{dof}\theta_p)]^T. \qquad (13) $$

The following table lists the natural frequencies [rad/s] and Figure 4 shows the normal mode shapes of the undamped linear system.

$$ \begin{aligned} \omega_0 &= 1, & \omega_1 &= 1.001, & \omega_2 &= 1.005, & \omega_3 &= 1.010, \\ \omega_4 &= 1.015, & \omega_5 &= 1.019, & \omega_6 &= 1.020 \end{aligned} $$

<figure><img src="image_5.png" /><figcaption>Figure 4. Mode shapes of the underlying linear system.</figcaption></figure>

### 4.3. Oscillator chain dynamics: non-linear system

We now study a cyclically symmetric chain of $N_{dof} = 12$ non-linear oscillators in the bistability zone. The parameters used are the same we used for the single oscillators above, except for $\eta_{\Delta} = 0.01$, which introduces a small coupling between the oscillators. The results obtained from time integration were used to derive initial conditions for the continuation algorithm. Figure 5 shows, in the left panel, the sum of the maximum potential energy of each
mass $\tilde{U}_{\max} = \frac{1}{2} \sum_{i=1}^{N} \tilde{x}_{\max,i}^2 = \frac{1}{2} \sum_{i=1}^{N} (\tilde{a}_i^2 + \tilde{b}_i^2)$ plotted against the linear damping coefficient $\xi_1$. Many solutions appear to be entangled, making it almost impossible to distinguish one from the other. Looking more closely at the overall structure created by the superposed solutions we can observe trajectories similar to snake and ladder branches [11][12][6][15] and twelve 'steps' (corresponding to the number of oscillators) can be identified. Each step is labelled with a red circle and the corresponding energy distribution is plotted in the twelve bar plots on the right-hand side of Figure 5. In each bar plot the mass number has been reported on the x-axis, while the bar height is computed as

$$ \tilde{u}_i = \frac{\tilde{x}_{\max,i}^2}{\max(\tilde{x}_{\max,i}^2)_{i=1,\dots,12}}, \quad (14) $$

such that the energy of each oscillator is normalized with respect to the one which is vibrating with the largest amplitude.

<figure><img src="image_6.png" /><figcaption>Figure 5. Left: bifurcation diagram for the non-linear oscillator chain. The subcritical Hopf bifurcation of the individual oscillator is indicated, and the complex snaking pattern linking the spatially homogeneous stationary static state with the state where all oscillators are vibrating fills the zone of bistability. Middle and Right: the average dimensionless energy of each mass for the 12 equilibrium solutions which are marked with a red circle in the snaking pattern.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. In each subplot from 1 to 12 the vibration shape, i.e. with the corresponding phase information, of the solutions in Figure 5 (left panel, red circles) is shown.</figcaption></figure>

The bar plots 1–12 in Figure 5 shed some light on the system dynamics: each step can be easily related to the number of masses on which the vibration is localised. From the bottom to the top, at the first step one mass is moving while the others are more or less motionless, at the second step two masses are moving and so on up to the twelfth. This resembles the usual snaking behaviour, where, for example in fluid dynamics, the steps can be related to the increasing number of convection rolls. Especially in the bottom part of our snaking (Figure 5, left panel) it is possible to see very close similarities with the classical snaking picture. On the other hand, when increasing the energy, and thus the influence of the non-linearities, the picture gets more and more distorted and only seems to be bounded by the very last branch in which all the masses are moving.

Consequently one might conjecture that the relative phases between the individual oscillators play a special role in our system. Figure 6 thus shows the vibration mode shape, indicating the relative phasing between the oscillators, in 12 bar plots numbered from 1 to 12, each-one related to one of the red circles in the bifurcation diagram (Figure 5, left panel). It seems that there is an additional sub-structure due to relative in- or out-of-phase behaviour within our localised vibration zones. Sometimes an oscillator is out of phase
with both of its neighbours, sometimes two adjacent oscillators go in phase, but out of phase to their respective other neighbours. At the present state of our study, we have not yet succeeded to obtain a deeper understanding of this result, but think it is at the very heart of the multiplicity of solutions observed. Further studies will need to follow.

At this point it also seems useful to have another look at the numerical accuracy of our computational approach. Figure 7 shows in each column results obtained from TI for four solutions using $\xi_1 = 0.1$ and different initial conditions. In particular from the left to the right the vibration is localised on one mass, five masses, eight masses and on all the twelve masses. In the first row the displacement time history is shown at the steady-state, in the second row the solution is shown in the phase space and in the third row the dimensionless maximum potential energy of each mass is given (TI results, dark blue). To assess the effect of the first harmonic truncation on the energy distribution among the different oscillators, HBM results are computed for the same parameter sets and using the TI solutions as starting point. The maximum potential energy, shown in Figure 7 (HBM results, pale yellow bars), proves that the energy is distributed exactly in the same manner and that the low-order approximation causes only a very small reduction in the potential energy of each mass (due to the loss of solution details) which seems acceptable for our purposes. Notice also that the shape of the solution with all the masses vibrating corresponds very well to the last mode of the linearised system, in which the motion is out-of-phase.
<figure><img src="image_8.png" /><figcaption>Figure 7. Each one of the four columns corresponds to one solution obtained with TI. In the first row the displacement time history, in the second row the space state trajectories ($\tilde{x}$ vs $\tilde{x}$), in the third row the dimensionless maximum potential energy of each mass in steady state condition is shown (dark blue). In the third row, the energy distribution obtained with single harmonic HBM is plotted (pale yellow).</figcaption></figure>

### 4.4. Isolas: closed solution branches

Figure 5 has suggested that our snaking picture in the bistability zone is more complex than the usual one. In Figure 8 (left panel) the snaking picture observed before is thus plotted again, but this time the ensemble of all the solutions is drawn in the background in pale gray, while four solutions are highlighted using thick lines. It turns out that they have a well defined shape resembling a figure eight. Hereafter we will refer to these branches as isolas, since they are isolated. Despite the whole picture seems really intricate, it is made up of a number of such isolas. Closer study shows that our isolas are actually not strictly speaking isolated in the sense that there are no other solutions branching off, but that there are further bifurcations and interlinking solutions and connecting branches involved. We will come back to this point below.
In obtaining the snaking picture we found and plotted 53 isolas and connecting branches. Following the branches of the isolas, localised solutions with a different number of vibrating masses are linked together, which is conceptually linked to the snaking phenomenon. In Figure 8 four isolas are selected with the peculiarity to have all the masses moving in out-of-phase manner. On the left panel eight points, lying on the straight line $\xi_1 = 0.15$, are marked with a red circle. The dimensionless vibration shape of each point is graphed in the correspondent bar plot on the right side. Focusing on the pair of points that belongs to the same isola, we can see that increasing the energy of the vibration leads to a higher number of masses involved: for example the first isola links solutions localised on one and three masses vibrating, the second isola links solutions localised on three and five masses vibrating, and so on. Compare the points 2 – 3 or 6 – 7: they belong to different isolas, nevertheless they are close in the bifurcation diagram. It can be seen that in fact those points have a slightly different vibration shape. Notice that the points 4 and 5 deserve some further attention: both of them involve five vibrating masses, but the solution has a different vibration shape (compare the bar plots of Figure 8). These considerations indicate how complicated the snaking phenomenon for our vibration system is, as not only the number of vibrating masses matters, but also the relative phases, i.e. the vibration shape plays a fundamental role.

Again, and just to double-check, in the bifurcation diagram (Figure 8, left panel) points obtained via time integration are denoted as red stars to show which branch of the isola is stable. Discrepancies between TI and HBM results are apparent and due to the 1 harmonic approximation, but small.
<figure><img src="image_9.png" /><figcaption>Figure 8. Left: The overall snaking pattern is kept in the background, while four isolas are highlighted. On each isola two points are marked with a red circle and a number. Stars indicate time integration results. Right: eight bar plots visualise the shape of the vibration for the corresponding solution marked.</figcaption></figure>

4.5. *Isolas and connecting branches superposition*

In this subsection we study in some more detail the bifurcation structures
of the isolas found. In Figure 9 the left panel shows on the background the
snaking picture (pale gray), while six solutions are emphasised to show on
which part of the pattern they lie. We chose isolas in the bottom part of
the pattern because they are less distorted than those in the upper part, and
the picture is thus easier to explain. Nevertheless any other choice of the
isolas would show the same main features. The six solutions are drawn in
pairs, using a solid and a dashed line, into the subplots (a), (b), (c) to show
to the reader how they are linked together. Looking at the subplot (a) one
can see that the left hairpin bends of the isolas are connected through an
independent branch which bifurcates from the isola. The same behaviour
can be observed in Figure 9(b) where the branch drawn with the dashed line
bifurcates from the isola drawn with the solid line. Figure 9(c) shows two
solutions which are very close to each other along one branch. All these six
solutions are then superposed to each other in Figure 9(d) to show how the
complex overall snaking pattern emerges. Looking at the six solutions as a
whole, it can be clearly seen that the snaking branches appear and even the
peculiar structure of the 'ladder' (in the sense that they connect two different
points of the same isola) is there, indicated by arrows in Figure 9(d).

In sum it almost seems that the usual snaking picture to be found in more
strongly non-linear, and perhaps more strongly dissipative systems, is kind
of broken into smaller elements, i.e. isolas and connecting branches. Further
work on clarifying these aspects is necessary.

<figure><img src="image_10.png" /><figcaption>Figure 9. Left: Bifurcation diagram in the plane $(\tilde{U}_{\max}, \xi_1)$. The snaking structure is left in the background, while six families of solutions are highlighted (thick line). On the right-hand side the six solutions are plotted in pairs (using one dashed and one solid line) into the subplots (a), (b), (c). In the subplot (d) all the six solutions are plotted together to show how they arrange in the overall snaking structure.</figcaption></figure>

4.6. *Vibration shapes for different solution branches*

Here an analysis of the vibration shapes is carried out to show how the
solution branches shown in Figure 9 relate to different vibration shapes. In
Figure 10 we plot in the left panel the solution branches again (from Figure
9(a)) in the plane $(\tilde{U}_{\max}, \xi_1)$. On the isola of solutions (Figure 10) eight
points labeled from 1 to 8 are marked with a red circle. We note that the
points 8-1-2-3 belong to a part of the isola on which the vibration remains
mostly localised on one single mass. After the solution number 3 a hairpin
curve brings us to another part of the isola on which the solutions 4-5 show a
vibration localised on three masses which vibrate in an out-of-phase manner.
Then, for the solutions 6 and 7 the vibration of two of the masses reduces
again and the vibration comes back to be localised one on just one mass.
We look now at the thick branch of solutions (Figure 10 left panel) which
connects two points of the same isola (thin solid line). The solution 11 shows
a vibration localised on one mass in a similar manner as for the solutions 7
and 8 on the previous isola. Moving towards the solutions 12-13-...-18 the
neighbouring mass starts vibrating in phase with the previous one and the
vibration remains localised on two masses (Figure 10) up to the intersection
between the connecting branch and the isola. To sum up: while the isola
comprises solutions localised on one or three masses, the connecting branch
corresponds to solutions localised on one or two masses.

<figure><img src="image_11.png" /><figcaption>Figure 10. Left: two solution branches, or isolas, are drawn (Figure 9(a)) in the plane $(\tilde{U}_{\max}, \xi_1)$. On each branch eight points are marked and labelled. For each point the corresponding bar-plot shows the shape of the vibration.</figcaption></figure>

Figure 11 shows in the left panel the branches presented in Figure 9(c). In
this case the two patterns have both a figure eight shape and they approach
each other in the middle zone. The subplots on the right-hand side of Figure
11 give the shape of the vibration, while in Figure 12 the dimensionless energy
of each mass $\tilde{u}_i$ is plotted for the points marked with squares/circles in Figure 11 (left panel). The isola of solutions drawn with a thin solid line connects solutions localised on two masses with others localised on three masses (Fig. 11 - Fig. 12). The thick solid line isola, instead, links solutions with three vibrating masses with others with four vibrating masses. Moving to the top of the isola, a smooth transition happens from one kind of localisation to the other. It is interesting to compare the solutions 4-5-6-7 with the solutions 18-11-12-13. The two branches are very close together but are not coincident. The solutions (18-4) (11-5) (12-6) (13-7) have the same number of masses involved in the vibration, but with a different shape. Figure 12 shows that even if the overall energy of the solutions (in pairs) is broadly the same, it is just distributed differently among the masses due to the different shape. For example if in the subplot 18 of Figure 12 the mass 7 is exchanged with the mass 8 the same distribution that appears in the subplot 4 is obtained, and the same reasoning can be done for the other three pairs, and so on for the whole two branches of the two isolas. Note that this is a feature that enriches the general picture of the snaking phenomena in vibrating systems as the appearance of the many different branches could be thought of to arise from some symmetries in the system. This could explain even why in the upper part of our snaking structure (Figure 5 left panel) the solutions appear more entangled: the larger the number of the masses involved in the vibration, the larger the possibilities to arrange them in different ways.

<figure><img src="image_12.png" /><figcaption>Figure 11. Left: two solution branches (Figure 9(c)) in the plane</figcaption></figure>
($\tilde{U}_{\text{max}}, \xi_1$). On each trajectory eight points are marked and labelled with a number. For each number the corresponding bar plot shows the shape of the vibration.

<figure><img src="image_13.png" /><figcaption>Figure 12. Dimensionless energy plotted for each solution marked with the corresponding number from Figure 11.</figcaption></figure>

**5. Conclusions**

In this work we have studied snaking bifurcations of a non-linear cyclically symmetric oscillator chain. Bistability has been introduced in a heuristic manner through non-linear velocity dependent forces. Solutions have been obtained by time integration and harmonic balance techniques. The bifurcation diagrams resulting in the bistability zone resemble typical snaking patterns, but also show marked differences. The solution branches are composed of isolas, which have a figure eight shape in the bifurcation diagram. When the isolas are put together, they picture a typical snaking pattern, which has been observed in many other fields of physics. Still, our findings suggest that the snaking behaviour in structural dynamics could be more complicated due to the superposition of different non-linear mode shapes:
solutions which have very different shapes present almost the same energy content and thus the corresponding solution branches overlap or are very close to each other in the bifurcation pattern.

Physically the snaking phenomenon is due to the weak nearest neighbor coupling among the oscillators and to the positive linear damping coefficient ($\xi_1 > 0$) which damps the small oscillations. In fact when one or few masses are in the high amplitude limit cycle the neighbors can be still stable under small oscillations and don't necessary jump on the high energy solution (limit cycle). We obtain therefore a multitude of solutions (from which, the snaking behavior) because of this bistability between low and high amplitude limit cycles.

From the present results it has become clear that more work on snaking in engineering structures is necessary. Future work will need to focus both on conceptual aspects, as well as on evaluating engineering relevance and impact. As for conceptual understanding, more work is necessary on understanding analogies and differences between the systems under study here and elsewhere. The larger number of states found in the present work needs to be clarified, and the relevance of symmetries, symmetry breaking and imperfect bifurcations, as observed here, needs to be understood better. As for engineering relevance, probably first another set of simplified and idealised model systems, ideally extracted from fluid-structure-interaction, friction-induced vibration, or similar fields, will need to be studied. Thereafter more realistic models, e.g. derived from larger scale computer assisted modelling systems (like finite element analysis or computational fluid mechanics) could be studied.

As for our personal impression about future impact, a deeper understanding of snaking phenomena in engineering systems could make engineers able to better and quantitatively more accurately predict localised non-linear vibration states, which are known to be a source of numerous issues in engineering and technology, like fatigue, strength, or noise. For these reasons we hope that this work might serve as a starting point to conduct further studies on localisation phenomena and snaking in vibration engineering.

## 6. Acknowledgements

A.P. is grateful to the Italian Ministry of Education, Universities and Research, which funded his PhD research. Part of the work has been supported by Deutsche Forschungsgemeinschaft in project HO 3852/11-1.
References

[1] Champneys AR. Homoclinic orbits in reversible systems and their appli-
cations in mechanics, fluids and optics. Physica D: Nonlinear Phenom-
ena, 112.1 (1998):158-186.

[2] Umbanhowar PB, Melo F and Swinney, HL. Localized excitations in a vertically vibrated granular layer. Nature, (1996) 382(6594), 793-796.

[3] Champneys AR, and Thompson JMT. A multiplicity of localized buckling modes for twisted rod equations. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, (1996.) 452. No. 1954.

[4] Thompson, JMT and Champneys AR. From helix to localized writhing in the torsional post-buckling of elastic rods. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. Vol. 452. No. 1944. The Royal Society, 1996.

[5] Lord, GJ, Champneys AR and Hunt GW. Computation of localized post buckling in long axially compressed cylindrical shells. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 355.1732 (1997): 2137-2150.

[6] Hunt GW, Peletier MA, Champneys AR, Woods PD, Wadee MA, Budd CJ, and Lord GJ. Cellular buckling in long structures. Nonlinear Dynamics (2000) 21.1 : 3-29.

[7] Grolet A and Thouverez F. Vibration analysis of a nonlinear system with cyclic symmetry. Journal of Engineering for Gas Turbines and Power, (2011), 133(2), 022502.

[8] Grolet A and Thouverez F. Free and forced vibration analysis of a non-linear system with cyclic symmetry: Application to a simplified model. Journal of sound and vibration, (2012), 331(12), 2911-2928.

[9] Niemela, JJ, Ahlers G and Cannell DS. Localized traveling-wave states in binary-fluid convection. Physical review letters, 64.12 (1990): 1365.

[10] Woodcraft AL, Lucas PGJ, Matley RG and Wong WYT. Visualisation of convective flow patterns in liquid helium. Journal of low temperature physics, (1999) 114(1-2), 109-134.
[11] Batiste O, Knobloch E, Alonso A and Mercader I. Spatially localized binary-fluid convection. Journal of Fluid Mechanics, (2006) 560, 149-158.

[12] Beaume C, Bergeon A and Knobloch E. Homoclinic snaking of localized states in doubly diffusive convection. Physics of Fluids, (2011), 23(9), 094102.

[13] Thual O, and Fauve S. Localized structures generated by subcritical instabilities. Journal de Physique 49.11 (1988): 1829-1833.

[14] Blanchflower S. Magnetohydrodynamic convectons. Physics Letters A 261.1 (1999): 74-81.

[15] Dawes JHP. Localized convection cells in the presence of a vertical magnetic field. Journal of Fluid Mechanics 570 (2007): 385-406.

[16] Burke J and Knobloch E. Homoclinic snaking: structure and stability. Chaos: An Interdisciplinary Journal of Nonlinear Science 17.3 (2007): 037102.

[17] Avitabile D, Lloyd DJ, Burke J, Knobloch E and Sandstede B. To snake or not to snake in the planar Swift-Hohenberg equation. SIAM Journal on Applied Dynamical Systems, (2010) 9(3), 704-733.

[18] Whitehead, DS. Effect of mistuning on the vibration of turbo-machine blades induced by wakes. Journal of Mechanical Engineering Science, 8.1 (1966): 15-21.

[19] Ewins, DJ. The effects of detuning upon the forced vibrations of bladed disks. Journal of Sound and Vibration, 9.1 (1969): 65-79.

[20] Liu JK and Zhao LC. Bifurcation analysis of airfoils in incompressible flow. Journal of Sound and Vibration, (1992) 154(1), 117-124.

[21] Lee BHK, Price SJ and Wong YS. Nonlinear aeroelastic analysis of airfoils: bifurcation and chaos. Progress in aerospace sciences 35.3 (1999): 205-334.

[22] Chen YM and Liu JK. Supercritical as well as subcritical Hopf bifurcation in nonlinear flutter systems. Applied Mathematics and Mechanics, (2008) 29, 199-206.
[23] Pereira DA, Vasconcellos RM, Hajj MR and Marques FD. Insights on aeroelastic bifurcation phenomena in airfoils with structural nonlinearities. Mathematics in Engineering, Science & Aerospace (MESA), (2015) 6(3).

[24] Weiss C, Morlock MM and Hoffmann N. Friction induced dynamics of ball joints: Instability and post bifurcation behavior. European Journal of Mechanics-A/Solids, (2014) 45, 161-173.

[25] Gräbner N, Tiedemann M, Von Wagner U and Hoffmann N. Nonlinearities in Friction Brake NVH-Experimental and Numerical Studies (No. 2014-01-2511). SAE Technical Paper.

[26] Yulin AV and Champneys AR. Discrete snaking: multiple cavity solitons in saturable media. SIAM Journal on Applied Dynamical Systems, 9(2) (2010), 391-431.

[27] Taylor C and Dawes JH. Snaking and isolas of localised states in bistable discrete lattices. Physics Letters A, 375(1) (2010), 14-22.

[28] Billah KY & Scanlan RH. Resonance, Tacoma Narrows bridge failure, and undergraduate physics textbooks. American Journal of Physics, 59(2), (1991), 118-124.

[29] MJ Rudd. Wheel/rail noise—part II: wheel squeal. Journal of Sound and Vibration, 46(3), (1976), 381-394.

[30] Nayfeh AH and Mook DT. Nonlinear oscillations. John Wiley & Sons, (2008).

[31] Nayfeh AH and Balachandran B. Applied nonlinear dynamics: analytical, computational and experimental methods. John Wiley & Sons, (2008).

[32] Samaranayake S and Bajaj AK. Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry. Journal of Sound and Vibration, 206(1) (1997), 39-60.