# Nonlinear Dynamics and Chaos of a Nonsmooth Rotor-Stator System

Authors: André Brandão; Aline Souza de Paula; Marcelo Amorim Savi; Fabrice Thouverez
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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.The literature presents several works dealing with the dynamics of nonsmooth systems. Usually, this kind of system presents an unusual complex behavior and their description involves many mathematical and numerical difficulties [1–3]. Savi et al. [4] discussed a one-degree-of-freedom (1DOF) system considering numerical and experimental approaches. Several engineering systems are related to nonsmooth nonlinearities including oil drilling [5–8], manufacturing [9], and rotor dynamics [10–13].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–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–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–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].From all cited literature, the mathematical models can be divided into two categories. The first one considers a rotor-stator 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 rotor-stator 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.


## 2. Mathematical Model

The rotordynamic system investigated in this paper is represented by the horizontal rotor-stator system presented by the schematic picture of Figure1. It consists of a rigid rotor disk with mass mr and 2DOF (xr,yr), mounted on a massless elastic shaft, with stiffness kr, and subjected to linear viscous damping cr; a viscoelastic suspended stator, with mass ms and stiffness ks, and subjected to linear viscous damping cs, with two translational DOF (xs,ys), 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.Figure 1
Schematic picture of the rotor-stator model.This Jeffcott rotor consists of a cylindrical rigid structure rotating with an angular velocity,Ω, and with eccentricity, e, 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: (1)x˙rx¨ry˙ry¨rx˙sx¨sy˙sy¨s=01000000-ωr-2ξrωr0000000001000000-ωr-2ξrωr0000000001000000-ωs-2ξsωs0000000001000000-ωs-2ξsωsxrx˙ryry˙rxsx˙sysy˙s+Bor, in other words, (2)X˙=AX+B,where B={Bd}+{Bc} is the force vector composed by an unbalanced component, {Bd}, which is a result of the unbalanced masses on the rotor, and the contact component, {Bc}, 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:(3)ωs=ksms;ωr=krmr;ξs=cs2ksms;ξr=cr2krmr.


### 2.1. Forcing Vectors

The rotor system dynamics have two possible modes: contact and noncontact. Figure2 presents a schematic picture of the contact mode, showing the contact between the rotor and stator. In the illustration, the forces represented by fn and ft 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 ft=μfn.Figure 2
Contact between rotor and stator.In order to establish the contact condition, an auxiliary variable,g, is defined as follows: (4)g=xr-xs2+yr-ys2-δ.Based on that, the forcing vector{Bc} is defined representing nonsmooth characteristics with contact and noncontact situations. When g≤0, no contact is verified and Bc=0. Otherwise, when g>0 the contact force component is given by (5)Bc=gkcg+δ0xs-xr0ys-yr0xr-xs0yr-ys+μ0yr-ys0xs-xr0ys-yr0xr-xs,where kc 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:(6)Bd=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag,where ag stands for the gravity acceleration and φ0 is the initial phase shift.Therefore, the forcing vector{B} is given by (7)B=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag,g≤0,B=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag+gkcg+δ0xs-xr0ys-yr0xr-xs0yr-ys+μ0yr-ys0xs-xr0ys-yr0xr-xs,g>0.


## 2.1. Forcing Vectors

The rotor system dynamics have two possible modes: contact and noncontact. Figure2 presents a schematic picture of the contact mode, showing the contact between the rotor and stator. In the illustration, the forces represented by fn and ft 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 ft=μfn.Figure 2
Contact between rotor and stator.In order to establish the contact condition, an auxiliary variable,g, is defined as follows: (4)g=xr-xs2+yr-ys2-δ.Based on that, the forcing vector{Bc} is defined representing nonsmooth characteristics with contact and noncontact situations. When g≤0, no contact is verified and Bc=0. Otherwise, when g>0 the contact force component is given by (5)Bc=gkcg+δ0xs-xr0ys-yr0xr-xs0yr-ys+μ0yr-ys0xs-xr0ys-yr0xr-xs,where kc 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:(6)Bd=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag,where ag stands for the gravity acceleration and φ0 is the initial phase shift.Therefore, the forcing vector{B} is given by (7)B=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag,g≤0,B=0eΩ2cos⁡Ωt-φ00eΩ2sin⁡Ωt-φ0-ag000-ag+gkcg+δ0xs-xr0ys-yr0xr-xs0yr-ys+μ0yr-ys0xs-xr0ys-yr0xr-xs,g>0.


## 3. 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 Xi+1 through two different ways: using a single time evolution, being Δt the time step; and using two consecutive time evolutions, being Δt/2 the time step. As shown by Conte and de Boor [33], the error can be estimated by(8)Ei+1=XΔt/2ti+1-XΔtti+12p-1,where XΔt/2ti+1 represents the state variables evaluated using two time steps of Δt/2 each, XΔtti+1 represents the one using one single time step of Δt, and p is the Runge-Kutta method order; thus, p=4.Note that, during transition between contact and noncontact modes,Xi+1 verifies the contact condition but its precedent Xi 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 Emax. 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 valueEmax, a smaller time step Δtr can be determined by the ratio between the actual error Ei and the desired error, as shown in (9)Δtr=ΔtEmaxEiα.This variable time step approach is actually an iterative process that is repeated until the conditionEi≤Emax 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.


## 4. Numerical Simulations

Numerical simulations are carried out in order to investigate the dynamical behavior of the rotor-stator system. Table1 presents parameters employed in simulations, the same used by Demailly [36].Table 1
System parameters.
ω
s
r
a
d
/
s
ω
r
r
a
d
/
s
ξ
s
ξ
r
δ
m
μ
m
r
k
g
m
s
k
g
e
m
k
c
N
/
m
E
m
a
x
m
540
140
0.05
0.02
0.006
0
50
50
0.001
10
8
10
-
9The 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, Figures3 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.Figure 3
Bifurcation diagram forx˙s by increasing Ω.Figure 4
Bifurcation diagram forx˙r by increasing Ω and Lyapunov exponents for selected forcing frequencies.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(xs,ys) and the other one associated with rotor (xr,yr). Figure 5 shows a period-1 behavior related to Ω=131.7 rad/s, while Figure 6 shows a period-2 behavior related to Ω=132.2 rad/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.Figure 5
Stator (a) and rotor (b) phase space and Poincaré section atΩ=131.7 rad/s.
(a)
(b)Figure 6
Stator (a) and rotor (b) phase space and Poincaré section atΩ=132.2 rad/s.
(a)
(b)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Ω=133rad/s. Figure 7 shows 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 Figure 8. The periodicity of the previous results, presented in Figures 5, 6, and 8, is confirmed by the Lyapunov exponents, whose values are all negative.Figure 7
Basin of attraction atΩ=133 rad/s.Figure 8
Coexisting attractors atΩ=133 rad/s: (a) period-5 orbit; (b) period-4 orbit.
(a)
(b)More complex responses are now in focus considering other kinds of behaviors. Initially, the rotational frequency is increased toΩ=135rad/s. Figure 9 presents 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}.Figure 9
Stator phase space and Poincaré section atΩ = 135 rad/s showing a quasiperiodic response.Figure 10
Detail of the Poincaré section atΩ = 135 rad/s.By increasing even more the rotational frequency toΩ=140rad/s, the system presents a chaotic response shown in Figure 11 in stator phase space together with its Poincaré section. Figure 12 presents only Poincaré section of rotor and stator subspaces. The Lyapunov exponents obtained are λ={+11.4;-3.4;-9.5;-18.2;-25.6;-30.9;-38.5;-57.0}, confirming the chaotic behavior.Figure 11
Phase space and Poincaré section atΩ = 140 rad/s.Figure 12
Poincaré sections atΩ = 140 rad/s: (a) stator; (b) rotor.
(a)
(b)


## 5. Conclusions

This paper deals with the analysis of a nonsmooth rotor-stator 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.



---

*Source: 8478951-2017-10-09.xml*