Modeling and Dynamic Analysis for a Motion of Mounted-Based Axisymmetric Rigid Body under Self-Excited Vibrations in an Attractive Newtonian Field

In this article, the mathematical modeling and dynamic analysis for a motion of a rigid body mounted on a vibrating base affected by a strongly nonlinear damping and an attractive Newtonian field are investigated. *e governing equations are derived, and specific conditions for asymptotic stability of the motion are stated. Results of analytical analysis are used to identify some distinct types of motion, including the cases that are apparently regular or chaotic. A reduction technique via averaging is used to obtain the amplitude equation and the response diagrams to identify regions where the jump phenomenon may occur. Moreover, a simple analytical relationship between the parameters of the system, describing whether the jump phenomenon will be possible, is obtained. Homoclinic bifurcation diagrams are used to argue apparently that the chaotic behavior is slowly approaching a strange attractor at specific values of system parameters. Results of numerical solutions using the fourth-order Runge-Kutta method are closely coincided with the analytical ones to verify all types of motion.


Introduction
In engineering applications, the vibrational motion is typically produced internally from sources of noise such as motors, bearings, and other moving parts in machines. erefore, it is necessary to take into account the most significant destabilizing source that can seriously degrade the operation and, in some cases, it leads to chaotic behavior attached by a catastrophic failure of mechanical system devices. Generally, the vibration is described to characterize the machines dynamic force in terms of a single dominant frequency called the operating frequency. is particular modeling, indeed, contains a generic assumption; namely, the excitation force is ideal, which means that the dynamics of excitations are not affected by the dynamics of responses. is assumption becomes more untenable as the operating frequency increases further into the audio or noise ranges, cf. [1].
Clearly, the sources of vibrations output in machines are considered due to the presence of several moving parts such as reciprocating pistons and rotating shafts. For instance, reciprocating pumps and washing machines both generate sinusoidal vibration type due to the reciprocating primary motion or the primary rotational motion, respectively. Indeed, this phenomenon generates dynamic forces which might be sinusoidal, random, or transient. erefore, adequately, these machines can be modeled by the motion of a rigid body mounted on a vibrating base with the presence of dampers, cf. [2].
In rigid body dynamics, the recent studies have focused on the damped motion of self-excited rigid bodies. is type of motion is described as a body free to rotate about a fixed point in a resistant medium acted upon a time-dependent torque vector. e external torques are often arising from internal reactions which do not appreciably change the distribution of their mass, cf. [3][4][5][6]. e interesting engineering examples of this motion are gyroscopes and dynamics of rotating missiles. However, a modern case study of a self-excited motion of a rigid body can arise from the attitude motion of a rigid spacecraft subject to thruster torques, cf. [7][8][9].
Due to the nature of the vibrational motion, only a few mechanical systems are dominated by motion in one dimension and thus single degree of freedom (SDOF) models can be used. But, in general, most of the real mounted-based systems are always modeled as a piece of equipment vibrating about its center of mass in two or more of degrees of freedom (DOFs). So, as consequences in the presence of all six modes' modeling, they are coupled resulting in six coupled equations to describe such motion, cf. [9,10].
Self-excited vibrational systems might be characterized by a principal feature which is resided in the damping behavior due to the dissipative forces.
us, even an infinitesimal disturbance causes sustained oscillations in the system; however, when the perturbed displacement becomes sufficiently large, the damping becomes positive and suppresses further increase in its amplitude.
In this work, the following specific form of damping function is used, which is considered a generalized form of van der Pol's damping.
Using the dimensionless small parameter (ε), the concept of homoclinic bifurcation and the route to chaos are necessary in the study due to the severe existence in complex systems, in particular, with strong damping forms. Homoclinic bifurcation takes place due to a parameter when the change of the parameter causes a homoclinic path to appear at a certain value and then it disappears again, cf. [11][12][13]. e generation of chaotic set can be arisen by such parameter to bifurcate at an intermediate value into a solution which circuits the origin point (for example, say twice). Further increasing shows further bifurcation has taken place which results in paths circuiting the origin many times. If the process continues with further period-doubling circuits, then the parameter sequence converges to a certain value beyond which all these periods are presented consisting of a chaotic attracting set. en, this configuration is known as a strange attractor which is developed parametrically through processes of the period-doubling, cf. [14,15]. Consequently, at each bifurcation, there exists an unstable periodic orbit that will be always dwelled in, so that the final attractor has to include unbounded sets of unstable limit cycles. is means that the periodic solution becomes closer to the destabilizing effect of the homoclinic paths which appear as the parameter increases. en, now, the chaos is generated by the effect of that parameter, cf. [16][17][18].
In common, two main methods are used to measure the existence of chaotic behavior in nonlinear (nonautonomous) systems: Lyapunov exponents and Poincare's maps. e method of Lyapunov exponents measures the exponential rates of divergence (convergence) of nearby orbits of an attractor in the state space. erefore, it is the most precisely powerful tool for the identification of motion in dynamical systems. In [19], an algorithm is developed based on tracing the evolution of an initial sphere of small perturbation to a nominal trajectory. In [20,21], a method to compute Lyapunov exponents from an experimental time series based on the QR methods is presented.
In this work, we study the dynamics of a rigid body mounted on a vibrating base with a generalized nutational damping, using simulation, numerical continuation, and analytical stability theory. We begin by stating the equations of motion, the stability conditions in terms of problem parameters, and then theorems to predict the existence of periodic solutions. Analytical methods and numerical simulation results are used to identify the distinct types of motion, including persistent oscillations that appear to be a stable, periodic, or chaotic in their nature. e theoretical results will be applied to the tackled application beside the study of bifurcation and routes to chaos of the governing system. Melnikov's method is used to transpire that the ranges of chaotic behavior are affected by the change of some certain parameters. e verification and compatibility of theoretical results with the numerical ones through bifurcation and chaos diagrams as well as phase plane trajectories are taken into consideration in the study. e rest of the article is summarized as follows: Section 2 presents the modeling of a damped motion of mountedbased rigid body as an SDOF model in attractive Newtonian field. Section 3 discusses the study of stability analysis in both autonomous and nonautonomous cases of motion. Section 4 discusses the existence of periodic solutions for some specific cases of motion. Section 5 presents a reduction technique via averaging to capture the jump phenomenon. Sections 6 and 7 discuss the study of homoclinic bifurcation and routes to chaos. In the last section, the conclusion is given.

The Modeling and Governing Equations
Consider that a rigid body of mass (m) is mounted on a vibrating base of a sinusoidal excitation frequency (ω) as shown in Figure 1 with two sets of coordinates: one is fixed in space and the other is moving with the body. e two coordinates are coincident when the body is in equilibrium under the action of gravity force alone. e motion of the body is described by giving it a displacement relative to the inertial fixed axes. Let us choose the axes OX, OY, and OZ to represent the fixed frame in space with unit vectors I, J, and K, respectively, and the axes ox, oy, and oz to represent the principal axes constructed on the body at the fixed point (o) with unit vectors i, j, and k, respectively.
Hence, the translational displacement of the center of mass (C.G.) with respect to the fixed axes (X, Y, Z) reads Mathematical Problems in Engineering en, we obtain the velocity of the center of mass: We assume that the principal axes rotate in space with the same angular velocity of the body (ϖ) defined by (in terms of Euler's angles) where θ(t), ϕ(t), and ψ(t) are the nutation, spin, and precession angles, respectively, and (·) refers to the derivative with respect to time.
Let K(t) be the Poisson vector with the following direction cosines with respect to the moving axes: e relations of the direction cosines with Euler's angles read e potential energy of a Newtonian force field is considered by the Newtonian attraction of a point P, of mass M, on the body mass (m). If P is located on the negative part of the axis OZ, then we have here, G is the universal gravitation constant and R is the distance between point P and the fixed point (O). In general, the fact that R is larger than the dimensions of the body is taken into account, cf. [6]. Hence, the potential energy reads where here, A, B, and C are the principal moments of inertia along the moving axes ox, oy, and oz, respectively, and J (...) is the third moment of inertia of the body. e total potential energy (U) of the mass with respect to the inertial frame OXYZ is Consider the following approximation: and then U � −mgR + mg x c sin θ sin ϕ + y c sin θ cos ϕ + z c cos θ − 3g 2R A sin 2 θ sin 2 ϕ + B sin 2 θ cos 2 ϕ + C cos 2 θ e kinetic energy reads en, the system admits the following Lagrangian (L) on the Lie group SO (3): where I c � I 1 , I 2 , I 3 refers to the principal moments of inertia about the equatorial and polar directions through the C.G. of the rigid body. For simplicity, the mass of rigid body satisfies the symmetric conditions and the excitation displacement vector is only on the fixed vertical direction as follows: en, where d(t) � d 0 sin ωt is the time-varying amplitude of the vertical support motion along the vertical fixed axis (OZ) with constant amplitude d o and forced frequency ω. Using Routhian's function (R), where p ϕ and p ψ are the corresponding generalized momenta to the cyclic coordinates ϕ and ψ, respectively, cf. [8,9,22,23]: and then the equation of motion reads )) + F 1 cos ωt + F 2 sin ωt. (24) e right terms in equation (24) represent the damping term (D) and the harmonic forcing terms of amplitudes F 1 and F 2 > 0 with circular frequency (ω) applied to the body on such nonlinear spring. e damping term (D) is assumed to have the following generalized van der Pol's damping form: Using equation (24) and if p ψ � p ϕ � p o , the governing equation becomes Equation (26) can be generalized as follows: where Consider the following notations: For simplicity, let us assume that en, equation (27) reads If θ � x and θ . � y, then equation (32) is converted to the following system: Use the following approximations: and then equation (32) reads and the approximate system of equation (35) reads (36b) If θ ≪ 1, then equation (35) has the following linear form under a proposed fractional damping: where L D q θ is Liouville's fractional derivative defined as (38) where Γ(.) is Euler's Gamma function and n − 1 < q < n, n � 1, 2, 3, . . .. Liouville's fractional derivative possesses the following fractional derivative properties:

Stability Analysis
e importance of stability analysis in dynamical systems lies in understanding how changes in design parameters might cause bounded or unbounded responses. Two methods are proposed to study the stability for equation (32): linear analysis to obtain the transition curves of the system taking into account the fractional derivative of damping and the nonlinear analysis to obtain the general conditions of the stability on the whole, cf. [24][25][26][27][28][29][30].

Linear Analysis.
e linear approximation of the nonlinear system by linearizing it at the fixed points is mostly useful technique to predict the geometrical nature of the equilibria and the prediction of stability domains separated by boundary curves. In general, equation (32) has the following equilibria for all parameter values: (i) (x * , y * ) � (0, 0), which describes a motion in which one of the body principal axes coincides with the local vertical axis (ii) (x * , y * ) � ( ± π, 0) which describes an inverted motion of the rigid body us, the linearized homogenous form of equation (37) at a general equilibrium point θ * � 0 π or -π for the variable en, consider that Hence, Due to the fact that the fractional derivative is a nonlocal operator, the choice of Liouville's definition allows us to obtain the periodic solution based on the results provided in [31,32]. Since the fractional derivative typically depends on the values between the lower limit t o and the point of interest (t), in this case, we are interested in looking for periodic solutions of equation (44) However, we look for the regions of stability using the harmonic balance method for the homogenous linear simplified form of equation (44) under the action of Liouville's fractional derivative. From Floquet theory, it is clear that this equation admits solutions of period π or 2π, cf. [33,34].
Typically, for given values of the parameters δ and ϵ, all solutions of this equation should be either stable, unstable, or periodic. e (δ, ϵ) plane (Ince's diagram) is normally divided into stable and unstable regions separated by the transition curves to represent the periodic solution curves. e aim here is to present these transition curves in order to be able to assess the influence of the fractional derivative term on the system. Now, we seek the stability domains of equation (44) separated by the transition curves using the harmonic balance method. is method was firstly presented to specify the periodic solution curves by employing operational matrices associated with a Fourier basis; and, as a result, the conditions for the existence of nontrivial periodic solution were obtained by setting Hill's determinant to zero. So, we Mathematical Problems in Engineering consider a formal solution in the following Fourier expansion: en, inserting equation (45) into equation (44) and collecting terms, we obtain the following recurrence relations: As shown in Figure 2, the transition curves of the fractional damped homogenous linear part of equation (32) are obtained, as well as subsequently the stability regions. It is vividly shown that a slight change in the fractional order (q) or the damping coefficient (ζ) affects the shape and location of the transition curves. Clearly, this effect can be characterized by the location of the lowest point on the transition curve, which specifies the minimum value of forcing amplitude (ϵ), necessary to produce unstable domains in the linearized form. It is clear that the location of the minimum point on the transition curve due to the fractional damping moves this lowest point in a horizontal direction, but due to the damping coefficient it moves in the vertical direction until the unstable region disappears. ese results are completely coincided with the fractional parametric excitation of Mathieu's equation in [35] and Ince's equation in [36].

Nonlinear Analysis.
In this part, we consider the study of stability on the whole for the governing equation (equation (32)) for both autonomous and nonautonomous forms.
Then, the sufficient condition for a stable motion is Proof. e autonomous form of equation (32), when c � 0, is equivalent to We assume that By taking Lyapunov function in the form we obtain Hence, we have the following conditions: ii. (H(x, y) − H(x, 0))y < 0 for y ≠ 0, iii. −   (32): Then, the sufficient conditions for a stable motion are Proof. e form of equation (60) is equivalent to the system where, By taking Lyapunov function in the form we obtain So, every solution x(t) exists in the future and x(t) ⟶ 0 and y(t) ⟶ 0 as t ⟶ ∞ if the following conditions are satisfied: i .h(x, y) is continuous and nonnegative on, R × R and is bounded when Apply the conditions; then, the conclusion holds. □ Theorem 3. e governing equation, equation (27), has an asymptotically stable solution if Proof. By considering that in the domain k ∈ R + , 0 ≤ t ≤ T, |x| < ∞, |y| < ∞, and x 2 + y 2 ≥ k 2 . Hence, we have en, it is easy to obtain the general conditions:

Mathematical Problems in Engineering (i) h(x, y) and f(x, t) are locally Lipschitzian in x and y, respectively, and e(t) is continuous on R (ii) f(x, t) and e(t) are periodic in t of period
Apply the conditions; then, the conclusion holds. □

Existence of Periodic Solutions
Let us assume that equation (32) has a periodic solution with period T � (2π/ω). Also, it is considered that, at the equilibrium point θ * � 0, the following assumption holds: Hence, the only solution of equation (72) is θ � θ * . Equation (32) is rewritten as follows (θ � x, θ . � y): where the RHS represents the absorption h(x, y) and the supply (e(t)) of the system. us, the following theorem gives for equation (73) the general conditions for the existence of periodic solutions.

Theorem 4. Equation (73) has at least one periodic solution of period (2π/ω) if (i) h(x, y) and f(x, t) are locally Lipschitzian in y and
x, respectively, and e(t) is continuous on R + (ii) f(x, t) and e(t) are periodic in t of period Proof. For the domain t ∈ [0, 2π/ω], let |x| < ∞, |y| < ∞, and x 2 + y 2 ≥ k 2 and k ∈ R + . Consider that there exists a continuous function such that Assume that the following identities hold: and then let the two functions V 1 and V 2 be as follows: Hence, it can be concluded that there exist in the interior of two domains en, the solutions x(t) and y(t) exist and are bounded for all t ≥ 0. en, the conclusion holds.
Proof. By using the polar coordinates one gets the following: After substitutions and using the conditions ω n � 1 and Ω � 0, we obtain and one particular solution is Corresponding to the limit cycle, It is noticed that is means that all paths approach the shown radius of the limit cycle ρ o especially when μ � ]. en, the conclusion holds.
Mathematical Problems in Engineering

Averaging Analysis
We assume that the solution of the approximate form of governing equation, equation (35), is still given by With time-varying coefficients a(t) and b(t), in other words, we consider equation (85) as a transformation from a(t) and b(t). Hence, by imposing a third condition or a third equation independent of equation (35) and equation (85), assuming that x(t) and _ x � y(t) have the same form as the linear case, and the following acceleration term is obtained: Consider the following condition: en, we have to transform the system equation (85) to equation (88) with equation (35) into the two unknowns a(t) and b(t) as follows: where Γa a 2 − 3b 2 + higher harmonics.

(90)
Consequently, by considering that the major parts of a and b are slowly varying functions of time, they change very little during the time period T � (2π/ω) and, in addition to the first approximation, they can be considered constant in the interval [0, T]. Hence, we average both sides of the system equation (89) for a and b over the interval [0, T] and obtain Consequently, the initial conditions can be given in terms of those for equation (35) by (94)

Relative Equilibria and Stability.
According to equation (93), it can be deduced that all equilibrium solutions are given by _ a � _ b � 0. In addition to the fact that the other paths corresponding to solutions of equation (35) are nonperiodic in general, let (a o , b o ) represent any of the three fixed points, (a 1 , b 1 ), (a 2 , b 2 ) and (a 3 , b 3 ), of equation (93). Hence, the linear analysis in the neighborhood of point and then the corresponding linearized system is where e general conditions of stability for every stable fixed point are Moreover, the boundary between nodes and spirals can be given by the following identity: e phase diagram computed for a particular case ε � 0 is shown in Figure 3 for F � 0 and F � 0.005, respectively. In Figure 3(a), point P 1 is a center and represents a stable oscillation, but in Figure 3(b) there exist three equilibrium points: P 1 and P 2 are centers and P 3 is a saddle representing an unstable oscillation.
For the particular case of ε � 0.005 and λ � μ � ] � 1, the phase plane is shown in Figure 4 for F � 0 and F � 0.005, respectively. In Figure 4(a), point P 1 is a spiral point and represents a stable oscillation, but in Figure 4(b) there exist three equilibrium points: P 1 and P 2 are stable spiral points and P 3 is a saddle representing an unstable oscillation. e phase diagram computed for a particular case where ε � 0.005, λ � −1, and μ � ] � 1 is shown in Figure 5 for F � 0 and F � 0.005, respectively. In Figure 5(a), point P 1 is an unstable point as it bifurcates to generate a stable limit cycle, but in Figure 5(b) there exist three equilibrium points: P 1 and P 2 are centers and P 3 is a saddle representing an unstable oscillation. (93), three are three different amplitudes at the equilibria. Consequently, due to the existence of the unstable node, the amplitude surface (ω, F, ρ) has unstable oscillations. en, the surface exhibits a fold catastrophe or a cusp which leads to the jump phenomenon, cf. [14]. Illustratively, the cubic amplitude equation in (ϱ) representing the surface of equilibrium points (Ξ(a, b) � Λ(a, b) � 0) is given by

e Jump Phenomenon. From equation
where Hence, there may be as many as three positive values of ϱ satisfying the equation of amplitude. is explicates that, for certain ranges of the problem parameters, there might be three distinct periodic solutions of the governing equation. Also, the responses are stable when (dF/dρ) > 0 and unstable when (dF/dρ) < 0.
Consider the following: Since G(0) � 0 and G(ϱ) ⟶ ∞ as ϱ ⟶ ∞, the amplitude equation must have at least one positive root. For some parameters' values, the amplitude equation has three different real roots corresponding to its cubic form of ϱ if has two distinct solutions; thus, So the solution of the amplitude equation is interested in ϱ 1 , ϱ 2 > 0; then, us, the amplitude equation has three distinct real roots if F is controlled by the overlap region.
To explain, we take the case of undamped motion ε � 0; then, we have and then, as shown in Figure 6, the overlap region is controlled by ��� 155 162 Hence, the cusps are exactly determined at In Figure 7, it is clear that the surface under the fold gets along with unstable oscillations which will not be obtainable or accessible using equation (100). e surface displays what is known as a fold catastrophe which produces the jump phenomenon controlled by parameter ε at the values of 0, 0.1, and 0.25. It is clear that, with the increase of parameter ε at the constant values of damping parameters λ, μ, and ], the jump is going to disappear.
As shown in Figure 7, it can be predicted that, for certain values of ω and F at the lower values of ε ≈ 0, there are three responses: the oscillations with greatest and smallest amplitudes are stable but the middle one is unstable. But when ε is increased, then there exists only a single stable response. e consequences of this fold are shown illustratively using Figure 8 and 9. So, consider that F is held constant at a value which intersects the fold, and the applied frequency (ω) is varied in such a way that it crosses the edges of the fold. At the small values of ε ≈ 0, as ω increases beyond, after some irregular motion, the amplitude jumps at the critical value; henceforth, it follows the smooth curve. On the way back, the response point moves along the upper curve. en, here it must go down to the lower curve and then back to the upper one crossing the jumping point of the amplitude.

Homoclinic Bifurcation
e existence of homoclinic bifurcation with the support of the other problem parameters, particularly α and β, is considered. e homoclinic bifurcation problem in equation (35) can be considered at certain values of the governing problem parameters. In general, Melnikov's method can be used to predict the critical value to the homoclinic bifurcation, but in most of the cases the perturbation method has a less accuracy when the perturbation parameter increases further, cf. [37,38]. Some critical values are detected for some special perturbed systems as shown in Table 1 at some specific values of the problem parameters, cf. [2,39,40].
To verify the theoretical results, system (xi) is considered as depicted in Table 1. System (xi) has the saddle point at (0, 0) with an associated homoclinic trajectory with solution (x 0 (t)) at ε � 0. So, using Melnikov's function, the homoclinic bifurcation occurs approximately where From Figure 10, it is shown that there exist three equilibrium points which become stable spirals if λ > λ c and unstable if λ < λ c . As λ touches the critical value (λ c � 0.81714) the path changes from heteroclinic spiral   saddle connection (λ > λ c ) to a homoclinic saddle connection at the critical value to heteroclinic saddle connection (λ < λ c ). is transition is due to the fact that parameter λ generates a homoclinic bifurcation at λ c .
A simple numerical search reveals a stable periodic solution which is shown in Figure 11. Hence, if c is increased further, the solution bifurcates at a specific value of c into a solution which wraps around the origin twice. Further increase of c leads to further bifurcation that takes place giving more periods of solution as vividly shown. e process continues with further period-doubling to have a chaotic attracting set which is known as a strange attractor. At each bifurcation, an unstable periodic orbit will be always resided in, so that the strange attractor must contain a bounded set of unstable limit cycles. For its final stage, the periodic solution becomes closer to the destabilizing effect of the homoclinic paths which appears as c increases. is analysis is confirmed by computing Poincare's map and Lyapunov exponents to detect period-doubling and chaos. From Figure 12, Poincare's map becomes bounded sets of returns without any obvious repetitions implying the existence of strange attractor. In this case, at least one of Lyapunov exponents should be positive and the total sum of Lyapunov exponents should be negative, which are typically satisfied in Figure 12. e solutions' figures as depicted show wandering solutions of an irregular oscillations type without any uniform pattern. en, it is inferred that such solutions are said to display a chaotic behavior.

Conclusion
is work is concerned with modeling of a motion of an axisymmetric rigid body under an external sinusoidal excitation and a Newtonian attractive field with a general van der Pol damping based on the nutational angle of Euler's interpretation. Two approaches are established to study the stability analysis: e first presents the linear analysis to obtain the transition curves of the system considering the fractional derivative of the damping term using Liouville's definition. e effect of the fractional order and the system parameters on the stability of the solution are presented using Ince's diagram. e second approach discussed the nonlinear analysis of the system. Sufficient conditions that guarantee stable motion of the system are established for both autonomous and nonautonomous cases. Consequently, a sufficient set of conditions that ensure the asymptotic stability of the solution are given. Also, other sets of conditions that ensure the existence of at least one periodic solution and a limit cycle are derived. In addition, averaging analysis to clear the relative equilibria with their stability and the description of the jump phenomenon are shown. e existence of homoclinic bifurcation based on system parameters using Melnikov's function is illustrated. e transition from homoclinic bifurcation to chaos is identified analytically and numerically using phase trajectories, Poincare's map, and Lyapunov exponents.

Data Availability
e data used to support the findings of this study are available from the authors upon request.