Abstract
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. The 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.
1. 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. Therefore, 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. This 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. This 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. Therefore, 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. This 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. The external torques are often arising from internal reactions which do not appreciably change the distribution of their mass, cf. [3–6]. The 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–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. Thus, 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.where the parameters , , and are real constants; .
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–13].
The 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. Then, 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. This means that the periodic solution becomes closer to the destabilizing effect of the homoclinic paths which appear as the parameter increases. Then, now, the chaos is generated by the effect of that parameter, cf. [16–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. The method of Lyapunov exponents measures the exponential rates of divergence (convergence) of nearby orbits of an attractor in the state space. Therefore, 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. The 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. The 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.
The rest of the article is summarized as follows: Section 2 presents the modeling of a damped motion of mounted-based 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.
2. The Modeling and Governing Equations
Consider that a rigid body of mass 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. The two coordinates are coincident when the body is in equilibrium under the action of gravity force alone. The motion of the body is described by giving it a displacement relative to the inertial fixed axes. Let us choose the axes , , and to represent the fixed frame in space with unit vectors , , and , respectively, and the axes , , and to represent the principal axes constructed on the body at the fixed point with unit vectors , , and , respectively.
Hence, the translational displacement of the center of mass (C.G.) with respect to the fixed axes readshere, is the excitation displacement vector. Then,
Then, 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 , , and are the nutation, spin, and precession angles, respectively, and refers to the derivative with respect to time.
Let be the Poisson vector with the following direction cosines with respect to the moving axes:
The relations of the direction cosines with Euler’s angles read
The potential energy of a Newtonian force field is considered by the Newtonian attraction of a point , of mass , on the body mass . If is located on the negative part of the axis , then we havehere, is the universal gravitation constant and is the distance between point and the fixed point . In general, the fact that is larger than the dimensions of the body is taken into account, cf. [6]. Hence, the potential energy readswherehere, , , and are the principal moments of inertia along the moving axes , , and , respectively, and is the third moment of inertia of the body. The total potential energy of the mass with respect to the inertial frame is
Consider the following approximation:
and then
The kinetic energy reads
Then, the system admits the following Lagrangian on the Lie group SO (3):where 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:
Then,where is the time-varying amplitude of the vertical support motion along the vertical fixed axis with constant amplitude and forced frequency .
Using Routhian’s function ,where and are the corresponding generalized momenta to the cyclic coordinates and , respectively, cf. [8, 9, 22, 23]:
and then the equation of motion reads
The right terms in equation (24) represent the damping term and the harmonic forcing terms of amplitudes and with circular frequency () applied to the body on such nonlinear spring.
The damping term is assumed to have the following generalized van der Pol’s damping form:
Using equation (24) and if , the governing equation becomes
Equation (26) can be generalized as follows:where
Consider the following notations:
For simplicity, let us assume that
Then, equation (27) reads
If and , 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
If , then equation (35) has the following linear form under a proposed fractional damping:where is Liouville’s fractional derivative defined aswhere is Euler’s Gamma function and , . Liouville’s fractional derivative possesses the following fractional derivative properties:
3. Stability Analysis
The 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–30].
3.1. Linear Analysis
The 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), which describes a motion in which one of the body principal axes coincides with the local vertical axis(ii) which describes an inverted motion of the rigid body
Thus, the linearized homogenous form of equation (37) at a general equilibrium point = 0 or - for the variable with the time scale ; then it reads
Then, consider that
Hence,
Then,
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 and the point of interest (t), in this case, we are interested in looking for periodic solutions of equation (44) in domain.
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 , cf. [33, 34].
Typically, for given values of the parameters and , all solutions of this equation should be either stable, unstable, or periodic. The plane (Ince’s diagram) is normally divided into stable and unstable regions separated by the transition curves to represent the periodic solution curves. The 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. This 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 consider a formal solution in the following Fourier expansion:
Then, 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 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. These results are completely coincided with the fractional parametric excitation of Mathieu’s equation in [35] and Ince’s equation in [36].
(a)
(b)
3.2. 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.
Theorem 1. Consider the autonomous form of equation (32).Then, the sufficient condition for a stable motion is
Proof. The autonomous form of equation (32), when , is equivalent toWe assume thatBy taking Lyapunov function in the formwe obtainHence, we have the following conditions:Apply the conditions on ; then, the conclusion holds.
Theorem 2. Consider the following form of equation (32):Then, the sufficient conditions for a stable motion are
Proof. The form of equation (60) is equivalent to the systemwhere,By taking Lyapunov function in the formwe obtainSo, every solution exists in the future and and as if the following conditions are satisfied:Apply the conditions; then, the conclusion holds.
Theorem 3. The governing equation, equation (27),has an asymptotically stable solution if
Proof. By considering thatTake that Lyapunov function isin the domain , , , , and . Hence, we haveThen, it is easy to obtain the general conditions:(i) and are locally Lipschitzian in and , respectively, and is continuous on (ii) and are periodic in of period (iii) for all and and is boundedApply the conditions; then, the conclusion holds.
4. Existence of Periodic Solutions
Let us assume that equation (32) has a periodic solution with period . Also, it is considered that, at the equilibrium point , the following assumption holds:
Hence, the only solution of equation (72) is . Equation (32) is rewritten as follows (, ):where the RHS represents the absorption and the supply of the system.
Thus, 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 if(i) and are locally Lipschitzian in and , respectively, and is continuous on (ii) and are periodic in of period (iii) for all and is bounded(iv)(v)There exist , , , and (vi)
Proof. For the domain , let , , and and . Consider that there exists a continuous function such thatAssume that the following identities hold:and then let the two functions and be as follows:Hence, it can be concluded that there exist in the interior of two domainsThen, the solutions and exist and are bounded for all . Then, the conclusion holds.
Theorem 5. In the autonomous case of the approximate governing equation , there exists a limit cycle under the following conditions:
Proof. By using the polar coordinatesone gets the following:After substitutions and using the conditions and , we obtainand one particular solution isCorresponding to the limit cycle,It is noticed thatThis means that all paths approach the shown radius of the limit cycle especially when . Then, the conclusion holds.
5. Averaging Analysis
We assume that the solution of the approximate form of governing equation, equation (35), is still given by
With time-varying coefficients and , in other words, we consider equation (85) as a transformation from to and . Hence, by imposing a third condition or a third equation independent of equation (35) and equation (85), assuming that and have the same form as the linear case,
and the following acceleration term is obtained:
Consider the following condition:
Then, we have to transform the system equation (85) to equation (88) with equation (35) into the two unknowns and as follows:where
Consequently, by considering that the major parts of and are slowly varying functions of time, they change very little during the time period and, in addition to the first approximation, they can be considered constant in the interval . Hence, we average both sides of the system equation (89) for and over the interval and obtain
Since and can be considered constant in the interval , and (and so and ) can be taken outside the integral signs in equation (91). Typically, the set of functions forms an orthogonal set over the interval . Then, we have, for all , from the inner product definition of two functions defined on the interval , that
Hence, equation (91) is reduced to
Consequently, the initial conditions can be given in terms of those for equation (35) by
5.1. Relative Equilibria and Stability
According to equation (93), it can be deduced that all equilibrium solutions are given by . In addition to the fact that the other paths corresponding to solutions of equation (35) are nonperiodic in general, let represent any of the three fixed points, , of equation (93). Hence, the linear analysis in the neighborhood of point is assumed to be
and then the corresponding linearized system iswhere
The general conditions of stability for every stable fixed point are
Moreover, the boundary between nodes and spirals can be given by the following identity:
The phase diagram computed for a particular case is shown in Figure 3 for and , respectively. In Figure 3(a), point is a center and represents a stable oscillation, but in Figure 3(b) there exist three equilibrium points: and are centers and is a saddle representing an unstable oscillation.
(a)
(b)
For the particular case of and = = = 1, the phase plane is shown in Figure 4 for and , respectively. In Figure 4(a), point is a spiral point and represents a stable oscillation, but in Figure 4(b) there exist three equilibrium points: and are stable spiral points and is a saddle representing an unstable oscillation.
(a)
(b)
The phase diagram computed for a particular case where , = −1, and = = 1 is shown in Figure 5 for and , respectively. In Figure 5(a), point is an unstable point as it bifurcates to generate a stable limit cycle, but in Figure 5(b) there exist three equilibrium points: and are centers and is a saddle representing an unstable oscillation.
(a)
(b)
5.2. The Jump Phenomenon
From equation (93), three are three different amplitudes at the equilibria. Consequently, due to the existence of the unstable node, the amplitude surface has unstable oscillations. Then, 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 is given bywhere
Hence, there may be as many as three positive values of satisfying the equation of amplitude. This 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 and unstable when .
Consider the following:
Since and 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 ; then,
Thus, the amplitude equation has three distinct real roots if is controlled by the overlap region.
To explain, we take the case of undamped motion ; then, we have
and then, as shown in Figure 6, the overlap region is controlled by
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). The 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.
(a)
(b)
(c)
As shown in Figure 7, it can be predicted that, for certain values of and at the lower values of , 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.
The consequences of this fold are shown illustratively using Figure 8 and 9. So, consider that 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 , 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. Then, here it must go down to the lower curve and then back to the upper one crossing the jumping point of the amplitude.
(a)
(b)
(c)
(a)
(b)
(c)
6. Homoclinic Bifurcation
The existence of homoclinic bifurcation with the support of the other problem parameters, particularly and , is considered. The 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 with an associated homoclinic trajectory with solution at . 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 and unstable if . As touches the critical value the path changes from heteroclinic spiral saddle connection to a homoclinic saddle connection at the critical value to heteroclinic saddle connection . This transition is due to the fact that parameter generates a homoclinic bifurcation at .
(a)
(b)
(c)
(d)
7. Route to Chaos
Considering equation (33), some special cases are looked at to allow the parameter space to vary where the chaotic phenomenon might occur. By using the numerical configuration, the critical value of at , , , , , and = = = 0 might be roughly obtained.
A simple numerical search reveals a stable periodic solution which is shown in Figure 11. Hence, if is increased further, the solution bifurcates at a specific value of into a solution which wraps around the origin twice. Further increase of leads to further bifurcation that takes place giving more periods of solution as vividly shown. The 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 increases.
(a)
(b)
(c)
(d)
This 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. The solutions’ figures as depicted show wandering solutions of an irregular oscillations type without any uniform pattern. Then, it is inferred that such solutions are said to display a chaotic behavior.
(a)
(b)
8. Conclusion
This 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: The 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. The effect of the fractional order and the system parameters on the stability of the solution are presented using Ince’s diagram. The 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. The existence of homoclinic bifurcation based on system parameters using Melnikov’s function is illustrated. The transition from homoclinic bifurcation to chaos is identified analytically and numerically using phase trajectories, Poincare’s map, and Lyapunov exponents.
Data Availability
The data used to support the findings of this study are available from the authors upon request.
Disclosure
This research was performed as part of the employment of the authors at University of Tanta and University of Kafrelsheikh, Egypt.
Conflicts of Interest
The authors declare no conflicts of interest.
Authors’ Contributions
All authors have read and agreed to the published version of the manuscript.