Applications of Jacobi’s Elliptic Functions in Forced Vibration in Nonlinear Rotor Dynamics

In this paper studied a rotary system with a nonlinearity, the equations of motion of which is one of the types of the Duffing ’s equations with multidegree-of-freedom, and also consider the advantages of using the elliptic function method to solve problems of this type. The system studied in this paper is also distinguished by the fact that in addition to the rotor vibrations in elastic supports with a nonlinearity, the vibrations of the foundation are also taken into account. The article presents a comparative analysis of the numerical Runge-Kutta-Fehlberg's 4-order method with an error estimate of 5-order, the approximate analytical Van der Pol ’s method and the elliptic function method proposed by the authors by comparing the obtained equations of motion of the system, as well as by comparing the constructed frequency response characteristic. From the results obtained, it follows that the method proposed by the authors can serve as a more accurate, general case of previously used approximate methods.


Introduction
Currently, technological equipment may consist of many high-speed, highprecision mechanisms, including a complex control system and rotating elements, i.e. rotor systems. Jeffcott [1], Stodola [2], Dimentberg [3], Kelson [4], Kramer [5], Tiwari. [6], Rao [7], Tondl [8], Adams [9], Yamamoto and Ishida [10], Muszynska [11], Genta [12] and others made contributions to the study of the dynamics of rotor systems. Their work formed the basis of those physical and mathematical representations that served the further development of this field. Focusing on various issues related to the peculiarities of rotor dynamics when taking into account various factors, we can distinguish the works of the following authors, which address individual issues, such as various methods for determining critical velocities, for example, using modal analysis [12]; construction of Campbell frequency diagrams for determining critical velocities taking into account gyroscopic moments [13,14]; conducting a complete analysis of rotor dynamics from determining critical speeds to obtaining the amplitudes of oscillations and orbits of the rotor motion under the influence of an imbalance in an ordinary and complex form [15]; studying the dynamics of rotors taking into account some features, for example, high-speed rotors with complex rotation [16]; ensuring the stability of movement of rigid and flexible rotors [17]; analysis of nonlinear oscillations of rotors caused by the presence of various factors, for example, the use of non-linear supports, which lead to the appearance of sub--and super--harmonic oscillations [18,19], self-oscillating modes [20], and other phenomena depending on the type of bearings used in the machine [21,22].
To date, there are numerous studies of the nonlinear dynamics of rotor systems [23][24][25][26][27][28][29][30]. It must be emphasized that despite the widespread occurrence of such systems, their dynamics have not been sufficiently studied due to the difficulties associated with the need to consider the combined action of factors, such as the influence of a liquid during partial filling, external non-conservative forces, linear and angular imbalances.
In the design of rotary machines, one of the most important mechanical components described by nonlinear models that determines the operability and reliability of the system [24][25][26][27][28] are elastic bearings, in our case, rolling bearings act as elastic bearings. Ignoring the nonlinear properties of bearings alters the results of both qualitatively and quantitatively [23,29,30]. This is explained by the fact that when analyzing linear rotor systems with rolling bearings, an approximate estimate of the stiffness and damping properties of bearings is most often used, when in reality the stiffness of a bearing substantially depends on loading, i.e. from the operating mode of the rotor system, from the geometry and size of the clearance in the bearing, from the size of the fit of the inner and outer rings in the bearing, etc. [31,32].
Bearing mathematical models that consider nonlinearity factors are distinguished by complexity and, first of all, by the loads that they consider. In many cases, the Hertz theory is used to describe the bearing model, which relates the radial loads acting on the bearing and the deformation at the points of contact between the rolling body and the bearing rings [11]. When describing the bearing model, it is assumed that there are no phenomena of slipping of bodies or rolling surfaces. Damping is usually incorporated in the formulation of equivalent viscous and linear friction [33]. Also, for the most complete description of the process, it is important to consider the influence a change in inertial parameters, imbalance and asymmetry of the rotor installation on the shaft, external friction, and various kinds of positional forces [33][34][35]. Such complications of the model when analyzing the dynamics make it possible to study the influence of the gap size, rotation frequency on the frequency spectra and amplitude-frequency characteristics for any rotor system on rolling bearings.
At the moment, the use of rolling bearings has many advantages compared to magnetic bearings and lubricated plain bearings. Rolling bearings are relatively simple to manufacture (the same plain bearings are more complicated to manufacture, moreover, the thickness of the oil film very much affects the dynamics of the rotor system. For example, when manufacturing it is worth considering that with an increase in the film thickness the amplitudes of self-oscillations also increase [6,10,36], and insufficient film thickness will lead to a violation of the thermal regime in the bearing, which leads to quick wear of the stud), they are more durable, and do not require special working conditions as magnetic bearings. In most cases, the approximation by a cubic power series is used to describe the nonlinear properties of rolling bearings [9,10,[38][39][40].
In this paper, we study the forced oscillations of a vertical rotor system, which is described by a system of inhomogeneous nonlinear (cubic) differential equations of the second order, which are known as the Duffing equations. The Duffing equation is one of the classical equations of nonlinear dynamics [41]. This equation describes a fairly wide range of phenomena of the world around us [42]. The The Duffing equation can conditionally be considered in certain cases as an equation with "weak" nonlinearity, and in other cases as an equation with "strong" nonlinearity. The case with "weak" nonlinearity allows approximate solutions of the equation, since as practice shows, the solution in this case is nearly periodic. In this regard, this case has been fairly well studied by such authors as Stoker [43], Rao [44], Nayfeh and Mook [47], Klotter [50], Timoshenko [49], Hayashi [45], Tondl [46], Magnus [51], Kauderer, Balachandran, Benaroya [48], Kelly. This case was solved with sufficient accuracy by various methods, namely by techniques of averaging, partial linearization, asymptotic, and graphical methods. For example, authors such as Stoker, Hayashi, Rao, on the assumption that oscillations occur relative to a stable state, solve this equation by substituting a periodic function [43][44][45].
Authors such as Tondl and Schmidt [46] use a different approach to solving this problem. They use an approximate method of slowly varying amplitudes, where the amplitude of the oscillations is considered as a function of time. To analyze the stability of stationary solutions, linearized equations in variations are used.
Other authors such as Nayfeh, Mook, Balachandran, Benaroya [47] use various asymptotic methods to solve the perturbed Duffing equation. For example, Nayfeh, Mook and Balachandran used the so-called multi-scale method, the idea of which is to switch from given variables to slowly changing variables. Benaroya used methods of expansion in a small parameter to solve this problem [48].
It should also be noted that authors such as Timoshenko and Klotter used the Ritz averaging methods [49,50], where it is assumed that an approximate solution of the steady-state oscillations can be represented as a series. Among the authors who obtained an approximate solution by the method of partial linearization, it should be noted Magnus who managed to obtain a relationship between the amplitude and frequency of the system by the harmonic balance method [51]. It should be noted that there are also exact solutions of the Duffing equation without damping, expressed through elliptic functions described by Cveticanin, Kovacic, Hsu [52,53,56] and Starossek [57]. Likewise, based upon the elliptic functions, some approximate methods were advanced, such as Poincare-Lindstedt and Krylov-Bogolyubov methods, generalized Galerkin method, and the multiple scales method.

Problem Statement
To construct the equations of motion of the system, a fixed coordinate system Оxyz is introduced (Figure 1). Let be that in the equilibrium state the geometric center of the shaft (rotor) and the center of gravity of the foundation coincides with the axis Oz. The coordinates in the shifted position of the center of the shaft (rotor) О1 are denoted by х1 and у1, and the coordinates of the center of gravity of the rotor are indicated by хs and уs. The coordinates of the center of gravity О2 of the foundation are indicated by х2 and у2, respectively. It is postulated that the rotor undertakes a plane-parallel motion, and there is no rotation of the foundation around the coordinate axes. To derive the equations of motion, we use the following notation. Let 0  denote angular velocity of the rotor (shaft); radial force and radial deformation in a ball bearing which are fitted respectively by using the Hertzian contact formulation and cubic polynomial, and the specific expressions are given as follows: the Hertzian contact, (1) and the cubic polynomial, where Fr represents the radial force (N), δr is the radial deformation (m), Cb is the contact stiffness (N/m 3/2 ), c0 and c1 are respectively, the cubic and linear coefficients of the cubic polynomial. It can be found that the Hertzian analytical formulation is in good agreement with the experimental data when the maximum deformation is under 140 µm, which is not large enough. However, the possibility that these two models are equivalent is in a wide range of deformation. Thus, the predicted maximum deformation value is 1000 µm based on the experimental Hertzian analytical formulation, and then the data are fitted by cubic polynomial [58]. In our case, we will call coefficients c0 and c1 are stiffness factors (or constants) of the rotor support (stiffness of the rolling bearing), c2stiffness coefficient of the foundation support; χ and χ0external friction coefficients acting on the rotor and foundation. The kinetic energy of the rotor is determined by the formula: where Jpolar moment of inertia of the rotor, mmass of the rotor. The kinetic energy of the foundation is expressed by the formula: where M is the mass of the foundation. The total kinetic energy of the system is given by The potential energy of the isotropic elastic field of a foundation support The potential energy of the isotropic nonlinear elastic field of the rolling bearing is determined by the formula We energy dissipation function reads The coordinates of the center of gravity of the rotor are given as follows: where еthe value of the rotor static imbalance (linear displacement of the rotor center of mass from its geometric center -О1Оs).
The total potential energy of the system is: . 12 П П П  (11) Now, using the Lagrange equations of the second kind, we write the equations of motion of the system: My c y c y y c y y y with the following initial conditions The system of equations (12) describes the coupled motion of an unbalanced rotor and foundation with non-linear supports. In the absence of foundation movement, system (12) will consist of two equations describing the equation of motion of the rotor in nonlinear supports, which have been studied in sufficient detail in the [54,55].
We introduce the following dimensionless parameters Then the system of equations (12) we obtain Consider the forced oscillations of the rotor and the foundation. Forced oscillations of the system are caused by rotor imbalance. With a harmonic disturbing force, after some transient process, periodic oscillations multiple of the frequency of the disturbing force (steady-state oscillations) are established. The first and third equations of system (15) are not related to the second and fourth equations and their solutions can be sought independently of each other. In the future we will use the solutions of the first and third equations. We write the first and third equations in the form:

Classic Approximate Solution
For the solve this problem and for the comparing the results that will be obtained using the elliptic function we use by classical averaged Van der Pol's method. To do this, we will find for the functions f1(τ) and f2(τ) in system (17) in the following where is A1(τ), A2(τ), B1(τ) and B2(τ) in slowly varying values for the period of the oscillations. Substituting (18) into the system of equations (17) and equating the left and right sides of the equations with the same functions, we obtain For stationarity of solutions to system (19), it is necessary that A1(τ), A2(τ), B1(τ) and B2(τ) be constant, then The amplitudes obtained from the system (20) (Fig. 2). Hence, geometrically, elliptic functions can be represented as where b = 1, and The function amp(u) is the inverse of the elliptic integral of the first kind, i.e.: The main properties of elliptic functions are

Application of Elliptic Functions for the Solution
To solve system (17), we assume that that is the perturbing force cos(ητ) there is a special case of the disturbing force of the form cn(ωfτ,k), that is, we assume that η ≡ ωf, and when k = 0, cn(ωfτ, 0) = cos(ητ). Given this statement, we rewrite the system (17) .
The solution of system (26) will be found in the form where kelliptic modulus, ωfcircular frequency of disturbing force, cn(ωf τ,k) and sn(ωf τ,k)elliptical cosine and sine.
Considering basic relationships of elliptic functions (25)  f where, changing the values ωf you can build the responses characteristics of the vibrations of the rotor and the foundation.
The E ≡ E(φ, k) the incomplete elliptic integral of the second kind, and K ≡ K(φ, k)

Results
To verify the results obtained by the elliptic function method, along with equations (33) and (34), were obtained and compared the results of expanding the elliptic cosine in a series (equation (36) (Figures 3 and 4).  were eliminated by the asymptotic expansion of variable and derivative of the Poincare's type, i.e. by the degree of parameter, characterizing the nonlinearity [47,42]. Among the many authors K.Magnus can be distinguished, who suggested the solution of this problem applying the method of harmonic balance, when the restoring force is linearized and the stiffness coefficient of that is a function of the amplitude [51]. Other authors, such as C.S.Hsu, L.Cveticanin, I.Kovacic, when searching for a solution to the nonlinear equation of motion, considered the more generalized method based on the application of Jacobi's elliptic functions, as a result of that the authors got the generalized system of equations, which degenerated with the ellipticity modulus equals to zero into the system of algebraic equations, that were got earlier by other authors [52,53,56]. In all of the above mentioned cases the maximum value of the amplitude, lying behind the field of the breakdown of the amplitudes, was found by superimposition the amplitude-frequency characteristics, obtained from the system of algebraic equations and excluding unstable stationary solutions. that possible work carried out in one cycle is zero, they suggested a graphical method to calculate this maximum value; this method allows us to determine the maximum value of the amplitude by intersection of the family of hyperbolas, corresponding to the fundamental resonance frequency (ωf = 1) with the curve of the amplitude-frequency characteristic of free vibration of the relevant system [49,50].
In our case the method of harmonic balance was applied in conjunction with the Jacobi's elliptic functions, that makes it possible to get the generalized system      it is increased by five times (case λ = 5) the first linear resonance occurs at ωf = 0.06 and the breakdown of the amplitudes is observed at ωf = 1.4 respectively, the values of the amplitudes will be f1 = 3.5, f2 = 3. In the case of the further increasing the ratio of the linear stiffness coefficients (case λ = 10) the first linear resonance is already observed at ωf = 0.07, and the breakdown of the amplitudes occurs at ωf = 1. 45. The values of the amplitudes of the rotor and the foundation will be f1 = 2.6, f2 = 2.

Conclusion
A method is developed for calculating the maximum value of the amplitudes and constructing the frequency characteristics of forced vibrations of the "rotorfoundation" system on rolling bearings with a nonlinear characteristic based on Jacobi's elliptic functions. 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, have been determined. The features of the linear and nonlinear behavior of a system with many degrees of freedom, caused by vibrations of the foundation, are shown.

Data availability statement
Data sharing not applicable to this article as no datasets were generated or analysed during the current study. Figure 1 Scheme of a rotor on rolling bearings Figure 2 Interpretation the motion of the center of mass of the rotor and foundation  Frequency response of the foundation using different methods     Frequency response of the rotor with different mass ratio Figure 10 Frequency response of the foundation with different mass ratio Figure 11 Frequency response of the rotor with different stiffness ratio Figure 12 Frequency response of the foundation with different stiffness ratio