Nonlinear Dynamics and Chaos of a Nonsmooth Rotor-Stator System

Rotor systems have wide applications in industries, including aero engines, turbo generators, and gas turbines. Critical behaviors promoted by the system unbalance and the contact between rotor and stator lead to important nonlinearities on system dynamics. This paper investigates the complex behavior presented by a rotor-stator system’s dynamics due to intermittent contact. A four-degree-of-freedom Jeffcott nonsmooth rotor/stator system is used to describe the rotor behavior, while a viscoelastic suspended rigid cylinder represents the stator. Numerical simulations are carried out showing rich dynamics that include periodic, quasiperiodic, and chaotic responses. Special attention is dedicated to chaotic behavior and the calculation of Lyapunov exponents is employed as a diagnostic tool.


Introduction
Rotor systems have wide applications in industries, including aero engines, turbo generators, and gas turbines.The unbalance of the system, represented by an eccentric mass of the rotor, produces forced vibrations that can become critical with the deterioration through use.Another critical situation for rotor systems is the contact between rotor and stator that lead to important nonlinearities on system dynamics.Vibration reduction is one of the main goals in rotating machinery researches.The mathematical modeling of this kind of system is complex due to several reasons.In general, nonsmooth nonlinearities are preponderant being associated with discontinuous characteristics as intermittent contacts or dry friction.
Concerning the rotor-stator behavior, some works are focused on the occurrence of full-rub, a phenomenon where the contact between rotor and stator is permanent and is maintained by the friction [14][15][16].In general, the vast majority of the works considered two-degree-of-freedom (2DOF) models, related to rotor displacement, where the rotor is modeled as a Jeffcott rotor that may have contact with a stator flexibly mounted [10,11,13,[17][18][19][20][21][22].Considering 2DOF rotor-stator systems, Varney and Green [23] analyzed a 2DOF system considering a rotor with asymmetric support stiffness.Muszynska and Goldman [24] and Chu and Zhang [25] presented dynamical analysis focusing on the presence of nonperiodic solutions, including chaos.Still considering 2DOF, some authors focus on the contact modeling [26][27][28].
2DOF models of rotor-stator systems were deeply investigated in the literature; in contrast, few works address the analysis of higher DOF models, as considered by Popprath and Ecker [29], Cole [26], Patel et al. [30], and Childs and Bhattacharya [31].

Mathematical Problems in Engineering
From all cited literature, the mathematical models can be divided into two categories.The first one considers a rotorstator connection represented by stiffness, and this connection includes noncontact situation.The stiffness varies from noncontact to contact conditions.This type of model can represent rotor-bearing-stator systems.The second category of models is related to independent rotor-stator behavior in noncontact mode.Therefore, contact establishes the coupling between rotor and stator.This type of model can represent rotor systems such as aero engines, gas turbine, and turbo generator.
This paper deals with the dynamical analysis of a rotorstator system employing a four-degree-of-freedom (4DOF) nonsmooth system where both the rotor and the stator dynamics are taken into account.It is assumed that there is not rotor-stator coupling in noncontact situation.A Jeffcott rotor model is used to describe the rotor system behavior, while a viscoelastic suspended rigid cylinder represents the stator.The model allows one to study the rotor motion when intermittent contact occurs, which outcomes in complex dynamical behavior.Dry friction is considered during the contact using the same approach proposed by Popprath and Ecker [29].The main goal of this work is to investigate the details of system dynamics using nonlinear tools to characterize system behavior from qualitative and quantitatively points of view.In this regard, it is important to highlight the use of Lyapunov exponents to diagnose periodic, quasiperiodic, and chaotic responses.

Mathematical Model
The rotordynamic system investigated in this paper is represented by the horizontal rotor-stator system presented by the schematic picture of Figure 1.It consists of a rigid rotor disk with mass   and 2DOF (  ,   ), mounted on a massless elastic shaft, with stiffness   , and subjected to linear viscous damping   ; a viscoelastic suspended stator, with mass   and stiffness   , and subjected to linear viscous damping   , with two translational DOF (  ,   ), is also considered.The 2DOF considered for the stator distinguishes this work from the majority in the literature that considers rotor-stator systems.Under these assumptions, the rotor-stator system is described by a 4DOF model where there is a gap  separating rotor and stator.
This Jeffcott rotor consists of a cylindrical rigid structure rotating with an angular velocity, Ω, and with eccentricity, , which generates unbalanced forces.The stator is represented by a rigid cylindrical structure surrounding the rotor.In this regard, equations of motion are written in state space as follows: or, in other words, where {B} = {B  } + {B  } is the force vector composed by an unbalanced component, {B  }, which is a result of the unbalanced masses on the rotor, and the contact component, {B  }, that takes into account contact forces.The next section discusses the definition of this forcing vector that is directly related to contact model.Parameters presented in (1) are defined as follows: 2.1.Forcing Vectors.The rotor system dynamics have two possible modes: contact and noncontact.Figure 2 presents a schematic picture of the contact mode, showing the contact between the rotor and stator.In the illustration, the forces represented by   and   are, respectively, the normal and tangential forces applied by the rotor on the stator in the case of a clockwise rotation.These forces are related to each other by the friction coefficient  in such a way that In order to establish the contact condition, an auxiliary variable, , is defined as follows:  Based on that, the forcing vector {B  } is defined representing nonsmooth characteristics with contact and noncontact situations.When  ≤ 0, no contact is verified and {B  } = {0}.Otherwise, when  > 0 the contact force component is given by where   represents the contact stiffness.The use of large values of this contact coefficient represents the perfect rigid bodies' assumption.For details, please refer to impact theory presented in Popprath and Ecker [29].The contact description is complex, being related to several nonlinear phenomena, as treated in detail by Bhushan [32].In the literature related to rotor-stator system, however, it is usual to consider a simplification where the normal component of the contact force has a linear dependence with the indentation (radial intrusion depth), as considered in the present work.
Besides that, unbalanced forces are described by the following vector: where   stands for the gravity acceleration and  0 is the initial phase shift.Therefore, the forcing vector {B} is given by

Numerical Procedure
Nonsmooth nonlinearities introduce several difficulties on the numerical procedure employed to solve equations of motion.Basically, the state space is split into subspaces related to contact and noncontact modes and the transitions between them.Therefore, a numerical approach needs to be developed in order to properly deal with these characteristics.
A variable time step fourth-order Runge-Kutta method is employed here.This approach allows one to reduce the computational effort using smaller time steps during contact mode and transitions.This strategy is similar to the one presented by Conte and de Boor [33], which controls the time step by monitoring the estimated error of the fourth-order Runge-Kutta calculation.The estimated error is obtained by calculating the state variables X +1 through two different ways: using a single time evolution, being Δ the time step; and using two consecutive time evolutions, being Δ/2 the time step.As shown by Conte and de Boor [33], the error can be estimated by where X Δ/2 ( +1 ) represents the state variables evaluated using two time steps of Δ/2 each, X Δ ( +1 ) represents the one using one single time step of Δ, and  is the Runge-Kutta method order; thus,  = 4.
Note that, during transition between contact and noncontact modes, X +1 verifies the contact condition but its precedent X  does not.The time step calculation acknowledges the contact forces on the second step, while the previous step does not.This difference leads to an increase on the estimated error and, therefore, demands a decrease on the time step.When the system is inside or outside the contact subspace, the time step is kept constant.Nevertheless, when a transition is verified, namely, when the system passes from contact to noncontact modes and vice versa, the variable time step strategy is effective.When a transition from noncontact to contact is verified, the time step is decreased until the error is equal or smaller than a desired tolerance  max .Otherwise, when the transition occurs, from contact state to noncontact state, the time step that was decreased in the beginning of contact returns to a standard value, fixed before the beginning of the simulation.
When the error is greater than the desired value  max , a smaller time step Δ  can be determined by the ratio between the actual error   and the desired error, as shown in This variable time step approach is actually an iterative process that is repeated until the condition   ≤  max is satisfied.The  exponent is a constant that determines the number of iterations that needs to be carried out until the error reaches the desired value.This value is adjusted in order to obtain the minimum calculation time, which corresponds to making the least possible number of iterations without decreasing the time step excessively.In the majority of the cases presented in this paper,  ≈ 1 presents good results, and usually 1 to 3 iterations are necessary to reach the desired error.
The diagnostic of chaotic behavior is evaluated estimating Lyapunov exponents.Algorithms that computes Lyapunov spectrum of smooth systems are well established in the literature, as the one proposed by Wolf et al. [34].This algorithm monitors the local divergence of nearby orbits considering a reference and a perturbed orbit.Essentially, the algorithm uses a linearized version of the equations of motion to treat perturbed orbit.Due to system discontinuities, special attention should be given since nonsmooth equations of motion are represented to a switch model between contact and noncontact modes.Thus, the procedure described by Müller [35] for Lyapunov exponents calculation of discontinuous systems is employed in this work.The author proposes that linearized equations have to be supplemented by transition conditions in order to deal with mode changes.

Numerical Simulations
Numerical simulations are carried out in order to investigate the dynamical behavior of the rotor-stator system.Table 1 presents parameters employed in simulations, the same used by Demailly [36].
The analysis is focused on situations where the rotation speed is close to rotor natural frequency, 140 rad/s.Under this condition, the system presents high values of rotor displacement, representing a critical operational condition.Bifurcation diagrams, phase spaces, and Poincaré sections are used to present results.
Initially, a global investigation of the system behavior near the critical operational condition is carried out.In this regard, Figures 3 and 4 present bifurcation diagrams related to stator and rotor velocity under the slow quasistatic variation of frequency.Figure 3 shows the bifurcation diagrams for stator's horizontal velocity by increasing the rotational speed, Ω.Some enlargements are presented in order to highlight the system behavior details. Figure 4 shows the bifurcation diagram for the rotor's horizontal speed.These diagrams show the complex behavior of the system, presenting bifurcations, and several kinds of responses including periods 1, 2, and 4 as well as nonperiodic behavior.Figure 4 shows, together with the bifurcation diagrams, all Lyapunov spectrum for 7 different values of forcing frequency.Concerning Lyapunov exponents, periodic conditions, where all exponents are negative, quasiperiodic conditions, where the greater exponent vanishes, and chaotic conditions, where the greater exponent is positive, should be pointed out.Multistability of solutions is also possible but they are not identified from these results.
Based on the bifurcation diagram analysis, different system behaviors are investigated by assuming distinct values of rotational frequency.Phase spaces and Poincaré sections projections are used to present results.Basically, two subspaces are considered: one associated with stator (  ,   ) and the other one associated with rotor (  ,   ). Figure 5 shows a period-1 behavior related to Ω = 131.7 rad/s, while Figure 6 shows a period-2 behavior related to Ω = 132.2rad/s.It is noticeable that the stator orbit is more  complex compared with the rotor motion.Moreover, it is important to observe that stator orbit significantly changes after the period doubling bifurcation shown in Figure 3, while the rotor orbit keeps the same elliptic shape.Even in cases of higher periodicity, such as period-4 or chaotic behaviors, the rotor orbits preserves its elliptic shape.As in the case shown in Figure 6(b), the two points of Poincaré section seem to be only one.The enlargement of phase space, however, shows two different points.From Figures 5 and 6 it is possible to observe the differences of phase spaces, where the stator presents more loop when compared to the rotor.This behavior can be explained by the fact that there is no influence between rotor and stator in noncontact condition.
In contact condition, for example, the forcing frequency is predominant in steady state response of both bodies, while in the noncontact situation the rotor response is mainly in   the forcing frequency while the stator responds in its natural frequency.
The coexistence of attractors is an important characteristic that commonly appears in nonlinear systems.This situation is critical since the system response can jump from one attractor to the other due to small perturbations.The rotor-stator system presents this kind of behavior for some forcing frequencies, for example, when Ω = 133 rad/s.By increasing even more the rotational frequency to Ω = 140 rad/s, the system presents a chaotic response shown in Figure 11 in stator phase space together with its Poincaré section.Figure 12

Conclusions
This paper deals with the analysis of a nonsmooth rotorstator system.The Jeffcott rotor model is used to describe the rotor behavior, while a viscoelastic suspended rigid cylinder represents the stator.System dynamics are described by a 4DOF nonsmooth system with dry friction forces representing contact.Rotor and stator behaviors are evaluated       varying the rotation frequency.The range of analyzed forcing frequency is near the linearized rotor natural frequency, which represents a critical operational regime.Results show that system presents a complex behavior in very narrow frequency range.Several bifurcations are observed showing that system response can abruptly change with very slight change in forcing frequency.Different responses are observed including periodic responses (period-1, period-2, period-3, and period-4) and quasiperiodic and chaotic behaviors.
Besides the richness of responses, the coexistence of periodic attractors is observed.It is important to identify this kind of situation since system behavior can jump from one response to the other with tiny perturbation, which can be critical to the system.This rich behavior shows the complex dynamics associated with rotor-stator system dynamics pointing to the necessity of a proper nonlinear dynamical analysis for system design.

Figure 1 :
Figure 1: Schematic picture of the rotor-stator model.

6 Figure 4 :
Figure 4: Bifurcation diagram for ẋ by increasing Ω and Lyapunov exponents for selected forcing frequencies.

Figure 7
Figure7shows the basin of attraction for this frequency by modifying rotor initial position and velocity.All other initial conditions kept unchanged.The black points are related to a period-4 steady state, while red points correspond to a period-5 orbit.These orbits' phase spaces and Poincaré sections are shown in Figure8.The periodicity of the previous results, presented in Figures5, 6, and 8, is confirmed by the Lyapunov exponents, whose values are all negative.More complex responses are now in focus considering other kinds of behaviors.Initially, the rotational frequency is increased to Ω = 135 rad/s.Figure9presents phase space projection together with Poincaré section.The Poincaré section is highlighted in Figure 10.Lyapunov exponents estimation furnishes the following spectrum, where the maximum null value indicates a quasiperiodic response:  = {0.0;−2.8; −6.7; −9.7; −25.8; −35.4; −36.9; −54.1}.