Nonlinear Dynamic Response of Ropeway Roller Batteries via an Asymptotic Approach

: The nonlinear dynamic features of compression roller batteries were investigated together with their nonlinear response to primary resonance excitation and to internal interactions between modes. Starting from a parametric nonlinear model based on a previously developed Lagrangian formulation, asymptotic treatment of the equations of motion was ﬁrst performed to characterize the nonlinearity of the lowest nonlinear normal modes of the system. They were found to be characterized by a softening nonlinearity associated with the stiffness terms. Subsequently, a direct time integration of the equations of motion was performed to compute the frequency response curves (FRCs) when the system is subjected to direct harmonic excitations causing the primary resonance of the lowest skew-symmetric mode shape. The method of multiple scales was then employed to study the bifurcation behavior and deliver closed-form expressions of the FRCs and of the loci of the fold bifurcation points, which provide the stability regions of the system. Furthermore, conditions for the onset of internal resonances between the lowest roller battery modes were found, and a 2:1 resonance between the third and ﬁrst modes of the system was investigated in the case of harmonic excitation having a frequency close to the ﬁrst mode and the third mode, respectively.


Introduction
Aerial ropeways are transportation systems which are becoming increasingly popular in recent years, not only in mountain regions with ski resorts and in sightseeing areas, but also in urban environments [1]. Two criteria can be used to classify ropeway systems, namely, the number of cables, and the operating mode. In this work, attention is focused on monocable ropeways, in which the vehicles are fixed on a single propelling cable. The cable travels continuously in a loop at a uniform speed and passes across the supporting towers, equipped with an ensemble of rollers, which facilitate a smooth cable transit across the tower where a change of slope occurs. Investigating the static and dynamic response of these multibody systems is challenging, and the literature on this topic is limited and mainly focused on the dynamic response of the vehicles under the ropeway operating condition [2,3]. Some works [4][5][6] investigated the dynamic behavior of carrying hauling ropes, dealing with the study of the effects of moving loads in an existing ropeway system, while [7] investigated the nonlinear coupling between the motion of the hauling cable and the swaying dynamics of the cabins in bi-cable circulating gondola ropeway systems. They show that these systems depend on the track inclination. A recent work [8] presented a computational approach for determining the initial tension of the carrying cable of a bi-cable ropeway under an in-service moving load, and to investigate the influence of nonlinear effects on the displacement state and axial forces. Furthermore, a theoretical model that solves the minimization of aerial ropeway vehicle oscillations, induced as the vehicle passes over a support, was proposed in [9], and then experimentally validated by the same authors in [10].
On the other hand, little attention was paid to the study of the local dynamic interactions taking place between the vehicle, the cable, the roller battery, and the supporting tower; few works have investigated these aspects [11]. A nonlinear model of compression roller batteries, used as a starting point for the studies proposed in the present work, was first developed in [12,13]. In such works, experimental tests were also carried out on an operating ropeway, so as to identify the modal characteristics of the system by adopting the Enhanced Frequency Domain Decomposition method, and to validate the analytical model by comparing the experimental dynamic responses with the numerical simulations. The model was further extended in [14,15] to include a passive vibration control system made by a group of linear or hysteretic vibration absorbers to mitigate the accelerations of roller batteries in cableways caused by vehicle transit.
To the best of the author's knowledge, there are no publications dealing with the characterization of the nonlinear dynamic behavior of the ropeways roller batteries or aimed at investigating their dynamic response to periodic excitations. The present work aims at filling the lack of knowledge regarding the characterization of the nonlinear dynamic features of roller battery systems and to investigate the effects of the periodic excitations induced by the moving cable that may give rise to resonance phenomena and nonlinear modal interactions. These phenomena were studied in the past in different mechanical systems, such as cranes [16][17][18][19], tethering [20], pendulum systems [21,22], piezoelectric beams [23,24], and plates [25,26], or, more generally, in slender structures possessing strong geometric nonlinearities [27][28][29]. Such nonlinearities can affect the dynamic response of highly deformable structural elements and can be suitably exploited for practical applications, as it was largely demonstrated in the literature [30][31][32][33]. Within this context, the asymptotic approach proved to be a suitable method to investigate the dynamics of nonlinear mechanical systems. This approach allows us to perform sensitivity analyses on the effect of the system parameters on the nonlinear dynamic response by determining closed-form expressions of the bifurcation response.
By employing the method of multiple scales [34][35][36], the asymptotic frequency response analyses were conducted in the present work, together with a deep discussion of the meticulously derived analytical treatment to highlight the scenarios of nonlinear interactions between the external excitation and the modes of the system. Furthermore, a thorough bifurcation analysis, carried out via path-following of the fixed points of the modulation equations shed light onto several interesting features of the roller battery singleand coupled-mode responses.

Nonlinear Dynamic Characterization
The proposed mechanical model of compression roller batteries is based on the formulation proposed, for the first time, in [13], which includes finite kinematics, for the description of the rollers and the vehicle motions, and obtains the nonlinear equations of motion via an energetic approach. In such a model, the mechanical system is represented by an assembly of four mechanical subsystems, namely, the hoisting beam supporting the roller battery, the roller battery, the elastic cable traveling across the ropeway and, finally, the cabin attached to the cable. The cable is modeled as a one-dimensional pre-stressed string. On the other hand, the roller battery is a mechanism made of eight rollers rigidly connected to each other by means of hierarchical balancers rotating around moving hinges and positioned symmetrically with respect to a connection point with the hoisting beam. The latter is modeled as an equivalent one-dimensional Euler-Bernoulli beam, while the cabin behaves like a pendulum with the support point fixed to the moving cable.
The whole modeling procedure, including all relevant details of the underlying complex mechanical behavior, such as the effects of the vehicle grip on the rollers or the periodic forcing caused by the interwire spacing on the outer layer of the monocable rope, was derived and extensively discussed in [13]. Therefore, the author of the present work will refer to the equations and mechanical parameters of the case-study roller battery reported in the above mentioned work. Since, in this paper, only the nonlinear dynamic behavior of the roller battery system is investigated and characterized, the presence of the cabin is neglected. Therefore, the mechanical formulation, reduced to the dynamic configuration considered in this work, is briefly summarized, and all the relevant equations are reported in the next sections.

Kinematic Descriptors
The kinematics are described in the plane containing the direction collinear with the cable configuration as well as the vertical direction. In particular, two fixed frames are introduced, namely, the frame e 0 1 , e 0 2 with origin in C and oriented as the horizontal and vertical directions, respectively. The frame (e 1 , e 2 ), centered in C and oriented along the roller battery in the inclined direction, thus rotated with respect to the fixed (global) frame by ϕ R (see Figure 1). The geometric parameters characterizing the reference configuration of the roller battery, the cable, the cabin, and the grip are shown in Figure 1. In particular, along the direction e 1 , d C denotes the distance between C and B k (k = 1, 2), d B is the distance between B k and A j , and d A represents the distance between A j and each of the two corresponding roller-cable contact points P i . Similarly, the distances h C , h B , and h A represent the distances along e 2 between the mentioned hinges. The parameters characterizing the reference geometry of the cable are the side lengths L S1 and L S2 (i.e., the distances of the two supports from the first and last rollers P 1 and P 8 , respectively). The roller battery degrees of freedom (DOFs) are the rotations θ A j (t), θ B k (t), and θ C (t), about the hinges A j , B k , and C, respectively. Therefore, the current configuration of the roller battery can be described by the vectors r P i (t), r A j (t), r B k (t), and r C (t) which provide, at time t, the positions of the cable contact points P i (i = 1, . . . , 8) and of the hinges A j (j = 1, . . . , 4), B k (k = 1, 2) and C, respectively, or by the displacement vectors of the above mentioned points, that is, u B k (t), u A j (t), and u P i (t), respectively.
To discretize the time-and space-dependent function v(s, t) representing the vertical displacement of the centerline hoisting cantilever beam, having length L hb , the Galerkin method is adopted by considering a suitable number of trial functions to approximate the dynamics of the beam-like subsystem. In particular, the lowest N b mode shapes the φ b (s) of the cantilever beam, which are employed to discretize the beam deflection in space, while q b (t) (b = 1, . . . , N b ) represent the time-dependent generalized coordinates, hence, the beam deflection can be written as v(s, t) ≈ ∑ On the other hand, due to the multiple contact points P i (i = 1, . . . , 8) between the cable and roller battery, the cable discretization is carried out according to the finite element (FE) technique. In particular, ne = 20 finite elements are used to discretize the cable length such that each contact point P i coincides with a finite element node while the cable support points are the boundary nodes. First-order Lagrangian polynomials are employed to approximate the cable transversal displacement and to describe the kinematic function w e (ζ, t) in the eth element. The latter can be written as w e (ζ, t) = Ne(ζ) Te X(t), where the subscript e indicates the element, ζ is the local arclength, Ne(ζ) is the 1 × 2 vector collecting the shape functions of each element, X(t) = [p 1 (t), . . . , p ne+1 (t)] is the vector of the nodal degrees of freedom, having size (ne + 1) × 1, and Te is the 2 × (ne + 1) extraction matrix. Finally, due to the cable fixed boundary conditions, it turns out that p 1 (t) = p ne+1 (t) = 0.
The unilateral contacts at points P i (i = 1, . . . , 8) between the rollers and the cable are modeled by introducing fictitious springs collocated at each point P i . These springs behave as internal elastic constraints whose stiffness K i is suitably tuned in order to simulate a quasi-rigid behavior.

Nonlinear Equations of Motion
The equations governing the motion of the roller battery system, including the interaction between the hoisting beam and the cable are derived via the Euler-Lagrange approach. To this end, the potential and kinetic energies of the beam and cable are first calculated together with the potential energy of the fictitious springs and the kinetic energy of the rollers to calculate the Lagrangian of the system, from which the nonlinear equations of motion of the roller battery are obtained. In particular, the potential energy of the hoisting beam and the cable can be expressed as: respectively, where ∂ s represents the partial differentiation with respect to the beam arclength s (with origin at the beam clamped cross section), EI hb is the hoisting beam flexural stiffness, while K c , incorporating the cable pretension N 0 , is the global stiffness matrix of the cable. The energy stored by the fictitious springs connecting the cable and the rollers is computed as follows: where ∆L i (i) = u P i (t) · e 2 − w e,i (t) is the relative displacement along the roller local direction e 2 (corresponding to the transverse direction of the cable) and w e,i (t) is the cable displacement at the node corresponding to the ith roller.
On the other hand, the hoisting beam and the cable kinetic energy contributions to the dynamics of the system are given by: respectively, where ∂ t and the overdot represent partial and total differentiation with respect to time t, respectively. ρA hb is the hoisting beam mass per unit length, while M c , includes the cable mass per unit length ρA c , is the global mass matrix of the cable. Finally, the kinetic energy possessed by the rollers can be calculated as: where m P i , m A j , and m B k are the equivalent masses of the roller battery lumped at P i , A j , B k , respectively, and m tb is part of the mass of the roller battery, acting at C. Finally, the stored and the kinetic energies of the overall system are given by: The simplifications made in the model formulation proposed in [13] include the assumption of considering the cable as a straight, prestressed string. Nevertheless, when crossing the roller battery system, the cable is subject to a small curvature due to the slight change in slop of the cable passing from the first to the last roller. Although this simplification does not imply relevant errors in the description of the dynamic behavior of this system (as demonstrated by the comparison with experimental results in [13]), the actual cable curvature generates a force acting on each roller transversally to the cable direction. Moreover, since the cable is made by an assembly of wires wound around its axis to study the effects of the interwire periodicity when the cable travels at the speed V c , harmonic forces orthogonal to the cable are considered to act on each roller. In particular, by denoting L per with the interwire distance of the outer cable layer, the periodic forces arising from the cable curvature will exhibit a circular frequency given by Ω = 2πV c /L per and can be expressed as F i (t) = F e 2 cos Ωt (i = 1, . . . , 8), assuming that the maximum amplitude F of the force is the same at each roller.
The motion equations are then obtained via a Lagrangian approach. In particular, by introducing the Rayleigh dissipation function, and by calculating the Lagrangian of the system L(t) = T (t) − U (t), the equations of motion can be written as: where l = 1, . . . , N s and N s = 7 + ne − 1 + N b is the total number of DOFs of the system, d l is the linear viscous damping coefficient associated with the lth DOF, d dt represents total differentiation with respect to time t, and x l is the lth generalized coordinate extracted out of the N s × 1 vector x(t) collecting the system degrees of freedom, namely, [13,14] for further details). Finally, in Equation (6), F l (t) represents the lth Lagrangian component of the periodic excitation F i (t) acting at the ith roller (i = 1, . . . , 8) and is calculated as:

Nondimensional Vector Form of the Equations of Motion
The equations of motion are then nondimensionalized by adopting the characteristic time t c = ρA hb L 4 hb /EI hb and the equivalent beam span L hb as the characteristic length, respectively. The latter can be used for the nondimensionalization of the systems geometric parameters, and for the unknown generalized coordinates corresponding to p 2 , . . . , p ne and q 1 , . . . , q N b , respectively. Furthermore, it turns out that f c = ρA hb L 2 hb /t 2 c can be used to nondimensionalize the forces acting in the system: in this case, the nondimensional maximum force acting on each roller can be calculated as λ = F/ f c . The numerical simulations are then carried out considering ne = 20 finite elements for the cable discretization and N b = 4 trial functions for the hoisting beam discretization; therefore, the system dynamics are described by N s = 30 DOFs. The geometric and mechanical properties of the roller battery considered in this work are taken in consonance with the data reported in [13].
The equations of motion reported in Equation (6) involve strongly nonlinear functions of the system DOFs (i.e., nonlinear combinations of trigonometric functions); nevertheless, their time integration can be quite easily performed by using well-established numerical integration schemes. On the other hand, their analytical treatment requires an asymptotic expansion of the nonlinear functions present in Equation (6). Since the system investigated includes only geometric nonlinearities, the most relevant effects of these nonlinearities are given by the approximating second-and third-order polynomial terms. Therefore, a third-order Taylor expansion of Equation (6) was performed, and the obtained system of nondimensional nonlinear equations of motion can be conveniently recast in vector form as: where x is the N s × 1 vector collecting the N s generalized coordinates, while M, C, and K are the N s × N s mass, viscous damping, and stiffness matrices of the system, respectively. Moreover, f 0 is the N s × 1 vector of the nondimensional amplitudes of the external excitation, while i o , d o , and g o (o = 2, 3), are the inertial, velocity-proportional, and stiffness N s × 1 nonlinear operators of oth order, respectively, arising from the system geometric nonlinearities. Finally, 0 is the N s × 1 null vector and the natural initial conditions are considered (i.e., x(0) = 0, andẋ(0) = 0, respectively).

Modal Characterization
Eigenvalue analysis is first performed to characterize the fundamental modes of the system. By linearizing Equation (7), i.e., setting , neglecting the forcing term and the linear viscous damping force (i.e., setting f 0 = 0 and C = O, being O the N s × N s null matrix), the homogeneous system of linear equations governing the free undamped oscillations of the system assumes the classical form: The N s eigenvalues of Equation (8) provide the natural frequencies ω s (s = 1, . . . , N s ) of the system and the corresponding eigenvectors, conveniently normalized, deliver the N s mode shapes of the roller battery. In particular, by referring to the N s × 1 vectors and natural frequencies of the mth and nth modes of the system, i.e., (ψ ψ ψ m , ω m ) and (ψ ψ ψ n , ω n ), respectively, the following relationships are satisfied: ψ ψ ψ m M ψ ψ ψ n = δ mn and ψ ψ ψ m K ψ ψ ψ n = ω 2 m δ mn , being δ mn the Kronecker delta. The lowest eight nondimensional circular frequencies of the system are reported in Table 1 and the corresponding mode shapes of the roller battery are displayed in Figure 2.  Interesting considerations emerge from the analysis of the ratios between the lowest frequencies of the system. Figure 3 reports the ratios between the frequencies of modes (2,3,4,5,6,7,8) with respect to the lowest frequency ω 1 (Figure 3a) and of the modes (3,4,5,6,7,8) with respect to the second frequency ω 2 (Figure 3b), respectively. As shown in Figure 3a, the frequency of the third roller battery mode is close to twice the frequency ω 1 of the lowest mode while the sixth and seventh modes have frequency three times ω 1 . The same integer ratios occur between the frequency of the seventh and eighth modes with respect to the frequency ω 2 of the second mode, i.e., 2:1 and 3:1 ratio, respectively, as shown in Figure 3b. When in a nonlinear dynamical system, the frequencies of two modes are in an integer ratio with each other, then the internal exchange of energy between the modes involved may occur, leading the system to dynamic interactions, i.e., the so-called internal resonances [34][35][36]. It is a matter of fact that the internal resonances may more easily be activated between modes having lower frequencies, such as in the case of the mechanical system investigated, shown by the first and third mode of the roller battery (see red bar and red line in Figure 3a). Therefore, when the frequency Ω of the cable-induced harmonic excitation is close to the frequency of the system, primary resonance may occur; furthermore, if the latter is a multiple or an integer ratio of the frequency of a different mode then internal resonance may also take place.
In the system investigated in this work, the cable-induced periodic excitation has a frequency dependent on the cable speed which never exceeds the nondimensional value of 2.8 × 10 −2 in the ropeway studied here; this value corresponds to a maximum nondimensional excitation frequency Ω = 7.5. This means that low cable speeds may generate periodic excitations inducing primary resonance in all of the lowest eight modes investigated. After characterizing the nonlinearity of the lowest modes of the system, the focus will be devoted to the analysis of the primary resonance of the lowest mode (i.e., Ω ≈ ω 1 ). Moreover, it will be shown that when Ω ≈ ω 1 , then 2:1 internal resonance with mode three may take place only for very high excitation amplitudes. On the other hand, the threshold for the activation of two-to-one internal resonance between mode three and mode one when the periodic excitation has frequency Ω ≈ ω 3 is calculated, and the nonlinear frequency response curves of the resonating modes are determined for a selected value of the excitation amplitude.

Asymptotic Solution of the Equations of Motion
The method of multiple scales [34][35][36] is adopted to first characterize the nonlinearity of the system modes and then to study the bifurcation scenarios due to direct harmonic excitations to obtain closed-form expressions of the stable and unstable branches of the FRCs, together with the continuation of the fold bifurcation points.
By letting denote a small nondimensional parameter representing the norm of the deviation of the current roller battery configuration from its reference state, a third-order expansion of the solutions of Equation (7) is sought in the form: where t 0 is the nondimensional fast time scale and t 2 is the slow time scale in which the nondimensional time derivative can be written as d/dt 2 2 , and ∂ 0 ∂ 2 = ∂ 2 /∂t 0 ∂t 2 , respectively. The perturbation treatment requires a preliminary rescaling of the viscous damping force and of the excitation term; therefore, since the nonlinearities of the system are such that resonant terms are generated at the cubic order. The linear viscous damping force is assumed to appear at the third order; this, analytically implies a re-scaling of the viscous damping matrix as C = 2 C, while the forcing term is rescaled so that f 0 = 3 f 0 . Substituting Equation (9) into Equation (7), and retaining only the terms up to order 3 , the following hierarchy of linear problems is obtained by equating to zero the coefficients of like powers of : order where D is a N s × N s matrix of constant terms, while g 1,2 , d 1,2 , i 1,2 , and i 2,1 are stiffness, velocity-proportional, and inertial third-order operators, respectively; which, in turn, are functions of the first-and second-order solutions.

Nonlinearity of the Lowest Normal Modes
By setting f 0 = 0 and C = O in Equation (12), it is possible to study the nonlinearity of the roller battery nonlinear normal modes by solving, in cascade, the hierarchy of perturbation problems given by Equations (10)- (12). In particular, to study the nonlinearity of the mth mode, the general solution of the -order problem, Equation (10), can be written as: where A m (t 2 ) andÃ m (t 2 ) are the complex and complex-conjugate amplitudes of the mth mode, respectively. The non-resonant terms are present in the second-order problem; therefore, after substituting the first-order solution (13) into Equation (11), the solution of the inhomogeneous 2 -order problem can be expressed in the form: where the vectorx 2,r (r = 1, 2, 3) is the linear combination of the N s modes of the system, that is:x 2,r = ∑ N s s=1 a r,s ψ ψ ψ s . The solution given by Equation (14) is substituted into the 2 -order problem and Equation (11) is then projected into the modal space through the N s -by-N s modal matrix Ψ, whose sth column is the vector ψ ψ ψ s (i.e., the sth mode), as: where is the vector collecting the secondorder inhomogeneous terms. By then equating the coefficients of terms having the same frequency in Equation (15), it is possible to calculate the 3N s coefficients a r,s and find the second-order solution x 2 .
First-and second-order solutions given by Equations (13) and (14), respectively, are then substituted into the 3 -order problem. The singularity generated by the secular terms exhibited in the third-order problem can be removed by enforcing solvability conditions [35]. To this end, by now introducing the general solution of the adjoint homogeneous problem of Equation (12): where cc stands for complex-conjugate, the solvability condition is enforced by requiring the orthogonality between (16) and the inhomogeneous term (right-hand side) in Equation (12) as: where: Equation (17), which annihilates the secular terms arising in r 3 , governs the modulation in time t 2 of the amplitudes A m andÃ m through to the following equations: where Γ m = Γ m (ω m ) is a strongly nonlinear function of the frequency ω m of the mth mode and of the systems mechanical parameters. After removing the secular terms, the solution of the inhomogeneous third-order problem can be calculated as: where the vectorx 3,r (r = 1, 2) is expressed, again, as a linear combination of the N s modes of the system, that is:x 3,r = ∑ N s s=1 b r,s ψ ψ ψ s . The solution given by Equation (20) is substituted into the 3 -order problem and Equation (12) is then projected into the modal space as: By then equating the coefficients of terms having the same frequency, it is possible to calculate the 2N s coefficients b r,s and find the second-order solution x 3 .
To study the nonlinearity of the mth normal mode, the complex amplitudes are expressed in polar form as A m ( , respectively, where a m (t 2 ) is the real amplitude and θ m (t 2 ) is the phase of the mth nonlinear normal mode. The polar form is then substituted into the modulation equations (19). Furthermore, the separation of the real and imaginary parts yields the following differential relationships in terms of the real amplitude and the phase of the mth mode: By defining the relative phase γ m (t 2 ) = σt 2 − θ m (t 2 ), with σ being a detuning parameter, it turns out that the fixed points of the dynamic system are obtained when ∂ 2 a m = 0 and ∂ 2 γ m = 0; therefore, the detuning can be calculated from Equation (21) as σ = Γ m a 2 m and the frequency of the mth nonlinear normal mode can be expressed in terms of the real modal amplitude a m as: The function Γ m is also known as the effective nonlinearity coefficient of the mth mode: when Γ m > 0, the nonlinearity is hardening, while, when Γ m < 0, it is softening.

One-Mode Projection of the Nonlinear Equations of Motion
The nonlinearity of the mth normal mode can be qualitatively investigated by neglecting the contribution of the remaining N s − 1 modes in the reconstruction of the second-and third-order solutions given by Equations (14) and (20). This approach, although simplified, creates a better understanding of the role of the system mechanical parameters into the effects of the geometric nonlinearity of the system. To this end, the solution of Equation (7) is expressed as x = ψ ψ ψ m q m (t), where q m (t) is the mth generalized coordinate. Equation (7) is then projected into the mth mode by pre-multiplying it by ψ ψ ψ m to obtain the mth nonlinear modal equation in the form: where c m = ψ ψ ψ m C ψ ψ ψ m is the modal damping coefficient, f m = ψ ψ ψ m f 0 is the modal component of the excitation term, and: are the second-and third-order inertial, velocity-proportional, and stiffness modal coefficients, respectively. The same asymptotic procedure, based on the method of multiple scales, can be used to investigate the nonlinearity of the mth mode through the nonlinear equation of motion (23), whose solution is sought in the form q m (t) = q m,1 (t 0 , t 2 ) + 2 q m,2 (t 0 , t 2 ) + 3 q m,3 (t 0 , t 2 ). Therefore, after rescaling the excitation term as f m = 3 f m and the modal viscous damping force as c mqm = 2 c mqm , the lowest three perturbation problems can be written as: To study the nonlinearity of the mth mode the excitation term is set to zero, i.e., f m = 0, and the modal viscous damping is neglected, i.e., c m = 0.
The solution of the first-order problem can be written as q m,1 = A m (t 2 )e iω m t 0 + A m (t 2 )e −iω m t 0 , while the solution of the second-order problem, which does not show resonant terms, is given in the form: At the third order, the singularity generated by the resonant terms is removed when the modulation equations, having the same form of Equation (19), hold; although, in this case, the nonlinearity coefficient of the mth mode has the simpler expression: After removing the secular terms by superimposing Equation (19), the solution of the inhomogeneous third-order problem can be written as: and can be substituted into the inhomogeneous third-order problem; then, by equating the coefficients of terms having the same frequency, it is possible to calculate the unknown coefficients of Equation (28), which turn out to beq 3,1 =q 3,2 =q 3 (ω m ), where: To study the nonlinearity of the mth normal mode and express the complex amplitudes in polar form, substitute their expressions into the modulation equations and separate the real and imaginary parts; this furnishes differential relationships in terms of the real amplitude and the phase of the mth mode in the same form of Equation (21). By introducing the relative phase γ m (t 2 ), the fixed points are obtained when ∂ 2 a m = 0 and ∂ 2 γ m = 0; hence, the frequency of the mth nonlinear normal mode can be expressed in terms of the real modal amplitude a m and through the effective nonlinearity coefficient given by Equation (27) as: ω nl m = ω m + Γ m a 2 m . The backbone curves of the lowest four nonlinear normal modes, that is, the nonlinear relationships between the modal amplitude a m and the detuning σ (where σ = ω nl m − ω m ), are shown in Figure 4a. The dashed lines in Figure 4a indicate the backbones calculated by adopting the one-mode projection, while solid lines refer to the solution obtained by using the full modal basis discretization. As shown in the figure, the differences between the curves obtained with the two approaches are more evident for the lowest skewsymmetric mode (red lines) and the lowest symmetric mode (blue lines). Nevertheless, as shown in Figure 4b,c, the shapes of the corresponding modes are less influenced by the modal discretization. Figure 5 reports the values of the coefficients of the nonlinear terms of Equation (23) for the lowest four modes. It is shown that all coefficients of the quadratic terms are negative. On the other hand, the inertial and velocity-proportional coefficients of the cubic terms are positive for all modes, while the coefficient of the stiffness cubic term is always negative. Moreover, the coefficients of the stiffness terms are orders of magnitude higher than the inertial and velocity-proportional coefficients. This indicates that the nonlinearity of the roller battery is mainly governed by the stiffness terms while the nonlinear effects provided by the inertial and velocity-proportional terms are negligible. This is further confirmed by the values of the nonlinearity coefficients Γ m reported in Table 2. In the latter are the values of Γ m of the lowest four modes calculated by considering the contribution of: (i) all types of nonlinearity, (ii) the inertia-associated nonlinearity (i.e., by setting c d2 = c d3 = c k2 = c k3 = 0), (iii) the velocity-associated nonlinearity (i.e., by setting c i2 = c i3 = c k2 = c k3 = 0), (iv) the stiffness-associated nonlinearity (i.e., by setting c i2 = c i3 = c d2 = c d3 = 0), respectively.
The results reported in Table 2, clearly show that the stiffness-related nonlinearity is predominant since the effective nonlinearity coefficient Γ m is an order of magnitude higher with respect to the associated nonlinearity of the inertia and velocity-proportional terms. On the other hand, it is shown that the velocity-proportional nonlinear terms confer to the system a hardening nonlinear behavior (i.e., Γ m > 0).

Cable Motion-Induced Resonances
In this section the nonlinear dynamic response of roller batteries subjected to periodic excitation provided by the forces induced by the moving cable will be investigated. In particular, the focus is first devoted to the analysis of the primary resonance of the lowest mode (i.e., when Ω ≈ ω 1 ); then, the two-to-one internal resonance between the first and third mode will be investigated when Ω ≈ ω 1 and when Ω ≈ ω 3 , respectively. To this end, the asymptotic approach presented in Section 4 based on the full modal discretization for the evaluation of the second-and third-order solutions (see Section 4.1), is adopted to calculate the nonlinear frequency response curves and the stability regions of the dynamic responses.

Primary Resonance
First-and second-order solutions of the hierarchical systems of equations, Equations (10) to (12), are given by Equations (13) and (14), respectively. The latter are then substituted into the 3 -order problem Equation (12) and the singularity generated by the secular terms can be removed by enforcing the solvability condition. To this end, after introducing the general solution of the adjoint homogeneous problem (i.e., Equation (16)), the orthogonality between (16) and the inhomogeneous term (right-end side) in Equation (12) is enforced as: To investigate the primary resonance of the mth mode, the frequency of the excitation term is expressed as Ω = ω m + 2 σ, with σ as a detuning parameter. Moreover, to make use of complex algebra, cos Ωt = 1 2 e i ω m t 0 +i σ t 2 + cc; then, at the third order, the singularity generated by the resonant terms is removed when the following modulation equations hold: where, C m = C m (ω m ) and F 0 = F 0 (λ, ω m ) are strongly nonlinear functions of the linear frequency of the mth mode. Note that, in case of the one-mode projection of the equations of motion, it turns out that C m ≡ c m , F 0 ≡ f m , and Γ m is given by Equation (27). By expressing the complex amplitudes in polar form and separating the real and imaginary parts, the modulation equations in terms of the real amplitude and phase read: By solving Equation (32) for ∂ 2 a m and ∂ 2 θ m and introducing the relative phase γ m (t 2 ) = σ t 2 − θ m (t 2 ), the fixed points can be sought to impose ∂ 2 a m = 0 and ∂ 2 γ m = 0, respectively. By solving the ensuing equations in terms of cos γ m and sin γ m , and enforcing the trigonometric identity, the frequency-response function providing the nonlinear relationship between the real amplitude a m and frequency detuning σ can be obtained in its implicit form as: The equation that allows for describing the bifurcation behavior of the system is given by G(a m , σ) = 0; moreover, from the stationarity of G(a m , σ) with respect to a m it is possible to calculate the loci of the fold bifurcation points as the four non-trivial solutions of the equation ∂ a m G(a m , σ) = 0, that is: with a F m (1) being the trivial solution (i.e., a F m (1) = 0). Finally, after removing the secular terms by means of Equation (31), the solution of the inhomogeneous third-order problem can be calculated as 3 m e −3iω m t 0 , where the vectorsx 3,r (r = 1, 2) are expressed, again, as a linear combination of the N s modes of the system, that is:x 3,r = ∑ N s s=1 b r,s ψ ψ ψ s . The third-order solution is then substituted into the 3 -order problem and it is then projected into the modal space as: By then equating the coefficients of terms having the same frequency, it is possible to calculate the 2N s coefficients b r,s and find the second-order solution x 3 .
Thereafter, the nonlinear equations of motion of the system expressed by Equation (6) are numerically integrated in time by using the built-in NDSolve algorithm present in the software Mathematica [37]. For each simulation, at a selected excitation frequency Ω, the maximum integration time was set to t max = 500(2π/Ω) , where (2π/Ω) is the excitation period, while the time step was ∆t = (2π/Ω)/40. At steady-state, the state vector (x, x) was saved and assigned as the initial condition for the time integration run with the subsequent Ω. Finally, backward and forward frequency sweeps were performed so as to calculate the numerical FRCs. Figure 6 shows the frequency response curves in terms of the vertical component (i.e., along the direction e 2 ) of the nondimensional displacement vectors u of points P 8 , A 2 , B 1 , and C, respectively, when Ω ≈ ω 1 . The comparison between numerical time integration (gray dots) and the third-order asymptotic solution obtained via the method of multiple scales (solid and dashed lines) shows a very good agreement in most of the range of stable solutions. The asymptotic analysis allowed us to investigate the unstable solutions (dashed lines) and calculate the threshold values of amplitude and frequency, where fold bifurcation occurs (i.e., at the boundary between solid and dashed lines). The analysis was performed for a value of the nondimensional maximum amplitude of the periodic force λ = 8.27 × 10 −3 , so as to show the softnening nonlinearity of the dynamic response. The value of the forcing amplitude is almost two orders of magnitude lower than the characteristic force of the system. This means that moderately high values of the forcing term can also induce a nonlinear resonance response of the system.  Figure 7a shows the FRCs in terms of the modal amplitude a m for two values of the modal damping coefficient C m for the case m = 1. In particular, C m was rescaled through a nondimensional coefficient ξ as C m = ξ C m and the cases of ξ = 1 and ξ = 0.5 are considered. This rescaling allowed us to perform a parametric study on the effect of damping on the nonlinear resonance response. In Figure 7a we report the regions collecting the unstable solutions in the plane (Ω, a m ) (i.e., blue and black shaded surfaces). The latter are bounded by the curves enveloping the fold bifurcation points and provided by Equation (34). The loci of the fold bifurcation points are then plotted in Figure 7b for

Two-to-One Internal Resonance
In this section, we investigate the phenomenon of the two-to-one internal resonance between two modes of the roller battery occurring when one mode of the system undergoes primary resonance due to the periodic excitation of the moving cable.
Without loss of generalization, the two resonant modes are indicated as mode m and mode n, respectively, and the corresponding frequencies have a ratio close to two (i.e., ω n ≈ 2ω m ). Two case-studies are investigated, namely, the case when the frequency Ω of the periodic excitation is near the frequency of the mode m (i.e., primary resonance of mode m, with Ω ≈ ω m ), and the case when Ω ≈ ω n (i.e., primary resonance of mode n), respectively. As mentioned in Section 3, 2:1 internal resonance is possible in the mechanical system studied in this work since the frequencies of the third and first mode have a ratio close to two.
In this case, the resonant terms are generated at the quadratic order, therefore, the viscous damping force and excitation need to be rescaled so as to appear at the second perturbation order; hence, we set C = C and f 0 = 2 f 0 , respectively. By now substituting Equation (35) into Equation (7), and retaining only terms up to order 2 , the following hierarchy of linear problems is obtained by equating to zero the coefficients of like powers of : order In the case of two interacting modes (i.e., mode m and mode n), the general solution of the -order problem (36) can be written as: where A m (t 1 ) and A n (t 1 ) are the complex amplitudes of the mth and the nth mode, respectively, while cc stands for complex conjugate. The first-order solution to Equation (38) is then substituted into the 2 -order problem. The singularity generated by the secular terms exhibited in the second-order problem can be removed by enforcing solvability conditions. By introducing the mth and nth general solutions of the adjoint homogeneous problem of Equation (37) we get: x * n = i ω n ψ ψ ψ n e −i ω n t 0 + cc, the solvability condition is enforced by requiring the orthogonality between (39) and (40) and the inhomogeneous term (right-end side) in Equation (37) as: To express the closeness of the two interacting frequencies, the detuning parameter σ 2 is introduced so that ω n = 2ω m + σ 2 . On the other hand, the primary resonance condition between the periodic force and the modes m and n is achieved when Ω = ω m + σ 1 and Ω = ω n + σ 1 , respectively, with σ 1 being a further detuning parameter. Thus, substituting the external and internal resonance conditions into Equation (41), it is possible to obtain the modulation equations for the amplitudes A m and A n (and their complex conjugates) for both cases of Ω ≈ ω m and Ω ≈ ω n , respectively.
Finally, after removing the secular terms, the solution of the second-order inhomogeneous problem can be calculated as: m e 2iω m t 0 +x 2,2 A 2 n e 2iω n t 0 +x 2,3 A m A n e iω m t 0 +iω n t 0 +x 2,4Ãm A n e −iω m t 0 +iω n t 0 +x 2,5 A m e iω m t 0 +x 2,6 A n e iω n t 0 +x 2,7 A mÃm +x 2,8 A nÃn +x 2,9 e iΩt 0 + cc, where the vectorx 2,r is expressed as a linear combination of the N s modes of the system, that is:x 2,r = ∑ N s s=1 a r,s ψ ψ ψ s . Note that, the terms proportional tox 2,4 ,x 2,5 ,x 2,6 ,x 2,9 , and their complex conjugates, respectively, are orthogonal to the general solutions (39) and (40) of the adjoint homogeneous problem of Equation (37). Therefore, although present inr 2 , these terms vanish when the solvability condition given by Equation (41) is enforced. The second-order solution is substituted into the 2 -order problem which is then projected into the modal space as: Ψ M ∂ 2 0 x 2 + Ψ K x 2 = Ψ r 2 . Finally, by equating the coefficients of terms having the same frequency, it is possible to calculate the coefficients a r,s of the second-order solution.

The Case of Ω ≈ ω m
In the case of primary resonance on mode m (i.e., Ω ≈ ω m ), the modulation equations of the amplitudes A m and A n have the following form: where α 1 and α 2 are positive coefficients, while C m >, C n > 0 are the components of the damping coefficients in the modulation equation of A m and A n , respectively, while F m > 0 is the component of the excitation amplitude in the modulation equation of A m . By expressing the complex amplitudes in polar form as A m = 1 2 a m e i θ m , A n = 1 2 a n e i θ n , introducing the relative phases γ 1 = σ 1 t 1 − θ m and γ 2 = σ 2 t 1 − 2θ m + θ n , and substituting the polar form into Equation (43). Separating the real and imaginary parts yields the four differential equations governing the slow modulation of the real amplitudes and the relative phases. Then, the fixed points can be calculated by setting ∂ 1 a m = 0 = ∂ 1 a n and ∂ 1 γ 1 = 0 = ∂ 1 γ 2 , respectively, and the bifurcation equations can be written as: C m a m ω m − 1 2 α 1 + α 2 ω 2 m a m a n sin γ 2 − 2F m sin γ 1 = 0, C n a n ω n + 1 4 α 1 + 1 4 α 2 ω 2 n a 2 m sin γ 2 = 0, σ 1 a m ω m + 1 4 α 1 + α 2 ω 2 m a m a n cos γ 2 + F m cos γ 1 = 0, σ 2 a m a n ω m ω n + 1 2 α 1 + α 2 ω 2 m ω n a m a 2 n − 1 8 α 1 + 1 4 α 2 ω 2 n ω m a 3 m cos γ 2 + 2F m a n ω n cos γ 1 = 0. (44) The system of four nonlinear algebraic equations (44) can be solved for selected values of the detuning parameters, σ 1 and σ 2 , and in terms of the nondimensional forcing amplitude λ (embedded in the coefficient F m ) to obtain the force response curves shown in Figure 8a. It is interesting to note that, in the case of primary resonance on mode m (i.e., Ω ≈ ω m ), a 2:1 internal resonance involving large amplitude oscillations of mode n arises for values of the nondimensional excitation amplitude λ orders of magnitude higher with respect to the value of the amplitude used in Section 5.1 (i.e., λ = 8.27 × 10 −3 ), to investigate the primary resonance of the mth mode. As shown in Figure 8b, only at very large values of λ is it possible to obtain frequency response curves involving comparable amplitudes of oscillation of the two resonating modes, a m and a n , respectively (i.e., coupledmode oscillations). In this circumstance, the modal amplitudes are two orders of magnitude higher, with respect to the case of primary resonance of mode m with no internal resonance. In the case considered (i.e., λ = 2.5), the FRCs of amplitudes a m and a n have a specular behavior with respect to the detuning σ 1 (i.e., the curves are symmetric with respect to the axis σ 1 = 0). Moreover, stable (solid lines) and unstable (dashed lines) responses are obtained in a short range of the detuning parameter around the value σ 1 = 0, while, in a larger range of positive and negative values of the detuning parameter σ 1 , multistable responses are found. Here the stability of the dynamic response is ascertained by the study of the eigenvalues of the Jacobian matrix of Equation (44). On the other hand, in the range of low values of the excitation amplitude λ, internal resonance induces negligible oscillation amplitudes of mode n, as shown in Figure 9a. In particular, at λ = 8.27 × 10 −3 (dashed vertical line) a n is more than one order of magnitude lower than a m , as shown in Figure 9b where it is reported that the ratio between the amplitude a m of the mode directly excited, and the amplitude a n of the internally resonant mode. The study conducted allows the claim that, although the third and first mode of the roller battery investigated are internally resonant, in the case of periodic excitation of the first mode, the amplitude of the forcing term must be enormously high to activate appreciable oscillations of the third mode. For the sake of rigor, it is worth noting that the ratio between the frequency of the third and first mode of the roller battery is not exactly equal to two. In particular, ω m = ω 1 = 0.376 and ω n = ω 3 = 0.734; this implies that, in the case-study investigated σ 2 = −0.018 (i.e., ω 3 = 2ω 1 − 0.018). Therefore, the force response curves and the FRCs are slightly different with respect to the perfectly resonant case (σ 2 = 0), although the dynamic behavior does not change appreciably. As shown in Figure 10a, the range of λ, inducing comparable amplitudes a m and a n in the coupled-mode oscillations, is the same of the case when σ 2 = 0, while in Figure 10b it is shown how the nontrivial detuning σ 2 implies a loss of symmetry in the FRCs, although the nonlinear behavior is qualitatively the same as the case when σ 2 = 0.
The force response curves shown in Figure 11a can be calculated by solving the system of four nonlinear algebraic equations (46) for selected values of the detuning parameters σ 1 and σ 2 and in terms of the nondimensional forcing amplitude λ (embedded in the coefficient F n ). Solid and dashed lines indicate stable and unstable responses, respectively, and the stability is ascertained by monitoring the eigenvalues of the Jacobian matrix of the Equation (46). Figure 11a displays a paradigmatic behaviour in the context of 2:1 internal resonances for a variety of nonlinear systems, the so-called saturation phenomenon. The response of the directly excited mode grows linearly up to a threshold excitation amplitude λ * past which the single-mode response becomes unstable in favour of a two-mode response. The directly excited response (i.e., the amplitude a n ) saturates at a constant value, while the amplitude of the internally resonant mode (i.e., the amplitude a m ) grows nonlinearly with the excitation amplitude. In Figure 11b,c the bifurcation scenarios are shown through the FRCs calculated at two different values of the excitation amplitude λ (i.e., at λ = 0.21 and λ = 8.27 × 10 −3 , respectively). These are obtained by path-following the fixed points of the modulation equations (46) and the stability of the system is monitored by checking the eigenvalues of the Jacobian matrix. Figure 11c, obtained for the value of λ < λ * , shows the FRC of the single-mode response (i.e., a m = 0, a n = 0) in which the only nontrivial amplitude is that of the mode directly excited by the periodic force (gray line).
The force response curves and the frequency response curves depicted in Figure 11 are obtained for σ 1 = 0 and for a perfect tuning among the frequencies of the resonant modes, i.e., σ 2 = 0, while Figure 12a shows the force response curves in the case of negative detuning (i.e., σ 2 = −0.018), which characterizes the system investigated in this work. As seen in the case illustrated in the previous section (i.e., Ω ≈ ω m ), when the primary resonance involves the higher frequency mode (i.e., Ω ≈ ω n ), the frequency response curves are symmetric with respect to the detuning parameter σ 1 , while losing their symmetry when σ 2 = 0 (as shown in Figure 12b). From the bifurcation analyses conducted for both cases (σ 2 = 0 and σ 2 = −0.018, respectively) it turns out that, as expected [34], at point A and C (A and C ), the single-mode response (a m = 0, a n = 0) becomes unstable. Fold bifurcation occurs at points B and B , respectively; between those two points, the two-mode response (a m = 0, a n = 0) shows multistability. Finally, points D and D represent the boundaries of the detuning parameter σ 1 , outside which, the response is always a single-mode stable solution.  (a) Force response curves of the amplitude a m (black line) and a n (gray line) when σ 1 = 0 and σ 2 = −0.018. (b) Frequency response curves of the amplitude a m (black line) and a n (gray line) when σ 2 = 0 and at high excitation amplitude λ = 0.21.

Conclusions
In this work, the nonlinear dynamic behavior of roller batteries was characterized, investigated, and discussed for the first time in the literature. An accurate analytical treatment of the equations of motion was proposed to highlight the nonlinearities of the system and the method of multiple scales was employed to deliver closed-form expressions of the roller battery nonlinear dynamic responses. The nonlinearity of the lowest normal modes of the system was clearly identified and showed to be of the softening type. The simplification of the asymptotic treatment, obtained through the one-mode projection of the equations of motion, shed light onto the effects of the system nonlinearities and to prove that the stiffness nonlinearities govern the dynamic behavior of the roller battery.
The effects of the periodic direct excitation deriving from the cable motion were investigated by studying the primary resonance of the lowest normal mode of the roller battery. It was shown that, for moderately high values of the excitation, the nonlinear response was multistable in a wide range of excitation frequencies. The analytical results were also validated through numerical time integration of the nonlinear equations of motion. Moreover, the instability regions, whose boundaries were analytically obtained in closed form, were parametrically investigated so as to discuss the role of the structural damping, and for demonstrating that higher values of the damping coefficient reduce the region where unstable responses occur.
Furthermore, the modal characterization of the system allowed us to discover the potential scenarios of modal interactions and to investigate the phenomenon of internal resonance between two modes having frequency ratios close to two (i.e., the third and first mode) both in the case of primary resonance of the lowest mode and of the higher mode, respectively. A stability analysis of the fixed points of the modulation equations for the two mode amplitudes and phases allowed us to describe the bifurcation scenarios both for the single-mode and the coupled-mode responses. It was shown that, although internal resonances were possible in the system, in order to cause large-amplitude coupled-mode oscillations, the amplitude of the periodic excitation must be orders of magnitude higher and sufficient to induce large oscillations in the primary resonance of the lowest mode.
Funding: This research received no external funding.