Novel Asymptotic Solutions for the Planar Dynamical Motion of a Double-Rigid-Body Pendulum System Near Resonance

The planar dynamical motion of a double-rigid-body pendulum with two degrees-of-freedom close to resonance, in which its pivot point moves in a Lissajous curve has been addressed. In light of the generalized coordinates, equations of Lagrange have been used to construct the controlling equations of motion. New innovative analytic approximate solutions of the governing equations have been accomplished up to higher order of approximation utilizing the multiple scales method. Resonance cases have been classified and the solvability conditions of the steady-state solutions have been obtained. The fourth-order Runge–Kutta method has been utilized to gain the numerical solutions for the equations of the governing system. The history timeline of the acquired solutions as well as the resonance curves have been graphically displayed to demonstrate the positive impact of the various parameters on the motion. The comparison between the analytical and numerical solutions revealed great consistency, which confirms and reinforces the accuracy of the achieved analytic solutions. The non-linear stability analysis of these solutions have been examined and discussed, in which the stability and instability areas have been portrayed. All resonance cases and a combination of them have been examined. The archived results are considered as generalization of some previous works that are related to one rigid body and for fixed pendulum’s pivot point.


Introduction
The non-linear dynamical systems are studied in many research works due to its importance in many fields, for example in industrial applications, biology, and medicines, see [1][2][3][4]. Perturbation approaches [5] can be used to obtain the considered work [10] through the effects of an external force on the dynamical motion of the spring pendulum.
Furthermore, other types of pendulums have been studied in several works e.g. [12][13][14][15]. In [12], the authors used a dynamic absorber that can be moved along the longitudinal of transverse direction of the excited pendulum. High efficiency of this absorber is found for very little damped systems. The influence of the fluid flow on the motion of simple pendulum and a spring one is examined in [13] and [14] respectively. On the other hand, the motion of a vibrating absorber pendulum with a rotating base is investigated in [14]. The stability analysis is performed on a linear system where the rotatory motion is adjusted according to the law of proportional control. The chaotic motion of an undamped spring pendulum is examined in [16] when its connected point moves in a circular trajectory. The extension of this work is found in [17] when the authors relied on the fact that the motion of the suspension point is on an ellipse. Some special examples for various motion of this point are presented. Authors of [18] have turned their attention to investigate the movement of the suspension point of a damped spring with linear stiffness on a Lissajous curve under the action of harmonically force and moments. In [19], the motion of a tuned absorber connected with simple pendulum is examined close to resonance. The MSM is applied to obtain the approximate solutions of [14,[16][17][18][19] and represented in various figures to show the significance of the applied forces and moments on the considered motion of each work.
The motion of a single damped rigid-body-pendulum (RBP) is examined in different research works such as [20][21][22][23][24][25], in which its suspension point is considered to be fixed in [20], while it was considered in motion on different paths as in [21][22][23][24]. The stiffness coefficients of the spring are constraint to be linear and non-linear in [20,21,23] and [22,24] respectively. The investigated motions are influenced by external forces and moments and studied analytically using the MSM. Resonance cases are categorized and examined to obtain the solvability conditions and the modulation equations (ME). The nonlinear stability of these equations is introduced in [24] to reveal the stability and instability ranges of the fixed points at the steady-state case. The numerical solutions of the RBP is investigated in [25] using the Runge-Kutta algorithms from fourth-order.
There are many theoretical and experimental studies discuss the motion of double pendulums. In [26], a modified midpoint integrator is used to analysis the planar movement of a double pendulum numerically, in which Poincaré sections and bifurcation diagrams are presented for specified values of energy. The internal resonance of this pendulum has been studied in [27]. The bifurcation structure of a parametric excited double pendulum at small amplitudes of oscillation is studied in [28] while the resonance and non-resonance cases of this pendulum are analyzed in [29]. The rotatory motion of a model consists of two pivotal pendulums about horizontal axes as investigated in [30]. It is supposed that the model has four relative equilibrium locations in which the stabilities of these locations are studied. The influence of the gravity force on the double pendulum's motion with a vibrating suspension point is investigated in [31].
This work concerns a novel problem of the planar motion of dynamical system consists of 2DOF double RBP, whose pivot point has been constrained to move in a Lissajous curve. The non-linear differential EOM are derived utilizing the second kind of Lagrange's equations. Higher estimations of the analytical solutions are achieved using the MSM. These solutions are compared with the numerical ones, which are obtained using the fourth-order Runge-Kutta method, to show the great accuracy between them. In light of the categorized resonance cases; the solvability conditions and the ME are obtained. Criterion of Routh-Hurwitz is applied to check the stability of the fixed points of the steady-state solutions. The non-linear stability analysis of the ME is examined and discussed through the stability and instability ranges of the plotted curves of the frequency responses. The time histories of the solutions are graphed as well as the curves of resonance to explore the significance of various values of the physical parameters on the systems' behavior. One of the most important applications of this work in our daily life is to utilize the vibrational motion for the robotic devices, pump compressors, rotor dynamics, and transportation devices.

Dynamical Modeling
The investigated dynamical model in this work consists of two pivotally connected double RBP of mass m j (j = 1, 2) from the points O j. The first pivot at O 1 is restricted to move in a Lissajous curve in the direction of anticlockwise, while the second one at O 2 is considered a connection point between the two bodies, see Fig. 1. Therefore, the kinematic coordinates of O 1 are where R x , R y , Ω x , and Ω y are known parameters.
It is taken into consideration that the geometric quantities e j and l 1 express the distance of the gravities centers z j (that are calculated from the rotation's centers) and the pendulum's link length, respectively. Moreover, let J j represent the moment of inertia, j are the angular position variables that represent the chosen generalized coordinates of the model, and g is the acceleration of gravity. (1)

3
To study the motion, we consider that the motion of the model is planar in addition to action of two k i n e m a t i c m o m e n t s M Here Ω 1 and Ω 2 denote the forcing frequencies of M 1 (t) and M 2 (t) respectively. Referring to above, we can introduce the potential and kinetic energies V and T as follows where dot represents the derivative with respect to time t.
According to the Lagrangian L = T − V , the governing EOM were obtained from the following second type of Lagrange's equations where ̇j are the generalized velocities corresponding to the generalized coordinates j and Q j are the generalized forces which can be written in the forms Introducing the dimensionless forms of the generalized coordinate, parameters, frequencies and time as follows Making use of Eqs. (2), (3), (4), and (5) we get the following dimensionless forms of the EOM It must me noted that the derivatives in the previous system of Eqs. (6) and (7) are considered with respect to the dimensionless time parameter . Moreover, This system has two second-order non-linear differential equations regarding 1 and 2 .

The Approved Method
This section presents the MSM that has been used to estimate the analytic solutions of the EOM and to categorize the resonances cases. Consequently, we can approximate the trigonometric functions of 1 and 2 up to third-order using Taylor's series, which are valid at the position of static equilibrium. Then Eqs. (6) and (7) will take the forms Now, we introduce the coordinates 1 and 2 as functions of small parameter 0 < << 1 as follows as where and are the new variables that can be sought as follows [32] where n = n t (n = 0, 1 , 2) are different time scales.
Based on these scales, the time derivatives have been transformed regarding to these scales as follows Furthermore, we propose that where B j ,Ñ j ,f j ,s 1 ,r 1x ,r 1y , and ̃1 are parameters of order unity.
Inserting Eqs. (10)-(13) into Eqs. (8) and (9), equaling the coefficients of the same power of to get the following three groups of partial differential equations (PDE) according to the orders of .
The first group of PDE of order ( ) The second group of PDE of order ( 2 ) The third group of PDE of order ( 3 ) The previous Eqs. (14)-(19) constitute a system of six non-linear PDE. The solutions of these equations can be gained consecutively. Therefore, the solutions of Eqs. (14) and (15) are where A j (j = 1, 2) are complex functions of the scales j that can be determined latter, while A j indicate to their complex conjugate.  (20) Introducing Eqs. (20) and (21) into Eqs. (16) and (17), then omitting terms that yield the secular ones to obtain the conditions for omitting these terms as follows Based on these conditions, the second-order solutions become where the letters CC are the complex conjugates of the preceding terms.
The conditions of the elimination of the secular elements from the third-order approximate Eqs. (18) and (19) have the forms Thus, the third-order approximations 3 and 3 can be written as follows (22) 2i It is worthy to mention that the functions A j (j = 1, 2) can be calculated from the eliminating conditions of secular terms Eqs. (22), (25), and (26).

Modulation Equations (ME)
This section's purpose is to characterise resonance cases, to gain the conditions of solvability, and to obtain the ME. It is known that if the denominators of the solutions Eqs. (27) and (28) approach to zero [33], then these solutions break down and the resonance cases are noticeable. As a result, it falls under the following categories: (i) The major external resonance appears when p 1 ≈ 1 , p 2 ≈ , (ii) T h e i n t e r n a l r e s o n a n c e a r i s e s w h e n ≈ 1 , p x ≈ 2 , p x ≈ 2 , p y ≈ 1 , p y ≈ .
If any of the preceding conditions are met, the system's dynamical performance will be extremely tough. It's worth noting that when oscillations escape from resonances positions, the AS in the previous section are still valid.
First, we'll look at two major external resonances that are presenting at the same time, in which the nearness of p 1 , p 2 to 1 , can be expressed by using the detuning parameters of 1 and 2 (which refer to the distance between vibrations and the stringent resonance [34]) as follows Therefore, we can introduce them in terms of O( ) as follows (29) Substituting Eqs. (29)-(30) into Eqs. (16)- (19), and eliminating terms that lead to secular ones to get the conditions of solvability of the second-order of approximation as in Eq. (22) and the third-order of approximation as follows Making use of conditions Eq. (22) into Eq. (31), to yield It is noted from Eqs. (22) and (32) that the conditions of solvability of the investigated model consisting of four PDE regarding the functions A j (j = 1, 2) , in which conditions (30) j =̃j; (j = 1, 2).
Eqs. (22) and (32) reveal that A j depend upon the slow time parameters 1 and 2 . Then we can represent A j as in the following polar form Let us consider the following modified phases [35] (33)   Substituting Eqs. (13), (33), and (34) into Eqs. (22) and (32), and separating real and imaginary parts to obtain the following four first-order ODE constituting the ME This system has been inspected to show that it has the solutions a j (j = 1, 2) and j which describe the modulations of the amplitudes and the phases, respectively. These solutions are graphed for various selected values of the system's parameters as displayed in Figs. 2, 3, 4 and 5 in accordance with the following data      It is notable that, we have closed symmetric trajectories to some extent which assert that the gained AS have a stable behaviors, which are predicted before, during the tested interval of time.

Steady-State Solutions for the Case of External Resonances
The purpose of the present section is to inspect the steadystate oscillations of the scrutinized dynamical system when the oscillations of the transitory processes have disappeared. Therefore, the conditions of the steady-state can be obtained . As a result, regarding the unknowns, a j and j we can obtain a system of four algebraic equations as follows: Omission of the modified phases j from this system produces the next two non-linear algebraic relationships between amplitudes a j and the frequency response functions described by the detuning parameters j in the following forms It should be noted that, the steady-state part is considered as a significant part in the procedure of stability examination. Therefore, we concentrate our investigation on a region (37) that close to the fixed points. To achieve this goal, let us consider the following forms of amplitudes and phases [37] where a 10 , a 20 , 10 , and 20 are the solutions of steady state, while a 11 , a 21 , 11 , and 21 are the corresponding perturbations. Inserting Eq. (38) into Eq. (35), after linearization we can obtain (38) a 1 = a 10 + a 11 , 1 = 10 + 11 , a 2 = a 20 + a 21 , 2 = 20 + 21 , To achieve the solutions of the given system Eq. (39), we can write the unknown functions a j1 (j = 1, 2) and j1 in exponential forms as q s e where q s (s = 1, 2, 3, 4) are constants and is the appropriate eigenvalues of these functions. In this case, the steady-state solutions a j0 and j0 (j = 1, 2) are stable asymptotically if the real parts of the roots of the following characteristic equations must be negative [38] where (40) 20 {−8f 1 f 2 ( 2 − 1) sin 10 sin 20 + f 1 a 20 cos 10 [3a 2 10 ( 2 − 1)   It is noted that Γ s (s = 1, 2, 3, 4) depend on some parameters like a j0 , j0 , and f j . The stability conditions at steadystate solutions can be written as follows using the criterion of Routh-Hurwitz [39] (42)

The Stability Analysis for the Case of External Resonances
In the current section, the approach of the non-linear stability is utilized to study the motion of the investigated dynamical system of the double RBP that subjected to the harmonic external moments M 1 (t) and M 2 (t) . In addition to simulating the equations of the non-linear system, the stability criteria are also performed. It is observed that the parameters of detuning 1 , 2 , and the frequency have a positive impact on the destabilization of the stability criteria. To sketch the stability areas or ranges of system Eq. (37), a specified procedure with various settings of the system was tested. The changes of the amplitudes a 1 , a 2 via time are plotted for distinct parametric regions. Therefore, Figs

Non-linear Analysis for the External Resonance Case
This section aims to study the characteristic of the non-linear amplitude for the system of Eq. (35) and explores their stabilities. Therefore, the accompanying transformation is taken into account [40] (43)

A Combination Case of Internal and External Resonance
We are going to study the case of combination of external and internal resonances. The closeness of , p 2 to 1 , is expressed through introducing the detuning parameters of 1 and 2 as follows Therefore, we can introduce them in terms of O( ) as follows Substituting Eqs. (45)-(46) for Eqs. (16)- (19) and removing terms that lead to secular ones, the following solvability criteria are obtained as in Eq. (22) and as in the following conditions (46) j =̃j (j = 1, 2).  (47) show how functions A j are dependent on 1 and 2 that can be stated in polar form as Let us consider the following modified phases [35] (47)

Conditions (22) and
Inserting Eqs. (13), (48) and (49) into Eqs. (22) and (47), and separating the real and imaginary parts, we obtain the ME of four first-order ODE as follow The solutions b j (j = 1, 2) and j of this system describes represent the modulations of the amplitudes and the phases, respectively. These solutions have been graphed for various selected values of the system's parameters as shown in Figs. 27 and 28 using the following data

3
A closer look at Figs. 27 and 28 reports that these figures are plotted at 1 = 0.11, 2 = 0.001, and f 1 = 0.002 when f 2 (= 0.001, 0.0015, 0.002) to reveal the time histories of the amplitudes b j (j = 1, 2) and the modified phases j . It is clear that when the value f 2 increases, periodic waves are produced, in which the amplitudes of the waves describing b 1 and b 2 increase and the number of oscillations remain stationary to some extent as seen in Fig. 27, while the oscillation's number of 1 increases with the increasing of time as drawn in Fig. 28b. Moreover, there is no significant change of the waves characterizing 2 with the change of f 2 values, and the behavior of the waves is decreasing with the time as observed from Fig. 28b.
The achieved AS of the rotation angles 1 and 2 through time are drawn in Fig. 29 (for the combination case of internal and external resonance) to examine the variation of these  Fig. 29b. These waves express that the motion has a steady behavior and the AS has a manner of stationary behavior. The phase plane figures that describe the stability of the system's motion have been portrayed in Fig. 30, which reveal the relation between the AS and their first-order derivatives at f 2 (= 0.001, 0.0015, 0.002). Closed trajectories are observed to indicate that the gained AS have a stable behaviors.

Verification of the Steady State Solutions
To examine the solutions at the case of steady state, we consider as studied previously in Sect. 5 the zero values of the first derivatives of amplitudes and modified phase i.e.,    Cancellation of the modified phases j from preceding system Eq. (51), two non-linear algebraic relationships in terms of b j , j and the frequency response functions are obtained in the forms Let us consider the following forms of amplitudes and phases [37] where b 10 , b 20 , 10 , and 20 are the solutions of steady state, while b 11 , b 21 , 11 , and 21 are the corresponding perturbations. Substituting Eq. (53) into Eq. (50), and making linearization we can obtain (53) The solutions of the preceding system can be obtained if we express the unknown functions b j1 and j1 (j = 1, 2) in an exponential form as C s e ℏ . Here C s (s = 1, 2, 3, 4) are constants and ℏ is the functions' eigenvalues. In such case, the steady-state solutions b j0 and j0 (j = 1, 2) will be asymptotically stable, if the roots' real parts of the following characteristic equations have negative values [38] where Taking into account the stability conditions at steady-state solutions by using the criteria of Routh-Hurwitz [39] are

The Stability Analysis for the Case of Combination Resonances
In the present section, the approach of the non-linear stability can be used to study the mentioned dynamical system of the investigated problem. The frequency and the parameters of detuning 1 , 2 play a significance role on the destabilization of the stability criteria. The stability diagrams have been drawn for the system Eq. (50) using various selected values of these parameters. The changes of the amplitudes b 1 , b 2 via time are drawn for various parametric regions. The response of the frequency of them as a function of 1 is plotted in Fig. 31, in which the different values of f 2 (= 0.001, 0.0015, 0.002) have been taken into consideration. It can be seen that, there are two fixed points corresponding to every value of f 2 . The stability and instability areas of these points are classified in Table 4 and described according to: 1. At f 2 = 0.001 , the unstable fixed points have been explored in the range −1 ≤ 1 ≤ 0.16 and the stable ones have been detected in the range of 0.16 < 1 ≤ 0.6 , in which they are followed by unstable fixed points in the range 0.6 < 1 ≤ 1.
Based on Fig. 37 and Routh-Hurwitz stability conditions, there exists six possible intersection points, four of them are unstable and the others two points (1.58528, −0.660867) and (1.60119, −0.675496) are stable. The stable and unstable points are marked with red and green colours, respectively. Now, we will look at the properties of the non-linear amplitudes of system Eq. (50) and their stability according to the following transformation [40,41] Substituting Eqs. (13), (57) into Eqs. (22), (47) and then separating the real and imaginary parts to obtain (57) n j =ñ j , q j =q j (j = 1, 2).
The amended amplitudes were evaluated in various parameter regions over the specified time interval, and the amplitude characteristics were displayed in phase-plane curves as in Figs. 38 and 39 using the following data The impact of the various selected values of f 2 on the paths of n 1 and q 1 is drawn in Fig. 38. Quasi-periodic decelerating waves are portrayed in parts (a) and (b) of this figure, in which the number of oscillations increases to some extent with the increasing of f 2 values besides decreasing of the wavelengths. On the other hand, the depicted curves of Fig. 39 reflect the change of the modified phases n 2 and q 2 with time when f 2 has the above values as in parts (a) and (b), respectively. The behaviors of the represented waves of n 2 and q 2 have periodic attitude, in which the amplitudes of these waves increase with the increasing of f 2 values. Moreover, the projections of the trajectories of the ME on the phase planes n 1 q 1 and n 2 q 2 are represented in parts (c) of Figs. 38 and 39. A closed symmetric trajectories curves are observed which confirm our prediction for the steady motion of the studied system.

Conclusions
A dynamical system with 2DOF double RBP is investigated as a novel model. The motion of its suspension point is restricted to be along a Lissajous curve in the presence of two subjected external harmonic moments. The general EOM are derived using the second kind of Lagrange's equations and solved analytically up to the third-order of approximation using the MSM. The solvability conditions and the f 2 (= 0.001, 0.0015, 0.002), f 1 = 0.002, 1 = 0.11 , ME have been obtained in light of the categorized resonance cases. Two mainly external resonances are examined together in addition to a combination case of both the internal and the external resonance. The numerical results of the EOM have been obtained using Runge-Kutta method from fourth-order to draw the solutions. It has been found that the numerical results agree with analytical ones which show the efficiency of the used perturbation method. Moreover, the time-dependent of the gained analytic solutions, besides both of the resonance curves, amplitudes and modified phases are plotted. The stability areas of the fixed points at the steadystate have been examined in accordance with the non-linear analysis approach of Routh-Hurwitz. The properties of the modified amplitude of the non-linear analysis of the external resonance cases are examined and discussed. This model is very distinctive for its use in many engineering applications, for example, the vibrational of the motion of the robot as a dynamic system and analysis the control aspects of flexible arm robotics, pumps compressors, rotor dynamics, transportation devices, shipboard cranes, and human or robotic walking analysis.
Funding Open access funding provided by The Science, Technology & Innovation Funding Authority (STDF) in cooperation with The Egyptian Knowledge Bank (EKB). No funding agency in the governmental, commercial, or not-for-profit sectors provided a specific grant for this study.
Data Availability All data generated or analysed during this study are included in this published article.

Conflict of Interest
There is no conflict of interest declared by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.