Computation of quasi-periodic localised vibrations in nonlinear cyclic and symmetric structures using harmonic balance methods

In this paper we develop a fully numerical approach to compute quasi-periodic vibrations bifurcating from nonlinear periodic states in cyclic and symmetric structures. The focus is on localisedoscillations arisingfrommodulationallyunstabletravellingwavesinducedby strong external excitations. The computational strategy is based on the periodic and quasi-periodic harmonicbalancemethodstogetherwithan arc-lengthcontinuation scheme.Dueto the pres- ence of multiple localised states, a new method to switch from periodic to quasi-periodic states is proposed. The algorithm is applied to two different minimal models for bladed disks vibrating in large amplitudes regimes. In the ﬁrst case, each sector of the bladed disk is mod-elled by a single degree of freedom, while in the second application a second degree of free- dom is included to account for the disk inertia. In both cases the algorithm has identiﬁed and tracked multiple quasi-periodic localised states travelling around the structure in the form of dissipative solitons.


Introduction
The emergence of localised vibrations in cyclic and symmetric structures is an important phenomenon due to potential problems induced by high cycle fatigue. Mechanical components with cyclic symmetry are very common in the aerospace industry, where e.g. bladed disks of aircraft engines [1], space antennas [2] and reflectors [3] are usually composed of ideally identical substructures assembled in a cyclic and symmetric configuration. In the linear framework localised vibrations arise due to the lack of symmetry induced e.g. by manufacturing variability or wear. This topic is well-known in the aerospace industry as a mistuning problem, and studies have mainly focused on effective modelling techniques [4,5], experimental characterisation together with model identification [6,7], and even the use of intentional mistuning to reduce the level of vibrations in real applications [8,9].
In real applications, structures may deviate from the linear behaviour e.g. due to frictional interfaces, large displacements, or fluid-structure interaction [1,[10][11][12][13]. In the presence of nonlinearities the emergence of localised vibrations might go much beyond mistuning. It is well-known, for example, that even perfect cyclic and symmetric structures can localise vibrations in the nonlinear regime due to the dependence of mode shapes on amplitude, and also due to bifurcations [14][15][16][17][18][19][20]. Moreover, the effect of inhomogeneities when a system vibrates in the nonlinear regime is a challenging research field, and the role played by nonlinearity in an inhomogeneous mechanical system experiencing energy localisation is an active topic [21,22]. This paper introduces a fully numerical approach to compute quasi-periodic states bifurcating from periodic vibrations due to nonlinearities. The focus is on localised states resulting from the emergence of modulation instability for travelling wave excitations in cyclic and symmetric structures. This investigation has been motivated by recent findings in Refs. [23,24], in which the authors show that localised states in the form of solitons bifurcate from travelling wave responses. The previous findings have been computed using an asymptotic method applied to a very idealised minimal model, and more flexible methods are required for the case of more complex and more realistic applications. The fully numerical approach proposed in this paper relies on harmonic balance methods (HBMs) for the computation of periodic and quasi-periodic vibrations. Moreover, due to the presence of several coexistent localised states, a method to switch from periodic to quasi-periodic branches near bifurcation points is introduced. The numerical approach is applied to two different minimal models: (1) a bladed disk with one degree of freedom per sector; and (2) a bladed disk with two degrees of freedom per sector. The computed results show that, in both cases, strong localised states are possible due to nonlinear envelope dynamics.
The paper is organized as follows. In Sec. 2 the fundamentals of the HBM is introduced, including a new technique to switch from periodic to quasi-periodic branches. Section 3 focuses on the investigation of two minimal models, and quasi-periodic localised solutions are computed using the approach described in Sec. 2. Finally, Sec. 4 describes the main conclusions and suggests directions for future investigations.

The harmonic balance method
We investigate nonlinear dynamical systems mathematically described by where M, C, and K are the n × n mass, damping, and stiffness matrices, respectively, whileẍ(t),ẋ(t), and x(t) represent the corresponding acceleration, velocity, and displacement vectors. In Eq. (1), f ext (t) is the vector of external excitations, while the vector f nl(x) models the displacement dependent nonlinear forces.
We assume that the system in Eq. (1) is vibrating under periodic or quasi-periodic 1 regimes. Therefore, x(t) is expanded in a two-frequencies quasi-periodic Fourier series as In Eq (2), h 1 and h 2 are the numbers of retained harmonics for the frequencies 1 and 2 , respectively, a 0,0 is the vector of static displacements, while a k 1 ,k 2 and b k 1 ,k 2 are the vectors representing the Fourier coefficients. It is important to note that when h 1 or h 2 is null the quasi-periodic expansion in Eq. (2) recovers the standard periodic HBM (see e.g. Refs. [10,25,26]). Moreover, the quasi-periodic expansion is only valid when the two frequencies 1 and 2 are incommensurable, otherwise 2 can be written as a harmonic of 1 and, consequently, x(t) is in the periodic regime (see e.g. Ref. [27]).
In the following we introduce a hyper-time domain with two independent variables 1 = 1 t and 2 = 2 t. The scalar product between two quasi-periodic functions f 1 ( 1 , 2 ) and f 2 ( 1 , 2 ) in hyper-time domain is defined as The aim of the HBM is to transform the set of nonlinear differential equations in Eq. (1) into a system of nonlinear algebraic equations. This system is obtained after substituting Eq. (2) into Eq. (1), and finally projecting the resulting equation onto the trigonometric basis using Eq. (3). This process leads to where L( 1 , 2 ) is the matrix representing the linear dynamics, z is the vector of unknown Fourier coefficients in Eq. (2), and 0 is the null vector. The matrix L( 1 , 2 ) is algebraically represented as [27][28][29] where the inner matrices L 0,0 and L k 1 ,k 2 are In Eq. (4), g nl (z) and g ext are the vectors projecting f nl (x(t)) and f ext (t) into the frequency domain, respectively. Usually, g nl (z) is not directly computed since displacement dependent nonlinear forces are, most of the times, defined in the time domain. Therefore, an alternating frequency-time (AFT) domain technique based on the hyper-time concept [27,28,30], as depicted in Fig. 1, is implemented. In order to illustrate this approach, suppose that the set of coefficients a 00 , a k 1 ,k 2 and b k 1 ,k 2 in Eq. (2) are given. Therefore, the displacement field x( 1 , 2 ) and the corresponding nonlinear forces f nl (x( 1 , 2 )) define surfaces in hyper-time domain. Finally, one can compute the Fourier coefficients of the nonlinear forces g nl (z) using the scalar product in Eq. (3). However, in practice, g nl (z) is obtained directly using a two-dimensional fast Fourier transform (2D-FFT). Similarly, the transformation from frequency to time domain can be implemented by means of an inverse 2D-FFT.

Harmonic selection
The number of harmonics defined by k 1 and k 2 in Eq. (2) can be drastically reduced. Without any additional information regarding the nonlinear forces f nl(x) , the vector of unknowns z can be reduced to roughly half of its original dimension due to the symmetry properties (odd/even) of the harmonic functions. In the literature, several different techniques have been suggested to decrease even more the dimension of z and, consequently, the computational efforts [27,31]. In this paper we are interested in localised vibrations induced by strong travelling wave responses. It is well-known in the literature that localised solutions bifurcate from travelling wave states due to modulation instability [32][33][34][35]. In this case, any general physical coordinate x n (t) can be written as where ( 2 t) is the modulating (or envelope) function dependent on 2 only, and f cw ( 1 t) is the carrier wave. One should note that in the case of a carrier wave with no static component the spectrum of x n (t), i.e. its frequency contain, cannot depend on 2 only. Therefore, the harmonics referred to k 1 = 1 can be removed from the solution basis. Fig. 2 displays two examples of retained harmonics using the full basis and the reduced one. In Panel (a), for h 1 = h 2 = 4, the number of retained harmonics reduces from 41, using the full harmonic expansion, to 36, using the modulation instability property. The reduction is more efficient when h 2 ≫ h 1 , as in Panel (b). In this case, for h 2 = 7 and h 1 = 1, the number of retained harmonics reduces from 23, using the full harmonic basis, to 15, when solutions localise due to modulation instability.

Stability of periodic solutions
Stability of periodic solutions is addressed using Floquet theory. Thus, a small perturbation h(t) is added to an already known periodic solution x p (t), and linear stability of x p (t) is computed based on the evolution of h(t). In this case, three different configurations are possible: and (3) if h(t) is constant then the system is marginally stable and a nonlinear stability analysis is required. It is possible to demonstrate (see e.g. Ref. [36]) that the linear stability of x p (t) can be computed from the so-called monodromy matrix . The most straightforward way to compute consists in time integrating the linearised system with periodic coefficientṡ where y(t) = {x(t)ẋ(t)} T , and This integration is carried out over one period T using 2n linearly independent initial conditions. For simplicity, the ith initial condition is usually assumed as The three possible scenarios about the stability of x p are given by the eigenvalues i of : (1) if max(| i |) < 1, the system is linearly stable; if max(| i |) > 1, the system is linearly unstable; and in the special case (3) when max(| i |) = 1, the system is marginally stable. In this paper, we are interested in unstable Floquet multipliers i with imaginary parts different than zero (Im( i ) ≠ 0). In this case, the system experiences a Neimark-Sacker (NS) bifurcation, and a branch of quasi-periodic solutions detaches from periodic states.

Numerical continuation
In the case of quasi-periodic oscillations arising due to NS bifurcations the second frequency 2 has to be treated as an unknown. This unknown can be added to the algebraic system by means of a phase condition. Physically, this phase condition pins down the envelope function by removing the intrinsic rotating symmetry of the problem. Mathematically, this constraint  can be included in Eq. (1) by setting one of the sinusoidal terms for the harmonic k 1 = 1 and k 2 = 1 as equal to zero. Therefore the new nonlinear algebraic system, with 2 as un unknown, is written as Fig. 3. Strategy to switch from periodic to quasi-periodic branches. In the first step, the NS bifurcation is identified from the Floquet multipliers. In the second step the unstable mode is used as initial conditions for Eq. (8) and the growth rate from the computed response is removed. The final result is added to the corresponding periodic regime, in the third step, and the final configuration in frequency domain is used as an initial guess for the solution of Eq. (11).
where a i 1,1 states the k 1 = k 2 = 1 cosine term for the ith degree of freedom. One should note that the harmonics k 1 = k 2 = 1 and the cosine term (instead of the sine term) have been chosen arbitrarily. In general, any non-zero quasi-periodic sine or cosine harmonic for any ith degree of freedom can be chosen as a phase condition.
In order to follow the solutions of Eq. (11), a numerical continuation scheme based on a tangent predictor and an arc-length corrector is implemented. The tangent prediction uses the Jacobian of Eq. (11), at a solution point, to estimate the next solution (see e.g. Ref. [37]). The arc-length corrector equation is added to the problem, where Z 0 is an already known previous solution and s is the step length. After adding Eq. (12) to the algebraic system in Eq. (11), the Fourier coefficients a k 1 ,k 2 and b k 1 ,k 2 as well as the frequencies 1 and 2 are treated as unknowns in the minimisation problem. This continuation scheme allows the periodic and the quasi-periodic HBM approaches to follow turning points (see e.g. Ref. [37]).

Branch switching
In this subsection we develop a strategy to switch from a periodic to a quasi-periodic branch. The most straightforward way to change from a periodic to a quasi-periodic branch uses results obtained from time integration [31,38]. In this case, solutions from unstable periodic regimes are used as initial conditions for a time-marching algorithm, and the system jumps to the quasi-periodic branch naturally. This approach works well if there is only one stable quasi-periodic solution within the analysed regime, and for unstable periodic states. In more complicated cases, when e.g. several quasi-periodic states coexist, it is difficult to guarantee that initial conditions from the unstable periodic state will lead the time integration scheme to the required quasi-periodic solution.
The new approach developed in this subsection does not require any stability condition for the quasi-periodic response, and also works for systems with several coexistent stable quasi-periodic solutions. Overall, the main idea is to use information from the monodromy matrix , evaluated at an unstable periodic solution, to generate an initial guess for the nonlinear algebraic problem in Eq. (11). The steps required are illustrated in Fig. 3. Firstly, a Floquet analysis is carried out for the periodic solutions. It is well-known from the linear stability analysis that, close to the bifurcation point, perturbations h(t) increase exponentially with growth rate = 1 2 Re{ln( i )} and pulsate with frequency 2 = 1 2 Im{ln( i )}, where i is the corresponding unstable Floquet multiplier. Therefore, this value of 2 is assumed as an initial guess for the algebraic system in Eq. (11).
The second step consists in time integrating the tangent system in Eq. (8) using the real part of the corresponding unstable eigenvector of the monodromy matrix as an initial condition. Since the system is unstable with respect to this perturbation, the initial condition will increase exponentially in time with growth rate , as shown in Fig. 3. The exponential part is then removed from the response, and the final configuration is a quasi-periodic signal with frequencies 1 and 2 . It is expected that nearby the NS bifurcation point perturbations will behave similarly to this linear approximation.
The final step uses information from the periodic state. This solution, which coexists with the quasi-periodic regime, is added to the previously calculated perturbation. The final state consists of a periodic state with a small quasi-periodic perturbation on top. The spectrum of this configuration is then used as an initial guess for the algebraic system in Eq. (11). It should be noted that the presented approach has two free parameters. The first one is the perturbation amplitude. Since initial conditions for Eq. (8) uses the unstable eigenvector of , these initial conditions may assume any arbitrary amplitudes values. The second free parameter is the value of 1 itself or, in other words, the distance from the computed solution to the actual NS bifurcation point. We have observed, in the simulations of Sec. 3, that far from the bifurcation point or with large perturbations amplitudes, the nonlinear algebraic problem has experienced convergence problems. On the other hand, when the system is solved very close to the bifurcation point, or with relatively small perturbations amplitudes, the algebraic system may converge to the undesired periodic regime. The compromise between these two parameters has been chosen empirically in the present investigation.

Numerical application
We apply the methodology described in Sec. 2 to compute quasi-periodic localised solutions bifurcating from nonlinear travelling waves in cyclic and symmetric systems. We then define i as the ith linear natural frequency, and i as the corresponding mode shape. Due to the cyclic symmetry, it is well-known that most of the eigenvalues come in pairs of 2i and 2i+1 having the same values (degenerated eigenvalue), while the corresponding eigenvectors 2i and 2i+1 are just shifted in space. Using the linear mode shapes of the cyclic structure, a travelling wave force f tw (t), also known as an engine order excitation, can be defined as f tw (t) = 2i cos( F t) + 2i+1 sin( F t), (13) where F is the force frequency. The external excitation has just one harmonic, and the values of i and i+1 are directly applied as the Fourier coefficients within the HBM implementation. Moreover, the travelling wave in Eq. (13) can be seen as the superposition of two standing waves, 2i and 2i+1 , in which the time evolution is with appropriately chosen phases.
In the following we study two different minimal models for bladed disks vibrating in nonlinear travelling wave regimes. Firstly, a physical system composed of a cyclic and symmetric chain of Duffing oscillators is assumed as a physical model for a bladed-disk operating in large displacements. A similar system has been already investigated in Refs. [23,24], where it has been shown that solitons bifurcate from unstable travelling wave responses. The second example introduces a more complicated minimal model composed of two degrees of freedom per sector. In this case, travelling waves with zero and one nodal circles are possible. In both cases several coexistent quasi-periodic states may branch off in the form of localised solutions.

Bladed disk with a single degree of freedom per sector
The system under investigation is depicted in Fig. 4. It consists of N s = 24 masses m = 1 kg, connected to each other by linear springs k c = 1 N/m, and attached to the ground by linear springs k l = 1 N/m and viscous dampers c = 0.01 Ns/m. Each mass is also connected to the ground by cubic springs k nl = 0.1 N/m 3 and subjected to external forces f n . This minimal model with geometric nonlinearities can, for example, be obtained when the Von Karman theory is applied to a system of clamped beams or plates that are cyclically connected to each other by massless springs [10,39]. The N s equations of motion for the system in Fig. 4 are written as where x 1 = x N s +1 due to the cyclic symmetry.
In the following analysis the external force is assumed as an eighth engine order. For low excitation amplitudes, the system vibrates in an almost linear regime, and the frequency response functions (FRFs) show stable solutions only. For stronger forces the system looses stability due to NS bifurcations, and quasi-periodic solutions detach from the periodic branch. Fig. 5 displays a typical result calculated for a travelling wave with force amplitude of 0.025 N. The plane wave solutions have been calculated with the periodic HBM assuming a series expansion with three harmonics (h 1 = 3). In Fig. 5, solid black lines identify stable travelling waves, while the black dashed line represents unstable solutions. The plane wave FRF shows the usual stiffening effect, as expected due to the cubic springs. However, no periodic bi-stability is observed. In fact, homogeneous states loose stability due to NS bifurcations before the first turning point (around 2.00 rad/s), and recover stability again only after the second turning point (around 2.02 rad/s). Inside the unstable periodic regime, two quasi-periodic branches detach from plane waves. The first NS bifurcation, computed from the Floquet analysis, is identified with the blue line in Fig. 5. The branch is followed with the quasi-periodic HBM assuming one harmonic for the carrier wave (h 1 = 1), and five harmonics for the envelope function (h 2 = 5).
The quasi-periodic branch detaches from plane waves as weakly localised solutions but internally changes its configuration along the branch to a strongly localised dissipative soliton. This feature is highlighted by the envelope functions plotted in the insets in Fig. 5. Moreover, around 2.01 rad/s another pair of unstable Floquet multipliers is computed, leading to the second branch of localised solutions identified by the red line in Fig. 5. These solutions are calculated with the quasi-periodic HBM assuming h 1 = 1 and h 2 = 5. The branch comprises, again, localised states moving around the structure in the form of dissipative solitons.
Panels (a) and (c) of Fig. 6 display the amplitude of each harmonic, for a specific degree of freedom, and for two solutions in the two different quasi-periodic branches of 5. These solutions show that sidebands components, equivalent to h 1 = 1 and h 2 = ±1, dominate the quasi-periodic dynamics. This is a typical feature of localised states emerging via modulation instability, and it has been extensively investigated within the nonlinear Schrödinger equation framework [40][41][42]. Panels (b) and (d) show, in dashed red lines, the same solutions but in time domain. The curves show, as discussed before, localised humps preserving their shapes while they move around the system. The same graphs also identify, in black solid lines, the corresponding results computed from time-integration. The comparison between the quasi-periodic HBM and the time marching results are in very good agreement.

Bladed disk with two degrees of freedom per sector
The physical model in Fig. 7 Minimal models similar to the one in Fig. 7 are often used in the literature for the investigation of bladed disk vibrations [16,43]. Firstly, we compute the system responses due to travelling wave forces of eighth engine order and zero nodal circles. In this case, the blades and the disc are vibrating in phase, and nonlinearities are triggered by the difference between the two  amplitudes. For low excitation forces, as expected, the system behaves similarly to the underlying linear model and solutions are always stable. However, for stronger excitations, the system looses stability and experiences NS bifurcations as in the previous example of Sub. 3.1. Fig. 8 displays a typical result, obtained assuming an external force with amplitude of 0.01 N. The results have been calculated assuming h 1 = 3 for periodic solutions, while h 1 = 1 and h 2 = 5 for the quasi-periodic responses. In this case, again, no periodic bi-stability is observed, and two NS bifurcations are identified. The first quasi-periodic branch, identified by the blue line in Fig. 8, shows a solution branching off from plane waves and internally changing its configuration to a strongly localised response with one hump before merging to the travelling waves again. The second quasi-periodic branch, in red, shows similar features of Fig. 5. In this case, again, localised solutions with two humps detach from plane waves and merge back to periodic solutions.
The results obtained from the quasi-periodic HBM approach are compared to time integration in Fig. 9. Panel (c) shows an example for the first quasi-periodic branch in Fig. 8, while Panel (f) illustrates the same quantities calculated from a solution in the second branch. In this case, as discussed before, each sector has the two degrees of freedom moving in phase. The results in time domain show humps of localised solutions moving around the structure preserving their shapes. The envelope function which modulates the disk and the blades have similar shapes, although they experience different amplitudes. This feature can  be confirmed from Panels (a) and (d), which display the spectrum for a blade, and Panels (b) and (e), which illustrate the same quantities for the corresponding disk displacement. The spectra for the blade and the disk are essentially the same, differing only in amplitudes. These findings can be confirmed from the results in the four insets in Fig. 8.
Finally, we investigate the system responses due to travelling wave forces of eighth engine order and one nodal circle. Physically, the disk and the blades vibrate out-of-phase within this regime. The overall behaviour is similar to the previous examples, in which branches of localised solutions bifurcate from travelling waves when modulation instability is triggered. Fig. 10 displays a typical result obtained for an external excitation with amplitude of 0.016 N. The periodic solution has been calculated assuming h 1 = 3, while the quasi-periodic branches have been computed assuming h 1 = 1 and h 2 = 7. Within this configuration, solutions are more damped than in the previous case, in Fig. 8, since dampers c 3 dissipate more energy. Moreover, nonlinearities play the role in lower amplitude regimes due to the out-of-phase motion, and localised solutions branch off for lower  displacement regimes. In Fig. 10, three branches of quasi-periodic states detach from travelling waves. The first branch, in blue, resembles a dissipative soliton with three humps moving around the structure. However, the other two branches, in red and green, show more complicated patterns, and the modulating functions have more irregular shapes.
The previous HBM solutions are compared to time integration in Panels (c) and (f) of Fig. 11. The two displayed results are obtained from initial conditions in the first (blue) and in the third (green) branches in Fig. 10. In both cases localisation arises due to modulation of travelling waves. In Panel (c), a regular hump moves around the structure preserving its shape. In Panel (f), the envelope function which modulates the travelling waves has a more irregular form, discussed before. Although this envelope function has an irregular shape, the final configuration still consists of a quasi-periodic motion with two incommensurable frequencies. The comparison between time integration and the quasi-periodic HBM shows good agreement. Panels (a) and (d) in Fig. 11 display the spectra for the blade responses, while Panels (c) and (e) in Fig. 11 show the same quantities for the disc. The amplitude spectra for the blade and the disk are, again, very similar, although the graphs do not show any information about the differences in phase.

Conclusions and outlook
In this paper we have employed the periodic and quasi-periodic HBMs to compute nonlinear localised vibrations in cyclic and symmetric structures. The investigation has been focused on models experiencing nonlinearities induced by strong engine order excitations. The strategy has been applied to two different minimal models: (1) a chain of linearly coupled Duffing oscillators; and (2) a more complicated bladed disk model with two degrees of freedom per sector. In both cases multiple quasi-periodic coexistent states may exist, and a new strategy to track these several quasi-periodic states have been introduced. The results computed in this paper corroborate previous findings, showing that cyclic and symmetric structures subjected to large travelling wave forces may localise vibrations due to envelope dynamics.
The present paper has studied localised vibrations arising due to smooth nonlinearities. However, bladed disks of aircraft engines are often subject to friction and impact induced e.g. by root joint, underplatform dampers, shrouds and damper wires. Therefore, future investigations will also have to deal with localised vibrations involving non-smooth nonlinearities. More work also needs to be done on model accuracy. Thus, reduced-order models from fully nonlinear finite elements will also be addressed.