Development of beam-to-beam contact detection algorithms for rotor-stator rubbing applications

In nuclear power plant turbosets, the design-basis accident consists of a blade-off on the low pressure turbine last stage. During the accidental shutdown, a severe rotor-casing interaction may occur at critical speeds due to large shaft line displacements originated by a high unbalance excitation. The contact between the shaft and the stator, also called the diaphragm in this study, induces an important angular deceleration rate and greatly modifies the turbogenerator dynamics including the amplitude of the loads in the bearings. Therefore the main objective is to verify that the designed turbine is capable of going through critical speeds without catastrophic consequences for the shaft line. To this end, a model of a turbogenerator has been developed to compute rotor speed transients by considering the rotating speed of the rotor as an unknown, which allows for the angular deceleration due to rubbing to be calculated in a more realistic fashion. Lagrange multipliers method is applied to compute contact forces. The diaphragm, which is a non-rotating bladed disks assembly, is modeled by curved and straight beams, and different assumptions for the contact detection are studied to find a compromise between CPU time and accuracy. Results of the numerical tool show that the contact forces are sensitive to the retained assumptions only when heavy rub occurs. Nevertheless, the rotating speed and bearing loads are computed with a satisfactory accuracy, even with an approximation on contact detection that saves CPU time. INTRODUCTION In turbomachinery, contact between rotating and stationary parts of the turbine, usually referred to as rub [1, 2], is known to be a serious malfunction. Rub is usually caused by unbalanced rotor, blade-off, misalignment of the rotor centerline, thermal unbalance, casing excitations, to name a few. In this study, the shaft-diaphragm interaction only is considered: indeed, in real turbosets, the blades-to-diaphragm contact [3] is negligible in comparison to the shaft-stator interaction to compute the angular deceleration due to friction and the torsional stresses in the coupling components. For a detailed overview of rubbing phenomena and its associated literature, the reader is referred to the very thorough survey given by Muszynska [4]. During the past decades, numerous authors have developed and enhanced mathematical models characterizing the rotorstator interaction. First mathematical models are the Jeffcott rotors [5] with the following assumptions: linear stiffness and damping, rigid disk on a massless shaft, Coulomb law as friction modeling, stator modeled by springs acting in the radial direction and sometimes with dampers. Moreover, due to its simplicity in programming, the penalty method is often used to compute contact forces. Then, with modern computers, the behavior of flexible rotors has been studied by finite element approach and/or modal synthesis techniques [6] allowing more realistic models. However these studies were limited to steady state analysis so that the angular velocity remains constant, which implies an in-


INTRODUCTION
In turbomachinery, contact between rotating and stationary parts of the turbine, usually referred to as rub [1,2], is known to be a serious malfunction.Rub is usually caused by unbalanced rotor, blade-off, misalignment of the rotor centerline, thermal unbalance, casing excitations, to name a few.In this study, the shaft-diaphragm interaction only is considered: indeed, in real turbosets, the blades-to-diaphragm contact [3] is negligible in comparison to the shaft-stator interaction to compute the angular deceleration due to friction and the torsional stresses in the coupling components.For a detailed overview of rubbing phenomena and its associated literature, the reader is referred to the very thorough survey given by Muszynska [4].
During the past decades, numerous authors have developed and enhanced mathematical models characterizing the rotorstator interaction.First mathematical models are the Jeffcott rotors [5] with the following assumptions: linear stiffness and damping, rigid disk on a massless shaft, Coulomb law as friction modeling, stator modeled by springs acting in the radial direction and sometimes with dampers.Moreover, due to its simplicity in programming, the penalty method is often used to compute contact forces.Then, with modern computers, the behavior of flexible rotors has been studied by finite element approach and/or modal synthesis techniques [6] allowing more realistic models.However these studies were limited to steady state analysis so that the angular velocity remains constant, which implies an in-1 crease of the driven motor torque in the case of rub.The equations of motion have been rewritten considering that the speed law is given in order to compute rotor speed transient [6].Anyway, the main limitation of these simulations lies in the fact that the angular deceleration which is originated by contact between the shaft line and the stator is not taken into account.To the authors' knowledge, only Dai [7] has investigated a rigid rotor model where an additional equation computes the instantaneous rotating speed depending on the contact forces.This model can reproduce the rotor response with partial as well as full annular rubbing.Finally, other works included the coupling between lateral and torsional vibrations for a steady state.Edwards et al [8] concluded that the torsional phenomenon should be included in rubbing modeling.
Therefore, the primary goal of the present paper is to give a better representation of the rotor-stator interaction during speed transients.To this end, a numerical tool [9] has been developed to compute the transient dynamical response of the shaft line with shaft-to-stator contact.In this model, the rotating speed of the rotor is considered as an unknown, which allows for the angular deceleration due to rubbing to be calculated in a more realistic fashion.In conjunction with the Lagrange multipliers approach where friction is considered, the governing equations of the coupled rotor-stator system are solved using a time-stepping procedure based on the explicit central differences scheme [10].
The paper is organized as follows.The first section presents the modeling of the two structures used in the numerical tool.The contact dynamics is investigated in the second section where assumptions are used during the contact detection.Finally, the different time integration results corresponding are compared.

STRUCTURAL MODELS
In this section, the FE modeling of the two structures and the general algorithm to solve speed transients are presented.

SHAFT LINE
A 1-D turbine system is built by modeling the shaft with straight spinning Timoshenko beams, bladed disks with circular rigid disk and imbalances with concentrated masses.An academic example is depicted in figure 1.The originality of the modeling lies in the fact that the angular position is considered as an unknown of the problem.Gyroscopic effects are taken into account which couple the torsional vibrations and the lateral ones.
More details concerning the model are given in [9].

Rigid disk
The kinetic energy E D c of an axi-symmetric rigid circular disk rotating at an angular speed P ' takes the form: where M D , I DX and I DZ respectively stand for the mass, the second moment of inertia and the polar moment of inertia of the disk.Referring to figure 2, .u;Â v / denotes the bending displacements in the .X; Z/ plane, .v;Â u / the bending displacements in the .Y; Z/ plane, w the traction along Z axis and ˇ, the torsion in Z direction.The derivative of x with respect to time is denoted P x.Within the small perturbation framework, the kinetic energy is rewritten as follows: Equation ( 2) must be exact up to the third order so that all second order terms are retained in equations of motion.It is worth noting that the axial vibrations only are independent.
Shaft The kinetic energy E c of the shaft corresponds to the integration of the disk energy along its longitudinal direction: where , S , I x and I p respectively refer to the mass density, the cross-section area, the moment of inertia and the polar moment of inertia of the beam.The potential energy E d of a spinning Timoshenko beam is given by the following integral: where E, G and k respectively stand for Young's, shear moduli and the transverse shear form factor.
Bearings In real turbines, oil film bearings are used to support the shaft.The dynamical study of a shaft line requires to describe the nonlinear behavior of the oil film.As a first approach, stiffness and damping coefficients of the oil film are approximated by linearizing Reynolds' equations with respect to the equilibrium position.Moreover a comparison of results between linear and nonlinear modelings [11] has shown that the behavior of the crushed oil film, when passing through critical speeds, can be accurately approximated by a very stiff bearing.The equations of motion derived from Hamilton's principle are nonlinear even without considering any contact.

DIAPHRAGM
A flexible diaphragm (see figure 3) has been developed to enrich the modeling of the interaction: the inner ring is modeled according to the theory of curved beams and the blades are considered as straight beams.
Inner ring: curved beams The inner ring, which initially belongs to the .x;y/ plane, is depicted in figure 4. It is very thin so that the Euler-Bernoulli theory can be used.Consequently, the deformations of the curved beams [12,13] are written as follows: (5) where u c and v c correspond respectively to radial and tangential displacements of the point that belongs to the centroidal line as The rotations in the 3-D space are introduced in figure 4. The path variable is denoted by s and R c is the average radius of the inner ring.Considering that the radial distance r of a point that belongs to the cross-section of the curved beam is negligible in comparison with the radius R c , " ss in equation ( 5) becomes: (6) Then considering the linear elastic behaviour of the casing, i.e.
ss D E c " ss , rs D G c rs and zs D G c zs (Hooke's law) where E c and G c are the Young's and shear moduli of the casing's material, the strain energy E c d for a curved beam element, ie an angular sector, is obtained by integrating on the elementary volume V c : where I cr , I cz , J c and S c respectively stand for the two moments of inertia, the polar moment of inertia and the cross-section area of the curved beam.The kinetic energy is equal to: with I cp denoting the moment of inertia.
Blades: straight beams The Euler-Bernoulli theory is also used for the blades which are slender structures.Their energies correspond to a particular case of equations 3 and 4: the angular speed is equal to zero and since Euler-Bernoulli theory is considered, the shear energy vanishes and the two usual kinematic relationships are used.Then, according to the notations introduced in figure 5, it comes: where I b z , I b y , I b x , J b and S c respectively stand for the two second moments of inertia, the polar moment of inertia, the crosssectional polar moment of inertia and the cross-section area of the rectangular straight beam.
A proper modeling of the connection between the curved beam elements and the straight beam elements is required.It is achieved through the following relationships coming from Euler-Bernoulli assumptions of curved beams: Equations of motions are derived in a very general fashion for a 3-D space but, in what follows, the diaphragm is represented by a planar model for the sake of simplicity.Equivalently, the conditions w c D 0 and Â s D 0 are considered.

CONTACT DETECTION
The forces of particular interest in this study are the contact forces acting between the shaft and the diaphragm where friction is considered.The following assumption for the contact detection is used: only one point of the rotor cross-section, which is supposed rigid, comes into contact with the diaphragm.

General theory
To describe the contact dynamics, the master-slave approach is used.Then, for any point of coordinates x belonging to the master interface, the gap function g.x/ between the two bodies can be stated as follows: where g 0 denotes the initial gap, n .m/ the unit vector tangent to the surface of the master body (ie outward pointing normal), u .m/.x/and u .s/.N x/ respectively the displacement of the point that belongs to the master interface and the one of the point of the slave surface whose coordinates N x minimizes the interpenetration.
When the Lagrange multipliers method is used, the contact conditions [14] can be seen as the Kuhn-Tucker optimality conditions [15]: all material point x that belongs to the master interface must satisfy the following equation: N 0, g.x/ 0 and N g.x/ D 0 where N stands for the positive contact force acting on the slave surface in the normal direction.The gap function must be positive because two bodies cannot interpenetrate: this is also known as the impenetrability condition.To these unilateral contact conditions, the Coulomb friction law is added when only sliding occurs: for which is the coefficient of friction, v T the tangential slip velocity and T , the contact force acting on the tangential direction.

Application to rotor-stator system
In this study, the master body is the rotor and the slave the stator.Then the coordinates N x of the casing must be determined.However, as the nonlinear detection of contact greatly increases the CPU time, different assumptions are made: 1. contact approximation (Hyp 1): the contact location on the casing (point C1 in figure 6) has the same angular position as the rotor geometric center; 2. exact contact (Hyp 2): N x, referred as point C2 in figure 6, is computed by minimizing the gap function g.x/ which is analytically calculated for any material point of the diaphragm by considering the exact deformed shape of the stator.
Once the contact location is obtained, the gap function is linearized to construct the contact constraint matrix in the normal direction.The contact constraint matrix in the tangential direction is calculated by considering that only sliding occurs [3].The radial contact force component N is computed by using the Lagrange multipliers method, which guarantees the impenetrability condition.The tangential contact force T , which is acting opposite to the direction of the relative sliding velocity, stems from the Coulomb friction law with an amplitude equal to j T j D j N j.
Finally, the friction torque takes the form C f ric D RF N .

GENERAL ALGORITHM
To the author's knowledge, the forward increment Lagrange method [10] gives satisfactory results for such contact-impact problems [3].Consequently this method is used in our study: a finite central differences scheme combined with Lagrange multipliers (prediction-correction algorithm) properly satisfies the contact detection and ensures the compatibility of the speed and the acceleration.Therefore, equations of motion are discretized in time according to central finite differences: with a D X or '.X n refers to the generalized displacements X in which bending, traction and torsion are stored at time t n and t to the time step.
Since the system of equations is non-linear, a specific time procedure is developed in algorithm 1.As the aim of this study is to compute the response of a low pressure turbine in abnormal conditions, the initial and boundary conditions are given as follows: operating at normal conditions, the turbine is suddenly disconnected after the blade-off.Simultaneously, only aerodynamical and Newtonian fluid friction slow down the shaft line [16].
Initialization of X and ' for t 0 and t 1 For n from 2 to n end do Prediction of (X Since the exact contact detection requires to update the contact location at each iteration of the nonlinear solver, the CPU time prohibitively increases.Therefore, in order to find a compromise between CPU time and accuracy, the contact location can be chosen to remain the one obtained during the prediction step (no update in the iterative process).

RESULTS
The shaft line mode used for simulations is depicted in figure 1.The initial gap between the shaft and the diaphragm is equal to 8mm.Convergence of results with respect to the time step and satisfaction of the impenetrability condition have been shown in [9].
The sensitivity to the proposed assumptions for the contact detection is now studied: the computed rotating speed and the displacements at the bearing location are depicted in figure 7 and 8. Table 1 summarizes the main results.From the beginning of the shutdown up to 3s, no contact occurs.The first contacts arise around 3:05s and correspond to a partial rubbing, as may be seen in figures 10 and 11 from the amplitude of the contact forces.During this phase of light rubbing, the results appear to be independent of the assumptions used for the contact detection, as shown in figure 13 for the available gap at the bearing node.
However, when heavy rub occurs, the sensitivity to the assumption on the contact detection becomes obvious.The simulation using the exact contact detection is considered as the reference calculation.Consequently, the contact approximation, referred to as "Hyp 1", gives approximated results on the evaluation of the bearing loads -see figures 10, 11 -and underestimates the contact forces, see table 1.If the exact contact location is not updated during the iteration process ("Hyp 2 no update"), the residual penetration increases, which leads to an over-estimation of the contact forces, see table 1 and figure 9.

CONCLUSIONS AND PROSPECTS
A mechanical system for rotor-stator interaction has been studied in this paper.The stator is a diaphragm and a flexible model has been proposed with curved beams representing the inner ring and straight beams representing the blades.A numerical tool has been developed to analyze speed transients for shaft-todiaphragm contact after a blade-off.Based on the Lagrangian multiplier method, different algorithms have been studied for the contact implementation in order to optimize the computation time.
In the case of partial rub, simulations have shown that results do not depend on the assumption made on the quality of the contact detection.However, since this study focuses on accidental shutdowns, full annular rubbing may occur.Results have shown that, when heavy rub occurs, the contact forces are then sensitive to simplifying assumptions for the contact detection.Nonetheless all these assumptions may be used to obtain other quantities of interest such as rotating speed or bearing loads with a good accuracy and saving CPU time.
Experiments on a test rig devoted to rubbing will be performed in the next future in order to validate this rotor-stator interaction modeling.The time integration results will be analyzed by means of the time-frequency tools.Finally, the diaphragm model will be extended by adding axial deformations and work is in progress to take into account the beam-to-beam contact in 3-D space.

Figure 1 .
Figure 1.SCHEMATIC OF THE ROTOR MODEL

Figure 3 .
Figure 3. SCHEMATIC OF THE DIAPHRAGM MODEL

Figure 7 .Figure 8 .
Figure 7. COMPUTED ROTATING SPEED WITH THE DIFFERENT ASSUMPTIONS FOR CONTACT DETECTION

Figure 11 .FERENTFigure 12 .Figure 13 .
Figure 11.CONTACT FORCES AT FIRST IMPACTS WITH THE DIF-FERENT ASSUMPTIONS FOR CONTACT DETECTION nC1 , ' nC1 ) Iteration to solve the nonlinear governing equations While (k Residual vector of equations k 2 Ä ") do If (No penetration) then