Investigation of Nonlinear Forced Vibrations of the “Rotor-Movable Foundation” System on Rolling Bearings by the Jacobi Elliptic Functions Method

: The paper considers a rotor system with a nonlinear characteristic. Its equations of motion are a kind of Dufﬁng class equations with multiple degrees of freedom. The paper shows the advantage of using the method of elliptic functions for solving problems of this type. This method enables us to take into account not only vibrations of the rotor installed in elastic nonlinear supports, but also vibrations of the foundation. A comparative analysis of application of the method of elliptic functions proposed by the authors is carried out by comparing the derived equations of motion of the system, as well as by comparing the obtained amplitude-frequency characteristics with the results obtained by the numerical Runge–Kutta–Fehlberg’s 4-order method and the approximate analytical Van der Pol method. The regions of resonant frequencies for superharmonic oscillations and bifurcation regimes are determined. It is concluded that the method proposed by the authors is a more accurate and general case than the previously used approximate methods.


Introduction
A significant contribution to the study of linear and nonlinear dynamics of rotor systems was made by such scientists as Rao [1], Yamamoto and Ishida [2], Tondle [3], Mushinskaya [4], Jeffcott [5], Dimentberg [6], Stodola [7], Kelzon [8], Tiwari [9], Cramer [10], Adams [11], Vance [12], Penny [13], and others, whose works are the basis of those physical and mathematical concepts that contributed to its further development. One of the main components determining the reliability of rotary machines is elastic supports [14][15][16][17]. In this paper, rolling bearings with a non-linear characteristic act as elastic supports. The neglect of nonlinear properties in the study of the dynamics of vibrations of mechanical systems usually leads to overestimated amplitudes and frequencies, which, accordingly, worsens qualitative and quantitative results, and causes unnecessary additional costs in the design of rotor systems [18][19][20]. In the linear theory of the dynamics of rotary machines, damping and elastic properties of bearings are not accurately specified, which is especially important when it is necessary to study high-speed and high-frequency processes. In fact, the rigidity of a rolling bearing depends on many factors, among which one can note the size of the gap and its geometry [21], the type of load and the mode of operation of the system [22][23][24][25], the number of rolling elements in the bearing, the size and type of fit of the rings, etc. [24,[26][27][28]. Often, Hertz's contact theory is used to analyze and describe the restoring forces of a rolling bearing. It connects the deformation at the points of contact between the rolling element of the bearing with the loads acting on the bearing in the radial direction, without taking into account the slippage of the rolling elements or rolling surfaces [19,20,28]. To describe the nonlinear properties of the restoring force of rolling bearings, the cubic power series approximation [4,18,28,29] is used, since this approximation is in good agreement over a wide range of deformation.
In addition, for the most complete analysis of the process, it is important to take into account the influence of such factors as the asymmetry of the installation of the rotor on the shaft, i.e., imbalance, variability of inertial parameters, various positional forces, external friction, mobility of the foundation, etc. The impact of foundation vibrations on the linear dynamics of the "rotor-bearings" system and, in particular, on rotor vibrations was studied in detail in [5,6,10,11,30,31].
In Refs. [30,31], mathematical models of the "rotor-bearing-foundation" system, in a linear formulation, with the rotor shaft located horizontally, are presented. In a nonlinear formulation, this problem was partially considered and analyzed by the method of complex amplitudes in [32,33]. A more complicated mathematical model of this kind makes it possible to investigate the effect of rotational speed on frequency spectra and amplitudefrequency characteristics for any rotor system on rolling bearings.
In this paper, the motion of the rotor system is described by a system of cubically nonlinear non-homogeneous differential equations of the second order, known as the Duffing equations. In most cases, nonlinear differential equations are divided into equations with "weak" and "strong" nonlinearity according to the value of the nonlinearity parameter. The case with "weak" nonlinearity allows an approximate solution of the equation, since, as practice shows, the solution in this case is almost periodic. The case of eigenoscillations has been fairly well studied with sufficient accuracy by many authors using various methods, and the case of forced oscillations has been less studied [34].
It should be noted that there are also works on finding exact solutions to the Duffing equation for some isolated cases, expressed in terms of elliptic functions, described in the works of such authors as Tsvetitsanin, Kovacic, Hsu [35,36], and Starossek [37]. The solution of the Duffing equation using the method of elliptic functions was first shown by Hsu in [36]. On the basis of elliptic functions, the approximate methods of Poincare-Lindstedt, Krylov-Bogolyubov, the generalized Galerkin method, the multiple scales method, and the harmonic balance method were developed [35]. An analysis of the sources allows us to conclude that the method of elliptic functions has many advantages when searching for both approximate and exact solutions of quadratic oscillators, Duffing, and Helmholtz oscillators [34,35]. The method of Jacobi elliptic functions is also used in the study of rotor systems, in particular, in [38], where the Jeffcot rotor is considered. A review of the literature for the period of over 20 years shows that there are practically no theoretical works on the dynamics of vertical rotor systems in this formulation, where the method of investigation is the method of elliptic functions and the corresponding experimental works. In this regard, in this paper, for the first time, the system "Rotor-Rolling Bearings with a Nonlinear Characteristic-Movable Foundation" is investigated, which, as will be shown below, allows us to obtain more generalized results based on the method of elliptic functions.

Problem Statement
A symmetrical vertical rotor of mass m, with polar moment of inertia J, having a static imbalance e, rotates on rolling bearings with a constant angular velocity Ω 0 . The foundation of mass M is installed on elastic supports with an equivalent linear stiffness coefficient c 2 ( Figure 1). It is assumed that the rotor performs a plane-parallel motion, and there is no rotation around the coordinate axes. The radial compliance of bearings occurs due to the deformation of the rolling elements and raceways at the points of contact.
Two Cartesian coordinate systems are introduced, fixed-Oxyz and mobile, rigidly connected with the geometric center of the rotor-O 1 ξη, where η is the polar axis, and the ξ axis is arbitrarily drawn through the eccentricity vector of the rotor mass. The motion of the system is considered relative to the fixed coordinate system Oxyz. The coordinates of the geometric center of the rotor are designated as O 1 (x 1 , y 1 ), and its center of mass as O s (x, y). The foundation mass center coordinates are denoted as O 2 (x 2 , y 2 ), χ and χ 0 are damping coefficients (Figure 1). The non-linear restoring force of a bearing can be described using the Hertzian contact formula: where F C is the radial force, δ r is the radial deformation, C b is the non-linear stiffness coefficient.
For a qualitative analysis of the equations of motion of the system, the restoring force of rolling bearings in the form (1) can be approximated by a power series: where c 0 and c 1 are cubic and linear coefficients of a cubic polynomial, respectively. The kinetic energy of the system is defined as: Considering that: The potential energy of the isotropic non-linear elastic field of the system is defined as: The energy dissipation function reads: The coordinates of the center of mass of the rotor are given as: where e is the value of the linear displacement of the center of mass of the rotor from its geometric center O 1 O s , i.e., static imbalance. Using the Lagrange equations of the second kind, we write the equations of motion of the system in the following form: Equation (8) describes joint movement of the movable foundation and the rotor with static imbalance installed in elastic supports with a nonlinear characteristic. In the case of a fixed foundation, Equation (8) will consist of two Duffing equations describing the equation of motion of a nonlinear Jeffcot rotor, which have been studied in sufficient detail.
We introduce the following dimensionless parameters: Then, through the system of Equation (8), we obtain: ..

Equations of Motion's Solution
Since we are looking for a quasi-periodic solution of Equation (10) and, accordingly, we consider the case with weak nonlinearity, we will use the classic proven Van der Pol method to compare the results that will be obtained using the Jacobi elliptic function method. Since there are damping forces in the system, we represent the functions f 1 (τ) and f 2 (τ) in system (10) in the following form: where derivatives A 1 (τ), A 2 (τ), B 1 (τ), and B 2 (τ) are amplitudes slowly changing over the oscillation period. Substituting (11) into the Equation (10) and equating the left and right parts of the equations for the same functions and harmonics, we obtain: ..
The first and second order derivatives of the amplitudes of the rotor and the foundation are equal to zero, since for the stationarity of the solutions of system (12) it is necessary that A 1 (τ), A 2 (τ), B 1 (τ)), and B 2 (τ) be constant in time, then: From the algebraic equations of system (13), by varying the dimensionless frequency η, and analytically solving cubic amplitude polynomials for each case, the amplitudefrequency characteristics, shown in figures below, were obtained.
As the elliptic functions are periodic functions on the complex plane, they can be considered as a generalized case of trigonometric functions. In the case of stable forced oscillations, the oscillation trajectory of the center of inertia of a linear oscillator can be depicted as a circle (Figure 2, case (a)). In the case of small nonlinearity, the trajectory of the center of inertia of our system, which is a nonlinear oscillator, will describe an ellipse ( Figure 2, case (b)). Therefore, geometrically elliptic functions can be represented as: where b = 1, and: where here amp(u) is an inverse function of the elliptic integral of the first kind: Taking into account the above expressions (14)- (17), we represent in the form (18) the main relations for elliptic functions used in the following sections in the form: For a more accurate determination of the amplitudes and resonant frequencies of forced oscillations, we use the Jacobi elliptic functions. To solve Equation (10) by the elliptic function method, we assume that the perturbing force cos(ητ) is a special case of the elliptic cosine cn(ω f τ,k), i.e., we assume that η ≡ ω f , and for k = 0, cn(ω f τ, 0) = cos(ητ). Taking these statements into account, we rewrite system (10) in the following general form: ..
where ω f is a circular frequency of the perturbing force, k is the elliptic modulus, and cn(ω f τ,k) and sn(ω f τ,k) are the elliptic cosine and sine. The solution to Equation (19), as damping forces are present in the system, will be sought in the following form: Taking into account the main relations of elliptic functions, we obtain: Substituting the solution in the form (20) into system (19), and, taking into account (21), we obtain: where is: Here K ≡ K(ϕ, k) is a complete elliptic integral of the first kind, whereas E ≡ E(ϕ, k) is an incomplete elliptic integral of the second kind.
Equating the coefficients for the same functions and harmonics of Equation (22), we obtain a system of algebraic equations that enables us to determine the amplitudes of forced vibrations in the form (23): By varying the cyclic frequency in the system of Equation (23), the amplitude-frequency characteristics of the forced oscillations of the rotor and foundation were obtained (Figures 3-12). As can be seen, for k = 0, the system of Equation (23) degenerates into system (13). Equation (13) can be obtained by such classical approximate methods as the Van der Pol method, the harmonic balance method, etc. This fact is a direct confirmation of the generalization of the method of elliptic functions.

Results
The amplitude-frequency characteristics for an industrial centrifuge with parameters m = 24 kg, M = 25 kg, e = 0.001 m, c 0 = 1.1 × 10 7 kg/s 2 , c 1 = 0.87 × 10 11 kg/m 2 s 2 , c 2 = 3.26 × 10 5 kg/s 2 , χ 0 = 6.59 kg/s, χ = 4200 kg/s were obtained by the following methods: from the algebraic system of Equation (13) using the Van der Pol's method (Figures 3 and 4, the blue curve), from differential Equation (8) using the 4th order Runge-Kutta-Fehlberg method (Figures 3 and 4, the dotted red curve), from differential Equation (8) at c 1 = 0 by the 4th order Runge-Kutta-Fehlberg method (Figures 3 and 4, the black curve) and compared with the amplitude-frequency characteristic obtained from the algebraic system of Equation (23) derived by the method of Jacobi elliptic functions (Figures 3 and 4, the green curve). As a result of calculations, three resonant regions were determined for the rotor and the foundation, the presence of which, as will be shown below, is due to the nonlinearity of the supports and the mobility of the foundation.
The curve of the amplitude-frequency characteristic in the first resonant region is typical of linear oscillations and is observed for all cases at ω f = 0.051. When calculating the amplitudes by the Van der Pol method (Figures 3 and 4, the blue curve), the deviation from the results of the numerical method (Figures 3 and 4, the dotted red curve) in the case of the rotor was 2.8%, and in the case of the foundation it was 9.1%. In cases of calculation by the method of elliptic functions, the amplitudes for the first resonance region are the closest to the results of the numerical method. The deviation of the amplitude values for the rotor was 1.14%, and for the foundation 6.4% (Figures 3 and 4, the green curve). The amplitude values are presented in Table 1.   The amplitude-frequency characteristic curve for the second resonant region is typical of nonlinear oscillations with a rigid characteristic. Nonlinear oscillations and bifurcation modes in this case are due to nonlinearity of the restoring forces of the elastic supports. In this case, the amplitude breakdown for each method occurs at different frequencies with different amplitudes. For example, for the 4th order Runge-Kutta-Fehlberg numerical method, the resonance was observed at ω f = 1.33, in the case of the Van der Pol method, the deviation from the amplitude of the numerical method for the rotor was 17.6%, for the foundation 33.87%, while the resonance occurred at ω f = 1.25. In the case of Jacobi elliptic functions, this breakdown can be observed at ω f = 1.3, the deviation from the amplitudes for the rotor is 8.66%, and for the foundation it is 31.7%. As can be seen in the latter case, the deviations of the values of the resonant frequency and amplitude from the numerical method are much smaller. Amplitude values are presented in Table 2. Nonlinear systems are characterized by anharmonicity, i.e., the presence of ultraharmonic and subharmonic resonances that are multiples of the frequency of the main resonance. For the linear case, the main resonance occurs at ω f = 1. In our case, the main resonance region, due to nonlinearity, is observed at 1 ≤ ω f ≤ 1.33. In order to determine additional resonances, the amplitude-frequency characteristics of the corresponding linear case (Figures 3 and 4, the black curve), the amplitude-frequency characteristics obtained by the numerical method (Figures 3 and 4, the red curve), and the amplitude-frequency characteristics of the analytical by the Van der Pol method (Figures 3 and 4, the blue curve) were constructed. Since the second resonance, located to the right of the main one at ω f = 1, is not observed in the linear case and is a multiple of the main one, it was found that it is an ultraharmonic resonance.
It should be noted that in all the above cases, the main disadvantage of using the numerical method for the nonlinear case is the impossibility of determining the amplitudes of the rotor and the foundation for the amplitudes lying behind the break point. In the case of using the Van der Pol method, large deviations are observed both for the amplitudes and for the resonant frequencies. Jacobi's method of elliptic functions is free from these shortcomings due to its generality.
To optimize the workflow using the method of elliptic functions, a parametric analysis was carried out by constructing the amplitude-frequency characteristics of the rotor and foundation for different values of the nonlinearity parameter ε (Figures 5 and 6), for different values of the damping coefficients ζ 1 and ζ 2 (Figures 7 and 8), for different values of the linear stiffness parameter λ (Figures 9 and 10), and for different values of the mass ratio of the foundation and the rotor µ (Figures 11 and 12). All variable parameters are reduced to dimensionless form according to (9).
An increase in the nonlinearity parameter up to one order of magnitude inclusive ( Figures 5 and 6: cases ε = 5 and ε = 10) leads to a decrease in the amplitudes of the rotor and foundation in the region of the first resonance. For example, for ε = 5, the rotor amplitude decreases by a factor of 1.26, while the foundation amplitude decreases by a factor of 1.08. Further, with an increase in the nonlinearity parameter, i.e., for ε = 10, the decrease in amplitudes is also more intense. For the rotor, in this case, the amplitude relative to ε = 1 already decreases by a factor of 2.27, and for the foundation, by a factor of 1.4. With a decrease in the nonlinearity parameter within one order of magnitude, the opposite picture is observed (Figures 5 and 6: cases ε = 0.1 and ε = 0.5). For example, for ε = 0.5, the rotor amplitude increases by 1.5 times, and the foundation amplitude by 1.13 times. In the case of a decrease in the nonlinearity parameter by one order of magnitude, i.e., when ε = 0.1, the increase in amplitudes is more intense and exceeds the initial values by a factor of 1.96 for the rotor and by a factor of 1.43 for the foundation. No shifts in the region of the first resonant frequencies are observed in these cases. The absence of a shift along the frequency axis of the first resonance is due to the linear nature of the oscillations in this region. In addition, as shown earlier, these resonant amplitudes arise solely due to mobility and oscillations of the foundation. In the region of nonlinear resonant frequencies, the magnitudes of the amplitudes do not depend on the variation of the nonlinearity parameter. With an increase in nonlinearity, the slope of the amplitude-frequency characteristic curve increases, and the amplitude breakdown point shifts in the direction of higher frequencies. For example, for ε = 5, the drop in amplitudes shifts to ω f = 1.4, whereas for ε = 10, a drop is already observed at ω f = 1.43. As the nonlinearity decreases, the slope of the curve correspondingly decreases, degenerating into a linear case. The frequency and amplitude values are given in more detail in Table 3 and in Figures 5 and 6. Variation of the nonlinearity parameter only affects the slope at the second resonance and the amplitude values, since the non-linearity is known to limit the magnitude of the amplitudes.   In the case of an increase in the damping coefficients of the rotor and the foundation (Figures 7 and 8: cases ζ 1 = 5 and ζ 2 = 5, ζ 1 = 10 and ζ 2 = 10), the amplitudes for both the first and the second resonant regions decrease. For example, at ζ 1 = 5 and ζ 2 = 5, the foundation and rotor amplitudes at the first and second resonances decrease by a factor of 1.5. In the case of an increase in damping by one order of magnitude, i.e., ζ 1 = 10 and ζ 2 = 10, the amplitudes of the foundation and the rotor at the first and second resonances already decrease by a factor of 2.15. Shifts in the region of both the first resonant frequencies and the second ones are not observed in these cases. As the damping decreases ( Figures 7 and 8: cases ζ 1 = 0.1 and ζ 2 = 0.1, ζ 1 = 0.5 and ζ 2 = 0.5), the amplitudes of the rotor and foundation increase correspondingly in the first and second resonant regions. For example, when the damping coefficient of the rotor and foundation is halved, i.e., at ζ 1 = 0.5 and ζ 2 = 0.5, an increase in their amplitudes by a factor of 1.5 is observed. In the case of damping decrease by one order of magnitude, i.e., at ζ 1 = 0.1 and ζ 2 = 0.1, the foundation and rotor amplitudes already increase by a factor of 2.25. No shifts of the resonant regions are observed in these cases either. It is seen that an increase in the damping coefficient of the supports generally has a good effect on the dynamics of the system, however, due to the presence of structural limitations and the risk of shock loads at large values of the coefficient, the variation in the damping coefficients should be limited. The frequency and amplitude values are presented in more detail in Figures 7 and 8, and in Table 4.   With an increase in the mass of the rotor (Figures 9 and 10: cases µ = 10 and µ = 5), an increase in the amplitudes of the rotor and foundation is observed in all areas of resonant frequencies. In the case of a fivefold increase in the mass of the rotor, i.e., at µ = 5, the rotor amplitude increases by a factor of 1.4, and the foundation amplitude by a factor of 1.1. With a further increase in the mass ratio, the amplitudes also increase, for example, at µ = 10, the difference in the amplitude ratios with the case of µ = 1 already reaches 2.25 times. With an increase in the mass of the foundation (Figures 9 and 10: cases µ = 0.5 and µ = 0.1), in all cases and in all areas of resonant frequencies, a decrease in the amplitude values is observed, since in these cases the foundation acts as an anti-weight and dampens the oscillations of the system as a whole. For example, the values of the amplitudes at µ = 0.5 for the rotor decrease in magnitude by 1.06 times, and for the foundation by 1.27 times. With a further increase in the mass of the foundation, the decrease in the amplitudes is more intense and is already less than for the case µ = 1 by 2 and 1.97 times. The values of frequencies and amplitudes for different values of the mass ratio parameter are presented in more detail in Table 5 and in Figures 9 and 10.   When varying the linear stiffness coefficient λ, both a change in the magnitude of the amplitudes and their shift along the frequency axis are observed. With an increase in the coefficient c 0 (Figures 11 and 12: cases λ = 0.1 and λ = 0.5), the first and second resonant frequency regions shift to the left along the frequency axis, towards its decrease. In the case when λ = 0.5, the first region of resonant frequencies shifts to the value of ω f = 0.04, while the drop in amplitudes characteristic of the second region of resonant frequencies is already observed at ω f = 1.25. With a further decrease in the linear stiffness parameter, i.e., at λ = 0.1, the first resonant region already sets in at ω f = 0.035, and the drop in amplitudes is already observed at ω f = 1.15. As the stiffness coefficient c 0 decreases, the reverse picture is observed. For example, at λ = 5, the first region of resonant frequencies shifts in the direction of higher frequencies to the value of ω f = 0.06, and the amplitudes drop already at ω f = 1.4. In the case of a further decrease in this coefficient, the first resonance occurs at ω f = 0.07, and the breakdown of the amplitudes occurs at ω f = 1.45. More detailed data on the values of frequencies and amplitudes are presented in Figures 11 and 12.  For high-speed rotor systems, when designing a rigid rotor rotating in elastic bearings, in order to reduce the amplitudes and vibration overloads, its first and second critical speeds are usually below the operating speeds, whereas in the operating frequency zone the selfcentering effect is used. In our case, as the range of operating speeds is proposed to take the frequency interval between the first and second resonance, the region between the second nonlinear and ultraharmonic resonance, or the interval lying behind the ultraharmonic resonance. In these areas, due to the self-centering effect, there is a significant decrease in the magnitudes of the amplitudes.

Conclusions
Thus, for optimal operation, it is necessary that in the stationary mode the operating speed of the system be in the interval between the first linear and the second nonlinear resonance, i.e., 0.051 < ω f < 1.301 (the first working region), or between the second and ultraharmonic resonances, i.e., 1.301 < ω f < 3 (the second working region), or in the region behind the ultraharmonic resonance (the third working region), where, due to the selfcentering effect, the values of the rotor and foundation amplitudes decrease.
An increase in the nonlinearity parameter up to ε = 10 has a positive effect on the amplitude values. In this case, it should be taken into account that the frequency of the amplitude breakdown is shifted towards higher rotor speeds.
The operation of the system is also positively affected by an increase in the damping coefficients (ζ 1 , ζ 2 ) and an increase in the mass of the foundation, as in this case the foundation acts as an anti-load and contributes to vibration damping.
An increase in the linear rigidity parameter λ leads to a decrease in the first resonance and a shift in the second resonance, which also contributes to a smoother transition through the resonant frequencies.
In all of the above cases, it should be taken into account that the limit of the maximum values of parameters of nonlinearity, damping, foundation mass and linear rigidity is limited by the design and technical features of a particular industrial centrifuge.
A technique that enables us to calculate the maximum value of the amplitudes and construct the frequency characteristics of forced oscillations of the "rotor-foundation" system on rolling bearings with a nonlinear characteristic based on Jacobi elliptic functions has been developed. The optimal parameters associated with the mass of the foundation, the coefficients of linear and nonlinear stiffness, as well as the damping coefficients, at which the magnitudes of the amplitudes have optimal values without the need for fundamental structural changes, are determined. The specific features of linear and nonlinear behavior of the system with many degrees of freedom caused by foundation vibrations are considered.