Dynamical analysis for the motion of a 2DOF spring pendulum on a Lissajous curve

This study examines the motion of a spring pendulum with two degrees-of-freedom (DOF) in a plane as a vibrating system, in which its pivot point is constrained to move along a Lissajous curve. In light of the system’s coordinates, the governing equations of motion (EOM) are obtained utilizing the equations of Lagrange’s. The novelty of this work is to use the approach of multiple scales (AMS), as a traditional method, to obtain novel approximate solutions (AS) of the EOM with a higher degree of approximation. These solutions have been compared with the numerical ones that have been obtained using the fourth-order Runge–Kutta algorithm (4RKA) to reveal the accuracy of the analytic solutions. According to the requirements of solvability, the emergent resonance cases are grouped and the modulation equations (ME) are established. Therefore, the solutions at the steady-state case are confirmed. The stability/instability regions are inspected using Routh–Hurwitz criteria (RHC), and examined in accordance with the steady-state solutions. The achieved outcomes, resonance responses, and stability areas are demonstrated and graphically displayed, to evaluate the positive effects of different values of the physical parameters on the behavior of the examined system. Investigating zones of stability/instability reveals that the system’s behavior is stable for a significant portion of its parameters. A better knowledge of the vibrational movements that are closely related to resonance is crucial in many engineering applications because it enables the avoidance of on-going exposure to potentially harmful occurrences.

In Ref. 1 , it is investigated how a spring pendulum (SP) moves in relation to its two controllable factors, which are the energy and the frequencies of the spring and pendulum.Additionally, the authors investigated the phenomenon, specifically the back-and-forth movement of the spring and pendulum.By treating a long spring as a physical pendulum and formulating the mass in terms of the spring constant and various spring lengths, the estimation of the mass is considered necessary in Ref. 2 to develop the resonance.For Reynolds number more than 10 4 , damped oscillations of SP model with a variable continuously diminishing mass are studied in Ref. 3 , in which the damping parameters are influenced by the mass loss rate.
5][6][7][8][9] the dynamical behavior of a few various vibrational pendulum models connected with energy harvesting devices is examined as one of the best and most effective examples of converting mechanical energy into electric energy.The study conducted in Ref. 4 investigates the vibrations of a two-degree-of-freedom spherical pendulum subjected to horizontal Lissajous excitation.By employing a mathematical model, the outcomes of numerical simulations are presented through visually appealing multi-colored maps, highlighting the behavior of the largest Lyapunov exponent.In a recent publication 5 , a groundbreaking design is presented, which encompasses a novel and sophisticated model of a double variable length cable pendulum.This model is accompanied by a meticulously designed experimental setup that incorporates elastic suspension and a counterweight mass for enhanced performance and accuracy.The investigation focuses on understanding the intricate dynamics that arise from the influence of varying lengths on the frequency and amplitude of vibrations.The study conducted in Ref. 6 explores the interplay between parametric excitation and self-excited vibration within discontinuous systems.Through the use of a separate electromagnetic harvesting device, the pendulum's structure is altered in Ref. 7 , where harvesting is dependent on the magnet in the coil oscillating.It is noted that the harvester's effectiveness at both energy gathering and vibration reduction has increased.In Refs. 8,9, 3DOF harvesting models have been examined.The model in Ref. 8 is composed of two linked components: a nonlinear damping SP combined with an energy harvesting device and a nonlinear Duffing oscillator, while the other one in Ref. 9 is formed up of two connected parts: the first is coupled to a piezoelectric transducer, which transforms stresses and oscillations into electrical power, whereas the second is a nonlinear damping SP.
A semi-analytical approach was used in Ref. 10 to study the periodic movements of a periodically driven nonlinear SP, and the relevant stability and bifurcation analysis of these movements.After providing a consistent magnetic field in one direction, the motion of a SP is assessed in Ref. 11 .The AS, resembling Foucault's pendulum, are also obtained when the heavy pendulum ball and delicate spring are taken into account.
The bifurcation phenomena at its state of equilibrium are discovered after examining how the magnetic field affects the stability of the SP.The implicit mapping approach is used in Ref. 12 to calculate semi-analytically the entire bifurcation dynamics of period-3 motions to chaos.The harmonic amplitudes that change with excitation amplitudes to examine the complexity of periodic motion have been obtained.For the purpose of reducing vertical disturbances, a rotating pendulum absorber is suggested in Ref. 13 .By changing the rotating speed, the pendulum absorber's characteristic frequency may be dynamically modified across a large range.5][16][17] .The AMS 18 is used to obtain the essential AS of the EOM, where the resonance situations are categorized and examined.
To evaluate the impact of the approximations of higher-order on the chaotic processes of a multi-DOF dynamical system with weak nonlinearity, Ref. 19 examines a harmonically stimulated SP in a circular trajectory with internal resonance.A parametric and externally excited 2-DOF weakly nonlinearly system is investigated in Ref. 20 .There is a noticeable energy transfer between modes of vibration, where the selected resonance instance and the resonance conditions have been analysed and determined.The authors of Refs. 21,22build on the behavior of kinematically nonlinear excited SP, where its pivot point travels along an elliptic route.Simultaneous resonances have been studied in view of the exposed resonances circumstances.
In Ref. 23 , two different approaches have been used to resolve the motion of the nonlinear SP problem.According to the asymptotic analysis, three-time scales are utilized to get the AS with a respectable small error.In Ref. 24 , a generic model of a nonlinear damped excited SP is examined, in which its pivot point has been constrained to move along an ellipse trajectory with a stationary angular velocity.The AS are obtained up to the third correction.In Ref. 25 , the quadratic polynomial approximation was used to create an approximate controlling system with trigonometric nonlinearities.The presented unique approach ensures that the trigonometric functions are approximated with adequate precision not only around a specified point but also throughout the entire predetermined interval, contrary to the approximation accomplished using Taylor series.Thus, the suggested approximation is considered an approach that ensures better predictions for resonance responses in nonlinear mechanical systems.An approximation differential system is used to analyze the pendulum damper, which is modeled as a severely nonlinear auto-parametric system with 2DOF.As a foundation for the in-depth analytical analysis, the nature of the numerical solutions (NS) in the resonance state is examined.The resonant solution's stationary and non-stationary properties are described in Ref. 26 .In Ref. 27 , the asymmetrical damping of a pendulum and its nonlinear properties, have been represented mathematically.Three distinct forms of the resonance domain were studied, and it was discovered that the excitation amplitude and the pendulum's dynamic characteristics had a substantial impact on each type's attributes.
In Ref. 28 , a few unusual states that can occur when a ball is moving in a sphere-shaped cavity acting as a passively tuned mass damper for thin engineered structures have been illustrated.Three non-holonomic restrictions are placed under horizontal additive kinematic excitation in a 6DOF system.The controlling differential system is determined using the Appell-Gibbs method 29 .In Ref. 30 , two viscous dampers and two linked nonlinear springs in series are used to analyze the forced planar vibrations of an attached particle at a supported point.The third-order correction law is proposed as the constitutive connection for the elastic forces of each spring.Three terms of Taylor series are used to simulate the resulting geometric nonlinearity from the pendulum's transverse motion.In Refs. 31,32, the frequency responses of a 2DOF nonlinear dynamical model that simulates the motion of a damped SP in an inviscid fluid flow are examined and discussed.This work's main objective is to investigate the motion of a 2DOF non-linearly damped SP system.It is assumed that two harmonically generated forces act in both the transverse and longitudinal directions, as well as a harmonic external moment at the pivot that restricts the pendulum motion to being on a Lissajous curve.The regulating EOM are derived applying Lagrange's equations from the second type.For a higher level of accuracy, the EOM are analytically solved using AMS.The accuracy of the analytical solutions is determined by comparing them to the numerical ones that were derived using the 4RKA.In regard to the removal of secular factors, the solvability criteria and the ME are found.In order to confirm that the fixed points at steady-state solutions are stable or not, the RHC are applied.The non-linear stability analysis is used to expose various regions of stability or instability.A graphical examination of numerous plots associated with separate time-history curves, resonances, and stability zones is used to show how the system behaves.In many engineering applications, a deeper understanding of the vibrational motions that are closely associated with resonance is essential because it reduces the possibility of being continually exposed to potentially damaging events.

Formulation of the dynamical system
The major objective of the present section is to derive the governing EOM of the examined dynamical system.This system consists of a non-linear damped SP with a massless normal length l 0 and stiffness k and k 1 .The suspension point O 1 of this pendulum is constrained to move on a route of a Lissajous curve in an anticlockwise direction, in which it moves independently and harmonically in two orthogonal directions, while the other point is connected with pendulum mass m and oscillates in a plane.It is taken into account that point O , two orthogonal axes, OX and OY , that are pointed, respectively, vertically and horizontally, are considered, see Fig. 1.
The coordinates (x, y) describing the kinematic motion of the point O 1 are x = R x cos x t and y = R y sin y t , where R x , R y , x , and y represent known parameters.The planar motion of the investigated system is considered under the action of two perpendicular harmonic forces F 1 (t) = F 1 cos � 1 t and F 2 (t) = F 2 cos � 2 t in the radial and transverse directions of the spring, respectively, in addition to an anticlockwise moment M(t) = M 0 cos � 0 t at O 1 .Here, F 1 , F 2 , M 0 and 1 , 2 , 0 are amplitudes and frequencies of these forces and moment.Let J(t) rep- resents the elongation spring after time t , and C j (j = 1, 2) are the viscous damping coefficients for longitudinal vibration and the swing one.These coefficients prevent the system from reaching a critical case in both vibration directions.
The following expression provides a foundation for formulating the system's Lagrangian where g denotes the acceleration of gravity, denotes the swing's angle, l = l 0 + J r , J r = mg k denotes the spring's static elongation, and dots are the differentiation regards t .Equation (1) can subsequently be used to derive the governing EOM using the second type of Lagrange's equations below where q stands for the system's generalized coordinates and Q q represents a non-conservative generalized force that may be expressed as follows (1) (2) d dt where the dot denote the derivatives regarding τ .The initial circumstances that constitute with the above EOM may be stated as follows pertaining

The proposed method
In the current portion, the ASM can be used to achieve the solutions of the EOM ( 6) and ( 7) analytically, categorize the resonance situations, and get at the ME as well as the requirements of solvability.Therefore, we start by looking at the system's vibrations in close proximity to its static equilibrium point 33 .Consequently, it is possible to approximate the trigonometric functions according to the following expressions Then, using a tiny parameter called 0 < ε << 1 , we may represent the damping coefficients c j (j = 1, 2) , forces' amplitudes f j , moment's amplitude m 0 , and the parameters r x , r y , α as follows In a similar vein, we suppose that the vibrations' amplitudes ℜ and are of order ε .Then, one can write The expressions for the functions R and ˜ according to the AMS 18

can be represented as follows
It is important to note that τ n = ε n τ (n = 0, 1, 2) expresses new time scales that are reliant on τ , where τ 0 is rapid time scale whereas τ 1 and τ 2 are the slow ones.Additionally, due to their tiny size, terms of O(ε 2 ) and higher orders have not been taken into account.To deal with the EOM ( 6) and ( 7) and the supposed solutions (11) and (12), we need to transform the time derivatives in (5) to be in relation to the time scales τ n , as follows Vol.:(0123456789) www.nature.com/scientificreports/Substituting (8) through ( 13) into ( 6) and ( 7) and equating coefficients of various powers of ε with zero to obtain the below sets of partial differential equations (DEs).

Scientific
Regarding order (ε) Regarding order (ε 3 ) The prior sets include six successively solvable second-order non-linear partial DEs, which emphasize the significance of the solutions of the first set (14).As a result, these equations' generic solutions take the below forms where i = √ −1 , D j (j = 1, 2) are unknown dependent complex functions on τ j and D j are their complex conjugate.
Substituting the ( 17) and ( 18) into the second set of partial DEs (15) and removing the produced secular terms to obtain Consequently, the second-order approximation can be expressed as follows.
where the symbol c.c. refers to the aforementioned terms' complex conjugate.
To calculate the next requirements of solvability, substitute ( 17)-( 21) into the third set of partial DEs ( 16) and then delete the apparent secular terms.( 14) www.nature.com/scientificreports/Consequently, the third-order approximation may be written as follows The circumstances of eliminating secular terms ( 19), (22), and ( 23) can be used to calculate the functions D j (j = 1, 2) .One may readily find the asymptotic AS ℜ and up to the third approximation according to the substitution of solutions [( 11), (12)], [( 17), (18)], [( 20), ( 21)], and [( 24), (25)] into (10).Now, it's important to highlight that the obtained AS remain acceptable when their dominators depart from zero 34 .However, resonant scenarios emerge when these dominators get closer to zero.As a result, one can categorise these scenarios as follows.
It should be emphasized that if any of the preceding resonance scenarios occur, the behavior of the researched system would be difficult.Therefore, it would be necessary to alter the employed approach.
To address this issue, we will look into two fundamental external resonances that occur at the same time.
These relationships demonstrate how closely p 1 to 1 and p 2 to ω .To achieve this purpose, it is important to provide dimensionless values σ j (j = 1, 2) that are known by detuning parameters (which specify the separation between the strict resonance and vibrations) as follows In light of this, we can express σ j according to the use of ε as To obtain the following solvability requirements for the second and third-orders equations, substitute ( 27) and ( 28) into ( 15) and ( 16), and then eliminate terms that yield secular ones.
A closer look at the aforementioned solvability requirements reveals that they combine to generate a system of four nonlinear partial DEs involving functions D j (j = 1, 2) .In addition, these functions are exclusively dependent on the slow time scale τ 2 , as explored in the first two requirements in (29).Hence, we provide D j in the polar form shown below (28) www.nature.com/scientificreports/where ψ j (j = 1, 2) and ãj are real functions that describe the phases and amplitudes of the solutions R and ˜ .The modelling procedures mentioned above show that the first-order derivatives of D j can be expressed as follows In the context of (31), one can transform the partial DEs in (29) into ordinary DEs.Introducing (10), ( 13), (30), and (31), as well as the next adjusted phases into (29), distinguishing between real and imaginary portions to obtain the below autonomous system of four first-order ordinary DEs with regard to a j and θ j It is obvious that the aforementioned system of Eqs. ( 33)-( 36) characterizes the ME for the two resonances that are being analysed concurrently.This system has been solved numerically to obtain the solutions a j (τ ) and θ j (τ ) , which are graphed in portions of Figs. 2, 3, 4 and 5 according to the following data of the used parameters (30)  Curves of these figures are calculated when the damping parameter c 1 (= 0.05, 0.07, 0.09), c 2 (= 0.07, 0.09, 0.12), and the frequencies ω 1 (= 5.07, 5.18, 5.29), ω 2 (= 2.98, 3.21, 3.5) have various values, as seen in potions (a), (b) and (c), (d) of these figures, respectively.These curves behave in a decaying manner, and they reach the stage of full stability at the end of the first fifth of the studied time interval when the aforementioned values are taken into account.It is noted that a 1 has been impacted by the change of the values of c 1 and ω 2 , as drawn, respectively, in Figs.2a,d, 4a,d.In the same context, the temporal history of the amplitude a 2 is influenced by the change of the values c 2 and ω 2 , as seen, respectively, in Figs.3b,d, 5b,d.A closer look at the other parts of Figs. 2, 3, 4 and 5, one can observe that they haven't any variation with the change of c j and ω j .The reason for the change or non-change is due to the mathematical composition of the equations of system (33)- (36), as the first and third equations are dependent on c 1 and ω 2 , respectively.Whereas they do not explicitly  depend on the variable c 2 and ω 1 .Similarly, the second and fourth equations of the same system is independent on c 1 and ω 1 , in which they are depend on c 2 and ω 2 .
The projections of the plotted curves in Figs. 2, 3, 4 and 5 in the planes θ 1 a 1 and θ 2 a 2 are drawn in portions of Figs. 6 and 7.The behaviors of these curves have the form of spiral curves that are directed towards one point, which means that the functions described by these curves are stable.This conclusion is consistent with the  Figures 8 and 9 present, respectively, the achieved analytical solutions ℜ(τ ) and �(τ ) to highlight the temporal behavior of these solutions while taking into account the prior values of the system's parameters.This behavior has the form of quasi-periodic waves.It must be mentioned that these have been impacted more by the various  values of the frequencies ω 1 and ω 2 than the damping parameters c 1 and c 2 .The accuracy of the analytical solu- tions is evaluated by comparing them to the numerical ones of the original EOM that were produced using 4RKA according to the plotted curves in Fig. 10 at c 1 = 0.05, c 2 = 0.07, ω 1 = 5.07, and ω 2 = 3.5 .The comparison demonstrates excellent agreement between both solutions.

Solutions at the scenario of steady-state
The purpose of the current section is to investigate the steady-state oscillations of the considered dynamical system.In essence, temporary processes will cease to exist under the impact of damping, and the steady-state oscillations will then be generated [35][36][37][38] .Therefore, we regard the left-hand side of the system of ME ( 33)- (36) as zero.As a consequence, the equations of this system have been transformed into the algebraic equations shown below (37)    www.nature.com/scientificreports/After removing the adjusted phases θ j (j = 1, 2) from the system of Eqs. ( 37)-( 40), the next nonlinear alge- braic equations regarding the parameters of detuning σ j and adjusted amplitudes a j are obtained.
One of the most important aspects of steady-state oscillations is to analyse the stability.To analyse such a scenario, the system's behavior will be evaluated in a relatively close neighbourhood area to the locations of fixed points.To accomplish this purpose, we consider the next substitutions in ( 33)- (36)   where a j0 (j = 1, 2) and θ j0 denote the solutions at the steady-state, while a j1 and θ j1 denote relatively minor disturbances in comparison to a j0 and θ j0 .Consequently, after linearization, one gets Remembering that a j1 and θ j1 are defined, respectively, above as unknown perturbed functions of ampli- tudes and phases in the preceding system.Then we are able to outline their solutions as a linear superposition of q s e τ (s = 1, 2, 3, 4) , where q s represent constants and expresses the eigenvalue of these functions.If the solutions a j0 and θ j0 are stable asymptotically, then the real components of the roots of the yielded characteristic equation of the system (43)-( 46) must be negative where (41) (47 To determine the requisite stability criteria for the solutions in a certain steady state, the following RHC 18 can be used

The stability analysis
This section explores the stability of the examined system using the non-linear stability approach of Routh-Hurwitz.It must be remembered that the system under consideration consists of a moving, nonlinear, damped spring pendulum in a Lissajous route, which is influenced by an external harmonic moment M(t) as well as two perpendicular forces F 1 (t) and F 2 (t) .The requirements of stability are applied alongside the simulation of the system's non-linear evolution.A number of variables, such as damping coefficients c j (j = 1, 2) , frequencies ω j , and detuning parameters σ j , have been discovered to have a vital influence in the stability criteria.
The stability diagrams of the system of Eqs. ( 33)-( 36) are obtained by performing certain actions with different parameters of the system.The variation of adjusted amplitudes a j with time is plotted for different para- metrical zones, and its characteristics are presented using phase plane paths.
Figures 11 and 12 have been drawn, respectively, in planes σ 1 a 1 and σ 1 a 2 to represent the frequency response curves (FRC) when c j and ω j have different values in addition to the value of the detuning parameter σ 2 which is computed according to the relation σ 2 = p 2 − ω .In more details, curves of Figs.11a and 12a summarise the effect of the different values of c 1 (= 0.05, 0.07, 0.09), on the generated curves.Examining of these figures reveal that each curve contains only one critical fixed point over the whole domain as tabulated in Table 1.The stability and instability zones are discovered in the ranges σ 1 < 0.065 , and 0.065 ≤ σ 1 , respectively.It is critical to note that the solid curves represent the range of stable fixed points, whereas the dashed lines depict the range of unstable fixed points.
According to the positive impact of the various values of the other damping parameter c 2 (= 0.07, 0.09, 0.12), Figs.11b and 12b are drawn to display the FRC at these values.As aforementioned, one critical fixed point is observed for each curve, in which stable and instable fixed points at c 2 = 0.07 are generated, respectively, in the ranges σ 1 < 0.065 and 0.065 ≤ σ 1 .Whereas, at c 2 (= 0.09, 0.12) one finds other regions of stability and instability at the ranges σ 1 ≤ 0.066 and 0.066 < σ 1 .The drawn FRC in Figs.11c and 12c show the good influence of vari- ous values of the frequency ω 1 (= 5.07, 5.18, 5.29) on the behavior of the stability and instability areas, in which there exists a single fixed point for each curve.It is observed that the areas of stability are found in the ranges σ 1 ≤ 0.065, σ 1 ≤ 0.056, and σ 1 ≤ 0.048 , while the instability areas of the fixed points are generated in the range 0.065 < σ 1 , 0.056 ≤ σ 1 , and 0.048 ≤ σ 1 .Other stability and instability regions have been plotted at different values of ω 2 (= 2.98, 3.21, 3.5) as seen in Figs.11d and 12d.The stable fixed points are found in the ranges σ 1 ≤ 0.046, σ 1 ≤ 0.058, and σ 1 ≤ 0.065, while the unstable ones occurs in the ranges 0.046 < σ 1 , 0.058 < σ 1 , and 0.065 < σ 1 .
According to the value of the detuning parameter σ 1 , which is calculated using the relation σ 1 = p 1 − 1 , Figs. 13 and 14 have been drawn, respectively, to depict the FRC in planes σ 2 a 1 and σ 2 a 2 when c j (j = 1, 2) and ω j have various values.The range of stable fixed points is shown by the solid lines, while the range of unstable (49) ones is shown by the dashed lines.These figures illustrate that each curve contains critical and peak fixed points, which are tabulated in Tables 2 and 3.

Non-linear analysis
The purpose of this section is to clarify the properties of the non-linear amplitudes of the system of Eqs. ( 33)- (36)  and look into its stability.Consequently, we take into account the below transformation 39,40 where U j and V j are, respectively, the amplitudes' real and imaginary components.
Separating the distinct parts that yield from the substitution of (50) into ( 33)-( 36) to obtain where (50) Table 2. Critical and peak fixed points for the curves of Fig. 13 when σ 1 = p 1 − 1.  www.nature.com/scientificreports/ The modified amplitudes have been justified over an entire period of time in distinct parametric zones based on the previously mentioned data of the used parameters, in which their properties may be displayed in the curves of phase plane, as shown in Figs. 13, 14, 15, 16, 17 and 18.  Figures 15, 16, 17 and 18 show how the new modified amplitudes u j and v j change over time τ according to the numerical solutions for the system of Eqs. ( 51)-(54) when c j and ω j have different values.Decay waves have been graphed in light of these values until they become nearly motionless at the end of the time period.It is noted that these curves behave in a stable manner, which can be asserted when the projections of these curves are plotted in a suitable phase plane.Therefore, curves in Figs.19 and 20 have been drawn to explore how the projections of u j and v j are plotted in the planes u j v j when the aforementioned values of c j and ω j are considered.The behavior of the resulting curves shows spiral patterns oriented to one point for each curve, indicating that this behavior is steady and free of chaos.

Conclusion
This work has focused on analysing the planar movement of a spring pendulum with two degrees of freedom that undergoes vibrations, in which its pivot point is confined to move along a trajectory that resembles a Lissajous curve.By utilizing the system's coordinates, the EOM for the system have been successfully derived using Lagrange's equations.The AMS technique has been utilized to obtain highly refined solutions for this system,  surpassing previous approximations.These solutions have been contrasted with the obtained NS through the 4RKA method to reveal the exceptional precision achieved with the employed perturbation approach.The classification of resonance situations and the development of ME have been accomplished, taking into account the solvability constraints.Therefore, the solutions for steady-state scenarios have been verified.The RHC has been utilized to evaluate and plot both stable and unstable regions.The obtained outcomes, including FRC and stability zones, are displayed and visually depicted to evaluate the beneficial impact of various physical parameter values on the behavior of the analysed system.Upon scrutinizing the stability and instability zones, it becomes evident that the behavior of the system remains stable for a significant portion of its parameters.Furthermore, the nonlinear stability analysis of the adjusted amplitudes has been examined to reveal their stationary behavior.

Figure 1 .
Figure 1.Dynamical sketch of the system.

Table 1 .
Critical and peak fixed points for the curves of Figs.11 and 12 when σ 2 = p 2 − ω.