Dynamics of a mechanical system with a spherical pendulum subjected to fractional damping: analytical analysis

In the paper, nonlinear vibrations of a system with three degrees of freedom having a spherical pendulum are considered. The system comprises a mass element suspended from a linear spring and a viscous damper, and a spherical pendulum swung from the mass element. It is assumed that the fractional viscous damping occurs in the viscous damper and at the pendulum pivot point. The viscoelastic properties of damping are assumed to be described using the Riemann–Liouville fractional derivative. The fractional derivative of an order of 0<α≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 0 < \alpha \le 1$$\end{document} is assumed. The nonlinear vibrations of the system near internal and external resonances are analyzed. The equations of motion of the analyzed system are solved using the multiple-scale method. The steady-state approximate solution is studied. The effect of a fractional-order derivative on the system vibrations is examined.


Introduction
The presented work is a continuation of previous works by authors dealing with the system of three degrees of freedom with a spherical pendulum [1][2][3]. In this work, a system consisting of a mass element and a spherical pendulum swung from the mass element is examined. The mass element is suspended from a linear spring and a damper. It is assumed that the damping in the system studied is described by the fractional Riemann-Liouville derivative [4] and this damping occurs in the damper and at the pivot point of the spherical pendulum. The considered system with a spherical pendulum can be used as a model of a real machine or its components, which operates in an energy-dissipating environment. In many scientific works, the systems containing a spherical pendulum are used to model the dynamics of certain types of structures, such as cranes [5][6][7][8][9][10], vibration absorbers [11][12][13], energy harvesters [14]. Thus, the dynamics of systems with a spherical pendulum is an absorbing issue of scientific research and has been studied in a number of researches [15]. A brief review of publications dealing with this issue is presented in the paper by Han et al. [16].
In previous works, the authors studied autoparametric systems containing a spherical pendulum with a viscous and magnetorheological energy dissipation system. Sado et al. [1] analyzed the influence of initial conditions on energy transfer between vibrating elements and the existence of chaotic motion in a system with a spherical pendulum. Sado and Freundlich [2] studied a dynamical behavior of a three-degree-offreedom system having a spherical pendulum. This system was controlled by a magnetorheological damper. The analyzed system consisted of a spherical pendulum suspended from a mass block which was suspended from a vertical linear spring and a magnetorheological damper. They investigated the influence of magnetorheological damper parameters on the system vibration close to internal and external resonances. These studies revealed that in addition to the regular behavior of the spherical pendulum, chaotic oscillations for all coordinates can arise near internal and external resonant regions. In all above-mentioned authors' works, the obtained equations of motion were solved using numerical methods.
It is well known that energy dissipation can have a significant impact on the dynamic behavior of a structure or its components. For this reason, various advanced methods for modeling damping in mechanical systems are being developed. One of these methods is the use of fractional derivatives to model energy dissipation. The use of fractional derivatives to model energy dissipation has increased very significantly over the past few decades [17][18][19]. Fractional derivatives are used to describe viscous damping because these derivatives allow a more accurate description of the damping phenomenon over a wider frequency range [20,21]. These derivatives have also been employed in modeling processes of energy dissipation in systems having pendulums [22][23][24][25]. Thus far, however, there has been small number of publications investigating dynamical systems with pendulum and fractional damping, especially with a spherical pendulum.
Rossikhin and Shitikova [22] investigated damped oscillations of two-degree-of-freedom system with a plane pendulum suspended from a mass element which was attached to a spring. They assumed that the system oscillates in a viscous medium whose damping properties are described by fractional derivatives. Additionally, the authors assumed small finite amplitudes of vibrations which allowed the use of multiple-scale method to solve the problem. They studied the impact of the damping described by the fractional derivative on damped free vibrations and the energy transfer in the system. Seredyńska and Hanyga [23] studied damped vibrations of a planar, inextensible and extensible pendulums in which the damping was described by a fractional derivative. This analysis was an example of the presented method for solving nonlinear differential equations with fractional damping. They determined the conditions of existence, uniqueness and dissipativity for a certain class of nonlinear dynamical systems including systems with fractional damping.
Hedrih [24] analyzed multi-pendulum systems with fractional-order creep elements. In this study, parallel pendulums were joined with creep elements, which were modeled using fractional-order derivatives. The governing equations of the system and its analytical solution for selected cases of the pendulum system were presented. The vibration modes of the systems with one and two pendulums having creep fractional elements were analyzed. The authors concluded that there is a mathematical analogy in descriptions between multipendulum systems and chain dynamical systems.
To the our knowledge, thus far only by the authors have performed the study of vibrations of a spherical pendulum with fractional damping [3]. In the aforementioned work, the authors assumed fractional damping only in the damper attached to the mass element [3]. The effect of a fractional-order derivative on the system vibrations was analyzed using numerical calculations. The impact of the fractional damping on the system vibrations waveforms and on the energy transfer in the system was shown. Thus, this study is a continuation of the authors' earlier work.

Description of the analyzed system
In this study, we consider a system with a spherical pendulum suspended from an oscillator excited harmonically by a force F z (t) = P 1 cos(ν 1 t) acting in the vertical direction ( Fig. 1). Additionally, the pendulum is excited harmonically in horizontal direction by forces F x (t) = P 2 cos(ν 2 t), F y (t) = P 3 cos(ν 2 t). The oscillator consists a linear spring and a fractional damper. Furthermore, it is assumed that there is also fractional damping in the pendulum pivot point. This damping is expressed by moments proportional to fractional derivative of order α. Thus, the analyzed system has three degrees of freedom. The motion of the spherical pendulum can be analyzed using various coordinate systems [9,[26][27][28][29]. In this study, the spherical coordinates presented by Leung and Kuang [9], and Aston [28] are employed to describe the motion of the pendulum. The following generalized coordinates z, θ , φ are assumed (Fig. 1). The position of the mass element  1 is defined by coordinate z, whereas the position of the pendulum of mass m 2 and length l is defined by the coordinates: z, θ , φ. The coordinate z is the vertical displacement of the body of mass m 1 measured from the static position of equilibrium. The angle θ is the angle between the vertical axis and the deflection of the pendulum in the plane xz. The angle φ is the angle between the deflections of the pendulum in the plane xz and pendulum. The selected coordinates are useful in the dynamic analysis of the spherical pendulum [28,29] and enable more interesting results to be obtained than using the classical spherical coordinates [1][2][3].
The position of a pendulum bob of mass m 2 in Cartesian coordinates is determined as follows (see Appendix A) where z st is the static deflection determined as follows where g is gravitational acceleration The kinetic energy T of the system can be expressed as The potential energy V of the system is expressed as In this study, a fractional damping characterized by the damping coefficient c α 1 and the order of the fractional derivative α 1 is assumed in the damper, while for the coordinates θ and φ the damping at the pendulum pivot point is described by the damping coefficient c α 2 and the order of the fractional derivative α 2 . Thus, the dissipation force R(ż (α) ), moments M(θ (α) ) and M(φ (α) ) ( Fig. 1) are determined by following expressions where z(t), θ(t) and φ(t) are the generalized coordinates, c α 1,2 are damping coefficients and d α dt α is a fractional derivative of the order α 1,2 .
In this analysis, the fractional Riemann-Liouville derivative [4] is used, which it is defined as where Γ (m − α) is the Euler gamma function [4], m is a positive integer number satisfying inequality m − 1 < α < m and t > 0. The fractional derivative order is assumed to be in a range of 0 < α ≤ 1.
Using the fractional dissipation function D = The dimensionless equations can be obtained by introducing the dimensionless time τ = ω 1 t and defining the following parameters Using parameters defined in Eq. (8), the equations of motion (7) can be transformed into a dimensionless form (where the overbars are omitted for convenience)

Method of solution
An approximate solution to Eq. (9) can be obtained using the multiple-scale method [31]. For small osculations in the vicinity of equilibrium position, trigonometrical functions can be expanded into Maclaurin series; thus, Substituting the approximated trigonometrical functions Eq. (10) into Eq. (9), the following system of equations is obtained The approximate solution of Eq. (11) for small vibrations can be expressed by expansion with different timescales as shown below [31][32][33] where T n = ε n τ (n = 0, 1, 2, 3 . . .) (13) are new independent variables, ε is a formal small parameter, T 0 = τ is the fast timescale and T 1 , T 2 are slow timescale [31,32].
Using the chain rule, the integer-order derivatives can be expanded in series of a small parameter ε where D n = ∂ ∂ T n A fractional-order derivative can be expanded in series of a small parameter as shown Rossikhin and Shitikova [18,34] are the Riemann-Liouville fractional derivatives with respect time T 0 .
Introducing additional small parameters [33], Substituting expression (12)-(16) into (11) and equating terms standing at the equal powers of a small parameter ε and limiting the approximate solution by terms of ε 2 , a following system of recurrent equations can be obtained for ε 1 Since further in the present analysis the expansions for displacements are limited by the expressions (18) of order 2 , we assume that the amplitudes A z1 , A θ1 and A φ1 are functions of time T 1 only. Therefore, the sought solutions to Eq. (17) are as below where A z1 , A θ1 and A φ1 are arbitrary complex functions of the timescale T 1 , and overbars denote complex conjugate functions.
In general, the Riemann-Liouville fractional derivative of the exponential function may be calculated according with method presented by Rossikhin and Shitikova [34,35], namely It can be shown that if the lower limit of the integral in the definition Eq. (6) is −∞ then the Riemann-Liouville fractional derivative of the exponential function has a form [4,35] D α + e iωt = (iω) α e iωt (21) where D α + is defined as [4,35] D α The improper integral in Eq. (20) may be omitted under certain circumstances, which are justified in the papers by Roshikhin and Shitikova [18,35]. Furthermore, Roshikhin and Shitikova [34] have shown that the improper integral in Eq. (20) does not affect the solution obtained by the method of multiple timescales when it is limited to the first-and second-order approximations. Thus, in further analysis the simplified derivative of the exponential function Eq. (21) is used.
Substituting solution to the fist approximation Eq. (19) into equations for the second approximation Eq. (18), the following system of equations is obtained where cc. stands for complex conjugate terms. Then, in order to eliminate the expressions that result in secular terms, we need to distinguish the following cases 2β = 1 , μ 1 = 1 , 1−β = β , μ 2 = β , μ 3 = β . In the analyzed system, the internal resonance occurs for β = 0.5, whereas the external resonances occur for μ 1 = 1 , μ 2 = β and μ 3 = β . All resonances should be analyzed separately. We are considering the internal resonance for β = 0.5 and the external resonance for μ 1 = 1 . Introducing detuning parameters σ 1 and σ 2 , we assume that p 2 = p 3 = 0 and The secular terms in Eq. (23) may be eliminated if Assuming that Noting that the amplitudes a z1 , a θ1 and a φ1 are the functions of time T 1 , and considering that T 1 = εT 0 , then substituting expressions (26) into system of Eq. (25) and separating real and imaginary parts of Eq. (25), we obtain the following equations −a z1 ψ 1 − 1 2 γ 1 a z1 cos πα 1 2 where We assume for steady-state solution that Noticing that and substituting Eqs. (29) and (30) into Eqs. (27) and (27) have the form It can be concluded that according to the assumptions made previously, four cases of the steady-state solution are possible.

4.1
The case with amplitudes a θ1 = 0 and a φ1 = 0 The first case is if amplitudes a θ1 = 0 and a φ1 = 0, in this case the pendulum does not vibrate. The system corresponds to a one-degree-of-freedom oscillator with the mass m = m 1 + m 2 and the amplitude a z1 is a z1 γ 1 cos Thus, solving Eq. (32), the amplitude a z1 and the phase angle Θ 3 can be calculated, namely Equation (33) shows that for parameter σ 1 = 0, amplitude a z1 does not depend on the order of the fractional derivative α 1 . This dependency for small values of damping coefficient γ 1 and for σ 1 = 0 is weak.

4.2
The case with amplitudes a θ1 = 0 and a φ1 = 0 The second case is if the amplitudes a θ1 = 0 and a φ1 = 0, thus a following system of equations may be formulated Solving Eq. (35), the amplitude a z1 and the phase angle Θ 2 can be obtained, viz Similarly, the amplitude a φ1 and the phase angle Θ 3 can be calculated using the equations (35) Equation (36) shows that if σ 1 = −σ 2, amplitude a z1 does not depend on the order of the fractional derivative α 2 . This dependency for σ 1 = −σ 2 and for small values of damping coefficient γ 2 is weak.

4.3
The case with amplitudes a θ1 = 0 and a φ1 = 0 The third case is if the amplitude a θ1 = 0 and a φ1 = 0, then The amplitude a z1 and the phase angle Θ 1 can be calculated, namely The amplitude a θ1 and the phase angle Θ 3 can be calculated using Eq. (38) Similarly as in the previous case, if σ 1 = −σ 2, amplitude a z1 does not depend on the order of the fractional derivative α 2 . This dependency for σ 1 = −σ 2 and for small values of damping coefficient γ 2 is weak.

4.4
The case with amplitudes a θ1 = 0 and a φ1 = 0 The forth case is if the amplitudes a θ1 = 0 and a φ1 = 0, then solving Eq. (29), amplitude a z1 and angle Θ 2 can be derived, namely Equation (29) shows that cos Θ 2 = cos Θ 1 and sin Θ 2 = sin Θ 1 . Taking this into account, after some mathematical transformations, we can find that Having expression for a z1 and tan (Θ 2 ), we can derive expressions for amplitudes a θ1 and a φ1 by solving Eq. (43). In this case, amplitude a z1 does not depend on the order of the fractional derivative α 2 . This dependency for σ 1 = −σ 2 and for small values of damping coefficient γ 2 is weak.

A case of the internal resonance for β = 0.and the external resonance for μ 2 = β
We are considering the internal resonance for β = 0.5 and the external resonance for μ 2 = β . Introducing detuning parameters σ 1 and σ 3 , and assuming that p 1 = 0 whereas p 2 = 0 and p 3 = 0 The secular terms in Eq. (23) may be eliminated if Assuming that the amplitudes have the same form as in section 4 Eq. (26) and noting that the amplitudes a z1 , a θ1 and a φ1 are functions of time T 1 , and considering that T 1 = εT 0 , then substituting expressions (26) into system of Eq. (45), and separating real and imaginary parts of Eq. (45), the following equations can be obtained where We assume for steady-state solution that a z1 = 0, a θ1 = 0, a φ1 = 0, Noticing that and substituting Eqs. (48) and (49) into Eq. (46), we obtain following expressions Equation (50) shows that several cases of the steadystate solution are possible. The first case is if amplitude a z1 = 0, in this case the mass element does not vibrate. Amplitudes a θ1 , a φ1 are equal and The second case is when all amplitudes are not equal zero. In this case, the relationship between amplitude a θ1 and a z1 is as below where w θ = a z1 cos (Φ 1 ) + 2 γ 2 β cos πα 2 2 The relationship between amplitude a φ1 and a z1 is expressed as where The next two cases occur when one of the amplitudes a θ1 = a φ1 = 0. Equation (50) shows that this case is possible when amplitudes p 2 = p 3 = 0; thus, the pendulum does not vibrate. Another cases occur when a θ1 = 0 and a φ1 = 0 or a θ1 = 0 and a φ1 = 0. In these cases, one of the amplitudes p 2 = 0 or p 3 = 0 correspondingly, and the movement of the pendulum is in one plane.
The dynamic behavior of the system can be also analyzed when one of the forces p 2 , p 3 is zero. For example, if p 2 = 0 then Eq. (50) shows that the amplitude a z1 is expressed as The amplitude a φ1 can be calculated using Eq. (53) and phase angle Φ 1 may be calculated from Having calculated a z1 , a φ1 and Φ 1 , amplitude a φ1 and phase angle Φ 2 can be determined using Eq. (50). Equation (54) shows that amplitude a z1 does not depend on damping coefficient γ 1 but it depends on damping coefficient γ 2 and detuning parameters σ 3 . If σ 3 = 0, amplitude a z1 does not depend on the order of the fractional derivative α 2 . This relationship for σ 3 = 0 and for small values of damping coefficient γ 2 is weak.
The graph shown in Fig. 3 shows that the amplitude a θ1 also depends weakly on the order of the fractional derivative, virtually all curves coincide. The graphs presented in Fig. 3 show that the order of the fractional derivative effects on the range of existing real solution of Eq. (43) for amplitude a θ1 , namely an increase in the order of the fractional derivative, decreases the range of the solution. The decrease in the range of the solution is more noticeable for a higher damping coefficient γ 2 .

Conclusions
In this paper, analysis of a nonlinear three-degreeof-freedom system with a spherical pendulum is performed. A fractional damping is assumed in the damper and in the spherical pendulum pivot point. The approximate analytical solution is obtained using the multiplescale method. Steady-state solutions for different combinations of external and internal resonances are studied. The analysis is performed for two types of excitation. The first excitation case assumes only excitation by a vertical force acting on the oscillator and for an internal resonance for β = 0.5 and an external resonance for μ 1 = 1.0, while the second case assumes excitation with a force acting on the pendulum in the horizontal direction (Fig. 1) and an internal resonance for β = 0.5 and an external resonance for μ 2 = β.
It is shown that the amplitude a z1 depends weakly on the order of the fractional derivative, as well on the damping coefficient. Similarly, the amplitudes a θ1 and a φ1 depend weakly on the order of the fractional derivative. The study can be extended to a transient analysis and analysis of other external and internal resonances.
Acknowledgements The authors would like to express their thanks to the Organizing Committee of the 16 th International Conference "Dynamical Systems -Theory and Applications, December 6-9, 2021" for recommending this manuscript for publication in "Nonlinear Dynamics." Funding The authors have not disclosed any funding.

Data availability
The data generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of interest
The authors declare that they have no conflict of interest.
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://creativecommons.org/licenses/ by/4.0/.

Appendix A
Relationship between Cartesian and generalized coordinates Eq. (1) can be obtained by analyzing the appropriate distances and the triangles OC B and AB O shown in Fig. 4. The triangle OC B lies in the plane Ox 2 z 2 . The triangle AB O lies in the plane perpendicular to the plane Ox 2 z 2 and passing through the pendulum O A; thus, the segment AB is perpendicular to the segment O B and the angle at the vertex B is a right angle. The angle θ lies in the plane Ox 2 z 2 and it is between the vertical axis Oz 2 and the orthogonal projection of the pendulum O A on the plane Ox 2 z 2 , i.e., the segment O B. The angle φ is the angle between the deflections of the pendulum in the plane x 2 z 2 and the pendulum [9].
It can be seen from the triangle AB O (Fig. 4) that The triangle OC B (Fig. 4) shows that