Determination of the global responses characteristics of a piecewise smooth dynamical system with contact

In this paper the global response characteristics of a piecewise smooth dynamical system with contact, which is specifically used to describe the rotor/stator rubbing systems, is studied analytically. A method to derive the global response characteristics of the model is proposed by studying each piece of the equations corresponding to different phases of the rotor motion, i.e., the phase without rubbing, the phase with rubbing and the phase of self-excited backward whirl. After solving the typical responses in each phase and deriving the corresponding existence boundaries in the parameter space, an overall picture of the global response characteristics of the model is obtained. As is shown, five types of the coexistences of the different rotor responses and deep insights into the interactive effect of parameters on the dynamic behavior of the model are gained.

Phase angle μ Coefficient of friction τ Non-dimensional time ω Rotating speed of the rotor ω 0 Backward whirl "natural" frequency of nonlinear coupled rotor/stator system ω 2 Natural frequency of the rotor system with zero clearance 1

Introduction
Rotor-to-stator rubbing is a malfunction in rotating machinery that degrades the machine performance and may lead to a catastrophic failure of a whole machine through dry whip in the worst case. A large number of research papers deal with the dynamics of rub-related phenomena in rotor/stator contact systems. Some papers studied the possible rubbing responses induced through the rotor/stator rubbing and tried to derive the relationship between the rubbing responses and the system parameters, as summarized in [1]. Through the numerical and experimental studies, many rubbing behavior in rotor/stator contact systems have now become well known, for instance, the jump phenomenon and the synchronous full annular rubs [2][3][4], the partial rubs in sub-and super-synchronous whirl [5][6][7], the partial rubs in quasi-periodic whirl [8,9], the chaotic motion [10,11] as well as the dry friction backward whirl which includes two regimes: of dry whirl and dry whip [15][16][17][18].
The numerical and experimental studies are useful in detecting the possible rubbing responses in rotor/stator contact systems and in capturing some trends of dependence of the responses on system parameters. But a huge amount of work is inevitable in these ways in order to determine the boundaries of different rubbing responses of the rotor/stator contact systems. Furthermore, it was usually found that different researchers might present quite different response cascades and the inherent reason was not well explained. For instance, it was observed in tests of [12,13] that the rotor under the unbalanced excitation would follow a normal behavior with the increase or decrease of the rotating speed, that is, the rotor follows the response cascade: no-rub motion → a synchronous full annular rub → no-rub motion, and vise versa. Jump phenomenon was detected in both run-up and rundown process. However, it was reported that when an additional outside disturbance (hitting on the shaft) was applied on the rotor during the process, even at a very low rotor speed, the dry friction backward whirl with super-synchronous whirl frequency would be triggered from a synchronous response. A quite different response scenario with the increase of the rotor speed was observed in the tests of [3,16] as: no-rub motion → a synchronous full annular rub → a partial rub → dry whip (heavy rub). So dry whip was triggered solely by the mass unbalance for the rotor/stator systems in this case, which is different from the case mentioned above where an outside disturbance was needed. This indicates the necessity to further explore the existence boundaries of different responses of rotor/stator contact systems so as to get an overall picture on the interactive effect of system parameters on the appearance of different response cascades. The recent analytical works [4,17,18] have provided some satisfactory explanation on the response characteristics of rotor/stator contact systems mentioned in above tests.
This paper aims to depict the global response characteristics of a specific rotor/stator contact system by determining the existence conditions of the most typical rotor responses of the piecewise smooth nonlinear system. As known, the piecewise smooth nonlinear model consists of three pieces of equations of motion governing three phases of rotor motion, namely, the phase without rub, the phase with rubs and the phase in self-excited dry friction backward whirl. As will be shown in this paper, the existence condition for the most often observed responses, i.e., synchronous no-rub motion, synchronous full annular rub, partial rub and dry friction backward whirl, can be derived solely from a single piece of the equations even though some responses, i.e., the partial rubs, run across two phases of the rotor motion. After incorporating the existence boundaries of all of the responses together into the parameter space, an overall picture of the response characteristics of the rotor/stator contact system is achieved. Deep insights into the interactive effect of different parameters on the system responses as well as the coexistence behavior are gained. Moreover, the analytical methods used in this paper can be also applied to the more complex rotor/stator systems, i.e., the model considering the stator dynamics and the contact stiffness at contact surfaces.
The paper is organized in the following ways. In Sect. 2, the equations of motion for each phase of motion of the rotor/stator contact system are first introduced. Then the equations for different phases of motion are studied individually in Sect. 3 in order to derive the solutions and the existence boundaries of the corresponding responses. Section 4 will illustrate the global response characteristics of the present rotor/stator contact system in the parameter space through examples. Finally in Sect. 5, the conclusions from this work are drawn.

Mathematical model
The rotor/stator contact system consists of a Jeffcott rotor and a stator that is modeled as an added stiffness (see Fig. 1). This model physically describes the rotor in contact with a non-rotating susceptible circumferential stator or a mechanical seal. The rotor has a weightless shaft carrying a disk with mass m at the middle of the span. The mass center of the rotor is located at a distance e from its geometrical center. The stiffness of the shaft is k s . The clearance between the rotor and the stator is denoted by r 0 . The stator, which is rigidly fixed, has an elastic contact surface modeled as a symmetrical set of radial springs with isotropic The whirl frequency of the shaft and the rotor radius at the axial position of contact are denoted by ω w and r disk respectively. In the present analysis the gravity is neglected.
The governing equations of the rotor/stator contact system can be written as where r = x 2 + y 2 and v rel represents the relative velocity between the rotor and the stator at the contact point. Θ is Heaviside function with Θ = 0 when r < r 0 and Θ = 1 when r ≥ r 0 . Equation (1) implies that only the sliding contact is assumed to occur between the rotor and the stator. As can be seen from (1) the equations of motion, in fact, consist of three pieces that govern three phases of rotor motion: (i) the phase without no-rub when Θ = 0; (ii) the phase with rubs Θ = 1 (but v rel > 0 all the time); (iii) the phase of self-excited backward whirl (Θ = 1 and v rel fluctuates around zero where the dry friction effect must be taken into account). The first two phases of motion are quite familiar and widely used in the previous works to investigate the dynamics of the rotor/stator rubbing. The last phase of motion exhibits the well-known rubbing response: dry friction backward whirl, which contains two regimes-of dry whirl and dry whipas defined in [16] and which is not really well explored.
Equation (1) will be first rewritten in the nondimensional form: 3 where the non-dimensional variables are defined as: The prime indicates the differentiation with respect to the dimensionless time variable τ = ω 2 t .
By defining W = X + jY and |W | = √ X 2 + Y 2 , (2) can be written in a complex form. Rather than writing a unified complex equation, we will write the equations into three pieces of equations according to the phases of motion.
The piece of equation governing the phase without rubbing (|W | < R 0 ) is given by The piece of equation describing the phase with rubbing (|W | ≥ R 0 ) can be written as The piece of equations describing the phase of selfexcited backward whirl (|W | ≥ R c with the sign of V rel alternating) can be written as As is seen from the above discussion, the dynamics of the rotor/stator system in different phases of motion are governed by different equations of motion depending upon the amplitudes of the rotor deflection. The equation in each phase of motion is smooth. So we can solve each piece of equation(s) to get the corresponding solution and then carry out the stability analysis to determine its existence boundary. As will be shown, the most often observed responses in the rotor/stator systems, such as the synchronous no-rub motion, the synchronous full annular rub motion and the dry friction backward whirl, can all be determined in this way. Although the quasi-periodic partial rub motion alternates between the phase without rubbing and the phase with rubbing, its existence boundaries can still be determined by investigating the piece of equation governing the phase with rubbing. After the existence boundaries of the solutions from all the pieces of equations are determined, the global response behavior of the system in the parameter space can be depicted.

The solution in the phase without rubbing
In this phase, the rotor is free from the contact with the stator, the synchronous no-rub motion is the only response and can be derived from (3). The periodic solution can be written as The periodic solution is unconditionally stable. However, the solution is physically valid only when its amplitude is smaller than the clearance, namely A < R 0 , due to the existence of stator. To find the parameter regions where the synchronous no-rub motion exists, let A = R 0 in (6) and we get the parameter equation defining the boundary of no-rub motion: Equation (7) yields two pairs of real roots, each with an equal magnitude and opposite sign. Since the rotating speed is supposed to be positive and the synchronous no-rub response always undertakes forward whirl, two positive roots, Ω l and Ω u , are chosen to define the boundaries of no-rub motion. That is, the synchronous no-rub motion exists when the rotating speed Ω is less than Ω l or greater than Ω u .

The solutions in the phase with rubbing
Synchronous full annular rub motion is a response by which the rotor is in the full annular contact with the stator and whirls forward with constant amplitude at a frequency equal to the rotating speed of the rotor. Equation (4) can be used to determine this response by supposing a proper form of the solution according to the characteristics of the synchronous full annular rub. Let us denote the periodic solution as W 1 = Be j (Ωτ +ϕ) , two nonlinear equations of B and ϕ are yielded after substituting W 1 into (4) and separating the real and imaginary parts of the obtained equation, To eliminate ϕ from the above equations, we square both equations of (8) and then sum them up to get a second order polynomial equation of B in the form One can solve for the positive real roots of B from (9) and then substitute them back into either of the two equations of (8) to get ϕ for the given set of parameters. From the condition with rubbing, B ought to be greater than R 0 . The response of the synchronous full annular rub corresponds to the stable periodic solution of (4). So the stability of the solution must be analyzed. To carry out the stability analysis on the periodic solution of the nonlinear system, the physical understanding of the response is of good use. It is known that a periodic solution in the inertial coordinate system will correspond to a fixed point in the rotating coordinate system at the same rotating speed as the period of the solution. To transform (4) into the rotating coordinate system, we use the transformation W 1 (τ ) = ρ r (τ )e jΩτ . Substituting it into (4) and making some manipulations, yields Equation (10) has a fixed point in the form of ρ r0 = Be jϕ with B and ϕ being defined by (8).
Further, the stability analysis on the fixed point, and thus on the stability of the synchronous full annular rub solution in the inertial coordinate system, will be carried out. To accomplish this task, the complex amplitude of the response is written in the component form ρ r = X r + jY r , where the subscript "r" indicates the variable in the rotating coordinate system.
By defining the state vector in the rotating coordinate system as U = {X r Y r X r Y r } T , (10) can be written in a time-invariant differential equation where A and g( U ) are respectively the coefficient matrix and the nonlinear vector term: (11) about the solution of the fixed point U 0 = {B cos ϕ B sin ϕ 0 0}, yields a linear, timeinvariant differential equation where δU = U − U 0 is the perturbation solution to the fixed point. J( U 0 ) is the so-called Jacobian matrix which is given in the following form: From Jacobian matrix given in (13), one can set up the characteristic equation of the solution, The real parts of the characteristic roots decide the stability of the fixed point, and thus the stability of the synchronous full annular rub solution. Since the stability boundaries are of interest, bifurcation theory is being used as is done in [4]. The saddle-node bifurcation condition, by which a zero characteristic root exists, defines the existence boundary of the solution. The Hopf bifurcation condition, whereby a pair of pure imaginary characteristic roots appears, gives the boundary between the quasi-periodic partial rub and the synchronous full annular rub.
Quasi-periodic partial rubbing is a response of the rotor/stator system, by which the rotor rubs the stator intermittently. There are generally more than one frequency component in the system responses during partial rubs. The equations governing the partial rubs should be a combination of (3) for the phase without rubbing and (4) for the phase with rubbing. Generally, we have difficulty in obtaining the analytical solutions of the piecewise smooth equations. Without the solution of partial rub, it seems impossible to determine their existence boundaries.
Fortunately, we can still sketch the existence boundaries of quasi-periodic partial rub based on the study of other responses and on the understanding of the physical reason for the loss of stability of the partial rub. For instance, it is easy to determine the boundary between no-rub motion and partial rub. Since the synchronous full annular rub bifurcates into the partial rub through Hopf bifurcation, the Hopf bifurcation condition of synchronous full annular rub solution will serve as the boundary between the synchronous full annular rub and the partial rub. As observed in experiment [3,14] and simulations [4], the partial rub ceases to exist after the dry friction backward whirl was triggered by the imbalance with the increase rotating speed. In [17] the physical reason for the onset of dry friction backward whirl was exploited. It is known that the rotor in resonance at a backward whirl "natural" frequency of the coupled nonlinear rotor/stator system is the reason for the onset of the dry friction backward whirl from the partial rub. From this understanding, a method can be devised to determine the approximate onset condition-another existence boundary of the partial rub.
First, the backward whirl "natural" frequency of the coupled nonlinear rotor/stator system is sought. This frequency corresponds to the mode that the rotor whirls backwards. After setting the right-hand side of (4) to zero and then substituting supposed solution, W = H exp(j ω 0 τ ), into it, we solve the frequency with a negative value and the amplitude (H should be greater than R 0 ) as Secondly, when the rotor runs at a rotating speed near the magnitude of the backward whirl "natural" frequency, (4) is transformed into a form that describes the dynamics of the rotor/stator system near the backward whirl "natural" frequency, where the left-hand side of the complex system has the natural frequency equal to the above backward whirl "natural" frequency. Thirdly, an approximated linear equation to (16) is set up after defining as a small parameter. The reason is that |H | is the amplitude of the rotor in the backward whirl rubbing without imbalance excitation and |W | is the amplitude of the rotor in the backward whirl rubbing with imbalance excitation. Since Ω < 1 in the vicinity of the critical rotating speed of the onset of dry whip, the magnitude of the excitation Ω 2 in (4) is thus a small quantity. So |W | should be only a small modification to |H | and the difference between the corresponding inverse relations should be also small. So (16) can be written approximately as: The approximated equation can be solved by, e.g., the Multiple Scale Method [19], since W at the rotating speed near the magnitude of the backward whirl "natural" frequency, namely Ω = |ω 0 | + εσ where σ = O(1), is a small detuning.
Assume that the problem depends on many time scales: The solution of (18) can be represented by an expansion having the form The expansion to O(ε 2 ) is carried out in this paper and the corresponding equations of ε 0 and ε 1 are . By this moment, the range of the excitation frequency (the rotating speed) Ω is not specified as done in the analysis procedure of the Multiple Scale Method at the Secondary Resonance [19]. The solution of (20) yields where a(T 1 ) and θ(T 1 ) are the real functions of T 1 and

Substituting (22) into (21) leads to
where the prime denotes the derivative with respect to T 1 . Now, substituting the above expression of Ω into (23) and setting the coefficient of terms with exp(j ω 0 T 0 ) on the right-hand side of (23) to zero in order to eliminate the secular term, leads to a complex equation Separating the real and the imaginary parts and then solving a and aθ , gives where C 1 = R 0 2(ω 2 0 +ζ 2 ) , C 2 = R 0 (1 + μ 2 )(ω 2 0 + ζ 2 ) × |A|, α = arctan ζ ω 0 + arctan μA r +A i A r −μA i . A r and A i are the real and imaginary parts of A.
Next, (25) is transformed into an autonomous system with T 1 not explicitly appearing by defining and substituting into (25) to get The steady-state motions occur when a = γ = 0, which corresponds to the singular points of (27). They correspond to the solution of |W 0 | with a and θ by (28) and (26) can serve as an approximation of |W | in (17). By using |H | = H in (15), (17) can be calculated for the given parameters and ε. Finally, the critical rotating speed for the onset of the dry friction backward whirl or the loss of stability of the partial rub can be solved after the deviation εσ to |ω 0 | can be determined. To do this, we can first set σ a reasonable constant number, i.e., σ = 1 as used in this paper or a value that depends on the system parameters as given in [17], then use a trial and error process on ε in order to find the most proper ε that makes D f (as given above) a largest negative value.

The solution in the phase of self-excited backward whirl
Dry friction backward whirl is a motion state of the rotor/stator system, by which the rotor is in continuous contact with the stator, rolling or slipping continuously on the contact surface and whirling backward at a super-synchronous frequency. Dry friction backward whirl is a self-excited vibration of the rotor/stator contact system due to the presence of dry friction effect at the contact surface. Although the mechanism of dry friction backward whirl is well understood, most of the research works used the governing equations, which do not include the dry friction effect to explore this phenomenon. A little progress has been made in analytically determining the characteristics and the parameter region of dry friction backward whirl until now. Below, we will only sketch the solution procedure.
To make the reader to better understand the meaning of (5) and the methodology we adopted to approximate the amplitude of dry friction backward whirl, the mechanism for a rotor/stator system that performs dry friction backward whirl is first illustrated.
When rotor-to-stator contact occurs, the friction force at the contact point tends to drive the rotor into backward whirl. In case that the rotor whirls backward, the friction force will apply in the same direction as the rotor whirl and continuously turns the input rotational energy into the energy for the lateral vibration of the rotor system. The vibration amplitude of the rotor is driven higher and higher. The relative velocity at the contact point will change its sign when the rotor deflection is sufficiently large to make the whirl speed exceed the circumferential speed of the rotor at the contact point, i.e., r = −ωr disk /ω w (refer to Fig. 1(b)) or ω b |W | ≥ ΩR disk (see (5)). At this time, the friction force also inverts its direction and is opposite to the rotor whirl to dissipate the energy. The rotor deflection will be driven down. When the rotor deflection decreases to make the relative velocity at the contact point to change sign again, the rotor deflection will be driven up again. In this way, the amplitude of rotor deflection at dry friction backward whirl is kept oscillating around the value with a zero relative velocity, that is, the rotor undertakes a dry friction induced self-excited vibration-dry friction backward whirl. Now we can use the knowledge of dry friction backward whirl to derive the solution from (3). Since the fluctuation of the response amplitude around the zero relative velocity is much smaller than the total magnitude of the rotor deflection during dry friction backward whirl, the amplitude of dry friction backward whirl can be well approximated by the one at zero relative velocity. From the second equation in (5), it yields By substituting (29) into first equation of (5), we obtain the equation governing the motion of dry friction backward whirl in the following form: To derive solution for dry friction backward whirl from (30), two important features should be taken into account. (a) Since the contribution of the forced vibration to the total response at dry friction backward whirl is generally too small to be considered, the supposed solution needs only to consider the frequency component with super-synchronous frequency. (b) Most importantly, when the governing equation changes into the form of (30), the solution corresponding to dry friction backward whirl should have a time increasing amplitude, that is, a positive real part in the exponent, because the dry friction effect does not explicitly take effect (its effect is taken into account in (29) to give a constant amplitude). So the solution can be assumed in the form After substituting (31) into (29), one can get two equations of α and ω b , where ω b indicates the whirl frequency at dry friction backward whirl and α the capability of the system to undertake self-excited vibration. When α is positive, self-excited vibration-dry friction backward whirl-can take place. When α is negative, the responses other than self-excited vibration occur. Therefore, α = 0 is the critical condition for the existence of dry friction backward whirl [18]. In this way the boundary for dry friction backward whirl in the parameter space can be determined.

Global response characteristics
In this section, two examples are presented to illustrate the overall picture of the global dynamic response characteristics of the rotor/stator contact system with ζ = 0.05, R 0 = 2.00 and R disk = 20R 0 . As can be seen from the above analysis the existence boundaries of no-rub motion, the synchronous full annular rub and the dry friction backward whirl are determined rigorously with very high accuracy. The numerical simulations have proved this point. For the onset condition of dry friction backward whirl from partial rub motion, some approximations are adopted. The validity of the results is also checked by the numerical simulation as given in [17]. So the global dynamical characteristics of the rotor-to-stator contact system will be well demonstrated in the following analysis. In the following figures, curves Ω l and Ω u indicate the rotating speed whereby the deflection of the linear rotor covers the clearance. Thus, the synchronous no-rub motion is in the rotating speed range less than Ω l and greater than Ω u . Curve HP and line SN are respectively the boundaries of Hopf bifurcation and saddle-node bifurcation of the synchronous full annular rub solution. Curve DW is the boundary by which dry whip is triggered by imbalance and the partial rub ceases to exist. Curve DF gives the existence boundary of the dry friction backward whirl. Also, the parameter regions marked by numbers represent different regions of coexistence, i.e., "0"-no-rub motion + dry friction backward whirl; "1"-synchronous full annular rub + dry friction backward whirl; "2"-partial rub + dry friction backward whirl; "3"-synchronous full annular rub + no-rub motion; "4"-no-rub motion + synchronous full annular rub + dry friction backward whirl. Figure 2 depicts the global responses characteristics on the parameter planes of Ω-μ with β = 0.04. It is interesting to note that the dry friction backward whirl exists already from a very low rotating speed when the friction coefficient is larger than a certain value which satisfies the condition given in (15) as μ > 0.102. When the friction coefficient is below this value, the only rubbing response is the synchronous full annular rub solution. Due to the coexistence of the no-rub motion and the synchronous full annular rub (see region "3"), the "jump" phenomenon will occur in this region. When the friction coefficient is above this value, one can find the regions that the dry friction backward whirl coexists with other responses. So it is quite possible to induce the more destructive dry friction backward whirl from a relatively mild rubbing response or even from no-rub motion when the rotor is subjected to a sudden outside disturbance. It is even worse that there is a region (enclosed by DW and Ω u ) where dry whip is the sole stable response. Thus, it is important to avoid a large contact friction, especially that which exceeds the one defined by the criterion in (15). The coexistence of three responses can be also detected in a very small region-"4". Figure 3 shows the global response characteristics on the parameter planes of Ω-β with μ = 0.15. It is easily seen from Fig. 3 that with the increase of the stiffness ratio (it can be taken as the reduction of the stator stiffness but keeping the rotor stiffness unchanged, namely, using a soft stator), both Ω l and Ω u increase but the latter increases more rapidly. This makes the rotating speed range with the synchronous full annular rub larger and the speed range with the no-rub motion smaller. This seems to be unprofitable to the rotor response. However, it is found from Fig. 3 that when the stiffness ratio is larger than a certain value, i.e., greater than about 0.3 in this case, the 9 Fig. 3 Rotor response characteristics on parameter plane of Ω-β with μ = 0.15, where ζ = 0.05, R 0 = 2.0 and R disk = 20R 0 . The regions of only three of the five typical types of coexistence are found in this case. Ω l and Ω u are the critical rotating speeds of no-rub motion. SN and HP are the saddle-node and Hopf bifurcation boundaries of the synchronous full annular rub. DF and DW are the existence and onset boundaries of dry friction backward whirl synchronous full annular rub will no longer bifurcate into the more severe partial rub. More importantly, the induction of the destructive dry whip through the imbalance will not occur. In fact, when the stiffness ratio is greater than about 2.4, the imbalance can no longer trigger the dry whip although the rotor undertakes the partial rub. From this sense, a larger stiffness ratio benefits the rotor response. Since the criterion given in (15) is met in the range of the stiffness ratio examined in Fig. 3, the dry friction backward whirl exists in the whole region right to line DF. This means that the dry friction backward whirl coexists with all other rotor responses we have examined. From the point of view of the global responses, there are only three types of coexistence regions in this case, except for region "3" and region "4" mentioned above.

Conclusions
In this paper, the global response characteristics of a piecewise smooth rotor/stator contact system are investigated. An effective method to get the global response characteristics of this specific model has demonstrated that the piecewise smooth equations of motion can be analyzed piece by piece so that the existence boundaries of the typical responses in each phases of motion, i.e., in the phase without rub, with rub and of self-excited backward whirl, have been determined. Then, they are incorporated into the same parameter space to get the global response characteristics of the rotor/stator contact system. The results show deep insights into the interactive effect of different parameters on the system responses as well as the coexistence behavior. The main conclusions can be summarized as follows.
From the point view of global responses, five typical types of coexistence are found for the piecewise smooth rotor/stator contact system in this paper. They are the coexistence of no-rub motion and synchronous full annular rub, of no-rub motion and dry friction backward whirl, of synchronous full annular rub and dry friction backward whirl, of partial rub and dry friction backward whirl and of no-rub motion, synchronous full annular rub and dry friction backward whirl. Except for the first type of coexistence, the other types of coexistence are not quite familiar and occur, in fact, when the criterion given in (15) is met. From this criterion and also from diagrams of the global response characteristics in the parameter planes, it is known that small friction on the contact surfaces and large stiffness ratio of the rotor-to-stator stiffness will benefit the rotor rubbing behavior by avoiding the occurrence of the severe partial rub or the destructive dry whip.