Mikhael Tannous

Patrice Cartraud

Mohamed Torkhani

David Dureisseix

Assessment of 3D modeling for rotor-stator contact simulations

come

Introduction

An accidental blade loss causes a rotor unbalance which leads to rotor-to-stator contact while the turbine is slowing down and crosses a critical speed.

Rotor-stator contact can lead to diversified consequences (cf. [START_REF] Muszynska | Rotor Dynamics[END_REF]). It may be permanent like a full annular rub with a forward synchronous whirl, or with a backward whirl. This later can be a rolling without sliding called a dry whirl (or a dry friction backward whirl) and can also be a rolling (with friction) with a sliding in the opposite direction of the spinning (dry whip). The intensity of the friction forces is the main factor behind the dry whip phenomenon [START_REF] Bartha | Dry Friction Backward Whirl of Rotors[END_REF]. Therefore, the friction coefficient [START_REF] Choi | Investigation on the whirling motion of full annular rotor rub[END_REF], an the rigidity of the rotor and of the stator are the main factors [START_REF] Choy | Non-linear transient analysis of rotor-casing rub events[END_REF] involved.

The rotor-stator contact can be intermittent, i.e., with rebounds. Rebounds can be periodic, quasiperiodic or chaotic [START_REF] Popprath | Nonlinear dynamics of a rotor contacting an elastically suspended stator[END_REF].

The eigen modes of the rotor-stator system may change due to the coupling that occurs during the interaction phase [START_REF] Muszynska | Rotor Dynamics[END_REF][START_REF] Chu | Stiffening effect of the rotor during the rotor-to-stator rub in a rotating machine[END_REF]. Local zone plastifications may also lead to eigen mode changes.

Thermo-mechanical effects are also observed. In fact, the rotor-stator contact energy is transformed to heat that causes local expansions leading to supplementary unbalance and reduces the rotor-stator gap distance [START_REF] Millecamps | Influence of thermal effects during bladecasing contact experiments[END_REF]. Generally, most of this heat is received by the rotor causing it to deflect (the Newkirk phenomenon) [START_REF] Goldman | Thermal bending of the rotor due to rotor-to-stator rub[END_REF]. Thermo-mechanical effects can lead to spiral vibrations. This has been highlighted on a Jeffcott rotor in Sawicki et al. [START_REF] Sawicki | Thermomechanical behavior of rotor with rubbing[END_REF] and on real machines in Bachschmid et al. [START_REF] Bachschmid | Thermally induced vibrations due to rub in real rotors[END_REF].

However, if the rotor is performing a free rotation after a blade loss, slowing down can occur. Roques et al. [START_REF] Roques | Modeling of a rotor speed transient response with a radial rubbing[END_REF] remains one of the few research studies in the literature to deal with the slowing down of a turbine then leading to rotor to stator contact when passing through a critical speed. The authors highlight the effect of the friction coefficient, and of the stator rigidity on a turbine slowing down behavior. An analytical model was used in Braut et al. [START_REF] Braut | Dynamic analysis of the rotor-stator contact due to blade loss[END_REF] to study the Jeffcott rotor slowing down due to its contact with a suspended rigid stator.

It is obvious that the rotor-stator contact problem is a complex and highly nonlinear problem, presenting both multiphysical (vibrations, contact, thermo-mechanical effects, etc.) and multi-scale (local deformations, etc.) phenomena that are complicated to model and to take into account correctly. In the literature, assumptions are made and it is difficult to assess their validity, even if comparisons with experimental results are provided. As mentioned previously, rotor slowing down and Dry Whip depend highly on the friction torque and the rigidity of the system. However, the first models in the literature to deal with rotor-stator contact problems are based on a simple Jeffcott rotor [START_REF] Childs | Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis[END_REF]. Moreover, the rotational velocity is constant, and is kept so along the rotor-stator contact phase by a compensating torque. These simplified models are described by differential equations and solved analytically, neglecting the gyroscopic effects [START_REF] Karpenko | Piecewise approximate analytical solutions for a Jeffcott rotor with a snubber ring[END_REF].

One of the main drawbacks of these simplified models is the rigid body assumption for the rotor and stator modeling. A way to improve the model accuracy while limiting the computational coast consists in using 1D models based on beam finite elements. This has been proposed recently in Peletan et al. [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF] for the analysis of an industrial EDF turbine with a harmonic balance method (see also [START_REF] Von Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF]). The rotor to stator contact has also been studied in Roques et al. [START_REF] Roques | Modeling of a rotor speed transient response with a radial rubbing[END_REF] with the same type of beam model. The transient response is calculated with a time integration scheme. Though industrial examples have been considered in these two previous papers, the approximations resulting from beam assumption and thus the beam to beam contact for the rotor to stator contact modeling are difficult to assess. The aim of this research study is to assess the need for more realistic models, i.e., 3D models for rotor-stator contact problems. This is an up-todate issue if we examine the literature on the blade/casing contact problems. In fact, blade/casing contact problems involve physical aspects and require modeling techniques that present important similarities with rotor-stator contact problems [START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines[END_REF]. Most of the research developments in the literature use beam models to simulate blade/casing contact interactions. However, the need for 3D models for a better understanding of blade/casing contact has recently appeared [START_REF] Legrand | Full three-dimensional investigation of structural contact interactions in turbomachines[END_REF]. Indeed, local plastifications, surface coating extraction, local thermo-plastic effects, etc. cannot be seen unless a 3D model is used. The contribution of 3D modeling to blade/casing problem simulations raises the question of the validity of simplified beam models for rotorstator contact problems. In fact, rotor-stator contact problems involve more complex modeling and physical aspects than those of a blade/casing contact problem. For instance, a mistaken friction torque estimation can lead to a dry whip behavior, in spite of a forward synchronous whirl. Such misleading conclusions can occur on simplified rotor-stator contact models since such models are not capable of a rigorous presentation of the contact surfaces and the physical aspects of surface to surface contact and friction.

This paper aims at illustrating the contribution of a 3D finite elements modeling of the rotor-stator contact interactions and consequently the limits of a simplified 1D model, by comparing the results of both models. Numerical examples are constructed to be consistent with industrial EDF turbines. However some data have been changed for confidentiality purposes.

1D and 3D simulations are carried out with the open source finite element software Code Aster [START_REF] Edf R&d | Code Aster: a general code for structural dynamics simulation under gnu gpl licence[END_REF].

The paper is organized as follows. Firstly, the contact algorithms and the time integration techniques adapted to the problem are exposed. Then, the rotor-stator system is presented, and three application examples highlight the main contributions of a 3D model to rotor-stator contact problems.

The contact problem and the solution techniques

The contact/friction problem is a highly nonlinear complex problem. Solving a contact problem requires first to solve a geometric problem relative to the position of the two contacting bodies, and then to solve the incremental contact problem that enables one to obtain the contact pressure from which depends the friction force [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. The computational cost and the stability of the solution depend also on the choice of the time integration technique. distance between the master and the slave surfaces, respectively Γ 1 and Γ 2 . The nodes of Γ 2 are projected on Γ 1 and the distance between the two surfaces is computed as

The geometric and optimization problem

g n = min x1∈Γ1 x 1 -x2 (1)
x 1 is a node that belongs to Γ 1 and x2 ∈ Γ 2 is the result of the following optimization problem that is equivalent to projecting x 1 on Γ 2 according to the outgoing normal direction n of Γ 2 :

x 1 -x2 = min x2∈Γ2 x 1 -x 2 (2)
Three main cases are distinguished: g n > 0: no contact, g n = 0: perfect contact and g n < 0: interpenetration. Although the third condition is not a physical one, it may occur for some contact algorithms such as the penalty method (cf. Section 2.2). The normal contact condition is the non-penetration of the two contact bodies known as the Kuhn-Tucker-Karush [START_REF] Minoux | Mathematical Programming: Theory and Algorithms[END_REF] formulation:

g n ≥ 0, p n ≥ 0, p n g n = 0 (3)
where p n is the contact pressure. The tangential behavior (friction) generates a tangential force F t that is often simplified by a Coulomb friction law by adopting a friction coefficient µ that depends only on the material and surface properties:

F t ≤ µp n (4)
In a master-slave formulation, the master surface is often the one that possesses the higher rigidity. However, if both surfaces are of the same rigidity, Noels [START_REF] Noels | Contributions to energy-conserving time integration algorithms for non-linear dynamics[END_REF] proposes a double face algorithm that alternates master and slave surfaces. The master-slave formulation is adapted to a rotor-stator contact problem and enables us to avoid interpenetration if both surfaces are meshed with quadratic elements.

In this paper, the 3D rotor-stator contact problems are based on a master slave formulation with a quadratic mesh for both the rotor and the stator.

Choosing the contact solution method

A penalty contact method is the only solution technique proposed in 1D contact modeling in Code Aster. For the 3D model, Lagrange multipliers method [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and Lagrangian methods[END_REF] can be used. However, in Code Aster, this latter leads to an unaffordable computational cost. It also has convergence problems if the contact conditions are severe. Therefore, penalty contact algorithms are used in both 1D and 3D contact problems. Using the same method makes the two modeling techniques consistent and more appropriate to the objective of comparing their corresponding results (see Section 5).

The penalty contact method transforms the contact problem to a non-constrained optimization problem. The problem size is not increased by contact condition; however, its major drawback is that it allows interpenetration. In fact, this method imposes a contact force that is proportional to the interpenetration value and of opposite direction: F P ext = K n |g n | where K n ≥ 0 is a penalty coefficient. The higher K n , the more accurate the method is. If K n → ∞, the penalty method is equivalent to the physical problem described in Fig. 3 but high penalty coefficient values lead to tangent matrix conditioning problems and to algorithm divergence. Therefore, interpenetration cannot be avoided. Similarly, the tangential behavior of the contact can be penalized by a coefficient K t ≥ 0.

For a non-adhesive contact problem, the fundamental dynamic equation, according to the penalty method, reads

Mü + C u + Ku + Γ1 (K n g n δg n + K t g t • δg t )dΓ 1 = F ext (5)
If the contact is adhesive, Fig. 5 becomes

Mü + C u + Ku + Γ1 (K n g n δg n + F t • δg t )dΓ 1 = F ext (6)
A parametric study should be performed to verify the independence of the solution from the penalty coefficients K n and K t .

Choosing the time integration technique

Implicit integration schemes in impact problems usually present small interest, since their unconditional stability that allows us to use large time steps is not a real profit, since small time steps are required in order to accurately account for the physical phenomena that are taking place. Thus, explicit integration techniques are more adapted to contact problems, with a fewer cost per time step.

Moreover, for nonlinear dynamics, and particularly, for contact problems "conventional" integration schemes, such as the Newmark integration scheme present conditional stability. Unconditionally stable integration schemes were developed for such applications. We mention the Constraint Energy Momentum Algorithm (CEMA) method proposed by Hughes et al. [START_REF] Hughes | Finite-element method for nonlinear elastodynamics which conserves energy[END_REF] as well as the Modified Energy-Momentum Method (MEMM) method (see [START_REF] Kuhl | Constraint energy momentum algorithm and its application to non-linear dynamics of shells[END_REF]). The EMCA method (energy momentum conserving algorithm) was introduced by Simo and Tarnow [START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF] and modified in Romero and Armero [START_REF] Romero | High-frequency dissipative time-stepping algorithms for the dynamics of nonlinear shells[END_REF] by the introduction of a numerical dissipation, which lead to the EDMC method (energy dissipation momentum conserving). This latter is still one of the best for contact problems.

However, standard finite element softwares have not yet integrated such integration schemes, since their formulation depends on each study case. A Newmark integration scheme is the unique choice for 3D contact problems in Code Aster. The Newmark integration scheme is not energy conservative for contact/friction problems [START_REF] Géradin | Mechanical Vibrations: Theory and Application to Structural Dynamics[END_REF]. The energy conservation is achieved in a frictionless contact problem by a non-penetration kinematic condition expressed on the velocity [START_REF] Moreau | Numerical aspects of sweeping process[END_REF], also called the persistence condition [START_REF] Chawlar | Energy consistent algorithms for frictional contact problems[END_REF]. Moreover, high frequency oscillations appear and may lead to divergence, when the Newmark scheme is used to solve contact problems. However, if the time step is small enough and numerical dissipation is introduced, the Newmark scheme is stable and conserves the angular momentum [START_REF] Hughes | Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics[END_REF][START_REF] Kane | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical system[END_REF] for contact problems. This was performed in this paper for 3D contact problems and resulted in a large computational cost.

Moreover, the integration schemes used for 1D rotor-stator contact problems, and that for 3D contact problems cannot be the same on Code Aster. In fact, the integration scheme that offers good stability for 1D contact problems is an explicit non-damped Runge-Kutta integration scheme [START_REF] Van Der Houwen | Stability of collocation based Runge-Kutta-Nystrom methods[END_REF]. It offers a good stability for contact problems and a time step adaption technique that enables us to use small time steps in the vicinity of severe contact conditions, and larger time steps for noncontacting conditions. In our study case, the initial time step is 5 × 10 -5 s and goes to 5 × 10 -7 s in the vicinity of contact conditions. However, the main drawback of this integration scheme is that it is memory consuming and 8 GB RAM were needed to solve the 1D contact problems.

The rotor-stator system under study

In this paper, an example has been constructed from an accidental blade off in a nuclear power plant turboset. The corresponding unbalance is introduced under the form of a local force which induces rotor vibrations.

The considered system is depicted in Fig. 2, the disk being described as a rigid body. The upper part of Fig. 2 shows the upper half of the rotor's longitudinal cross section. The length and the radius of the rotor are 20.15 m and R = 1.075 m respectively.

Table 1: Dimensions and materials

The stator consists of two rings linked by a number of blades. In order to simplify the model, and as done in the thesis of Roques et al. [START_REF] Roques | Modeling of a rotor speed transient response with a radial rubbing[END_REF], the exterior ring is not modeled and the extremities of the blades are fixed. Two stator models are considered and are shown in the lower part of Fig. 2. For the first one, the stator is a simple ring, and for the second one a more realistic model is used, i.e., a ring that is suspended on 8 blades1 . The first stator model makes it easier to analyze the rotor-stator contact problem (controlling the interpenetrations, the rotor displacements, observing local deformations on the stator, etc.) and to make suitable conclusions, before considering the more complex and realistic example of the suspended stator.

The stator's center is located at the rotor's axis (on point N c in Fig. 2). The dimensions of the stator as well as the rotor's and stator's materials are found in Table 1.

A 6.8 × 106 N unbalance force that follows the rotor's rotations is applied to the simply supported rotor at point N b . It is triggered by the loss of a 100 kg mass at a 2.75 m distance from the rotating axis, when the rotational velocity of the latter is 1500 rpm.

The rotor rotational velocity is a data of the problem. Its time evolution is governed by Eq. ( 7) in Section 4.2, which leads smoothly to a constant value. Due to the rotating unbalance, the rotor vibrations are amplified until a rotor-stator contact is detected.

One current limitation of Code Aster is that the rotational velocity is a data of the dynamic problem. Its evolution law is not an unknown of the mechanical problem. Therefore, rotor-stator contact will not lead to a deceleration of the rotor, but rather, will result in stator local deformations, rotor torsions, etc.

Models of the rotor-stator contact problem in an industrial code

Few results in the literature deal with 3D rotor-stator contact problems. This may well be due to the somehow high computational cost engendered by such models. To reduce this computational cost, the 1D and 3D dynamical model coupling was the main concern in Tannous et al. [START_REF] Tannous | A beam to 3D model switch for rotor dynamics applications[END_REF] but no 3D contact modeling was addressed at that time.

This paper proposes a first 3D contact model between the rotor and the stator, implemented within an industrial FE code, namely Code Aster, whose features have been enhanced to address a complex 3D rotor-stator problem, consistent with the previous 1D contact model.

Contact model suited to the 1D geometry description: the impact law

The penalty approach is used for both the 1D and the 3D contact problems in the scope of this paper. In this section, we present the geometric approach that is available for 1D contact problems in Code Aster. It is a simple and efficient 1D modeling approach (cf. [START_REF] Alarcon | R5.06.03 Modelization of the Shocks and Friction in Transient Analysis by Modal Recombination[END_REF]), based on a node-to-node contact formulation.

The equations of motion that take into account the gyroscopic effects are settled in the 1D model rotating frame. The involved impact law has been improved recently to take into account the computation of the frictional torque between the rotor and the stator. Dedicated to the 1D problem, this model is presently one of the best 1D contact approaches, due to its computational efficiency and accuracy.

It is illustrated in Fig. 3. Considering a simple stator (a ring without blades, whose inner radius and thickness are respectively R and ), its beam model is a set of curved beam elements along a circle whose radius is R + /2. Concerning the rotor model, each cross section is represented in the beam model by a single node. It may potentially impact all the different nodes of the stator. The precision of the model indeed directly lies on the number of nodes on the stator, as well as does the computational cost. Each contact condition takes into account a rigid disk around the potential contacting nodes to represent the correct distance between the impacting parts. The disk positions are adjusted along the computation when the rotor and stator beam elements are deforming. Nevertheless, the model does not take into account the cross-section deformation.

An explicit integration scheme of the Runge-Kutta family is the optimal choice that has been selected in Code Aster for these 1D rotor-stator contact problems.

3D modeling of the rotor-stator contact problem within Code Aster

Contrary to the previous model for which the rotation velocity is prescribed by the user, in the 3D model, it is taken into account by prescribing movement of the end sections of the rotor. In this study, the rotor is subjected to a null rotation acceleration at t = 0, that is afterwards increased in such a way that it produces a given velocity ω at time t m with a null acceleration (to prevent transient effects):

ω(t) = ω × 3 -2 × t t m t 2 t 2 m (7)
These end sections are far enough from the contact area so that the different possibilities in the prescribed movement do not influence the solution close to the contact locations.

Contact conditions are dealt with a master-slave formulation with penalty, and the only available integration technique for 3D impacts in Code Aster is a damped Newmark scheme. This last choice may not be optimal in our case since it leads to high computational costs. Implementation of other schemes is a development perspective.

Model consistency

The two previous models appear to be consistent since the only differences concern the contact methods and their solution algorithms. Once converged, their solutions therefore can be compared. Since the study case involves a rotor which is slender enough, Fig. 2, the beam element model is expected to be accurate.

The 3D rotor is meshed with quadratic elements involving 1400 nodes (see Fig. 7). This mesh provides sufficient precision and reduces the difference between the beam and the 3D rotors to its maximum. The same 3D rotor model is used in the subsequent application examples. Note that a finer mesh (tests have been performed on a 2000, 3000 and 4000 node 3D rotor meshes) does not lead to a higher similarity between the beam and the 3D models.

Both solutions exhibit a similar behavior. As Table 2 illustrates it, a small but inevitable difference exists between the natural frequencies of the beam rotor model and those of the 3D rotor model. Their small values nevertheless allow us to conclude that the two rotor models are equivalent. Fig. 4 depicts the eleventh natural mode of the rotor.

Moreover, the 1D and 3D rotors have similar transient behaviors when submitted to a dynamic loading: Fig. 5 shows the orbits of point N c in the (x, y) plane due to an unbalance loading without a rotor-stator contact. It is obvious that the two rotors exhibit a similar behavior. This is also the case for the velocities2 and accelerations. The 1D and the 3D rotor models are therefore consistent.

However, the industrial stators are not slender enough to be modeled by beam elements. The ring itself, with the dimensions of Fig. 2, can present a 10 percent frequency difference, even on the lowest natural frequencies, between the 1D stator model and even the mostly refined quadratic 3D rotor mesh. More details can be found in the PhD thesis of Tannous [START_REF] Tannous | Development and evaluation of coupled 1D and 3D numerical approaches for rotor-stator contact modeling[END_REF]. Therefore, the 1D contact model will consist of a rotor modeled by beam elements, while the stator is modeled by 3D quadratic elements. The contact between the 1D rotor and the 3D stator is solved via the impact law (cf. Section 4.1) and is called a 1D contact model. The 3D contact model will consist of a contact between a 3D rotor and the same 3D stator as the one used for the 1D contact problem. The contact is solved via a penalty contact algorithm with a master-slave formulation.

For the 3D contact model, three sections of the rotor are subjected to contact with the stator. For the 1D contact model to be consistent with the 3D one, three contact sections are chosen on the rotor (three different shapes are used to point out the three different beam nodes in Fig. 6). Each one consists of a beam node that may enter in contact with a series of stator nodes situated at the same z-axis position and having the same shape as their rotor contact node in Fig. 6.

According to the impact law, the rotor nodes are affected by rigid disks of radius R (the same radius as the 3D rotor) while the stator nodes possess zero radius rigid disks.

Moreover, since the impact law and the penalty contact method (for 3D contact models) use penalty coefficients that influence the values of the contact and friction forces, the same penalty coefficients are used for both the 1D and the 3D contact models. Note that, since in Code Aster the rotational velocity is a data of the dynamic problem and since no direct and reliable access to the contact forces and torques is possible, then no slowing down due to the rotorstator contact is studied, and the comparisons between the 1D and the 3D contact models are performed on the basis of the rotor orbits and the stator behavior. The main contributions of a 3D model are highlighted.

Application examples

Three study cases are considered. In the first one, the rotational velocity of the rotor is 240 rpm and the stator consists of a simple ring (without blades, see Fig. 2). The stator external surface is fixed. Such a contact case is relatively simple to analyze. Interpenetrations are easily detected, since the rotor displacements are limited to the gap size (8 mm) added to the stator deformations, which are practically negligible in such a case.

We then model a more complex and realistic rotor-stator contact case at the same rotational velocity, in which the stator is suspended by blades. The rotor-stator contact and friction lead to blade deformations and cause stator rotations. The third and last study case is dedicated to a rotor-stator contact, when the rotor velocity is close to its critical velocity.

In all three study cases, the rotor starts rotating with null velocity and acceleration and reaches its constant rotating speed at t m = 0 : 01 s according to Fig. 7.

Rotor contacting a stator without blades

We consider the 3D rotor-stator system in Fig. 7 whose dimensions are those of Fig. 2. The stator consists of a simple ring with no blades, and with fixed external surface. Quadratic 3D elements are used3 : the rotor counts 1400 nodes and the stator counts 259 nodes. Master slave contact formulation is considered, with 136 nodes on the stator contact surface (slave surface) and 80 nodes on the rotor contact surface (master surface).

For the 1D model, we recall that the stator model is the same as in the 3D model. But the rotor is discretized with beam elements and counts 41 nodes.

Since penalty contact algorithms are used, a parametric study is performed, in the following, to show the negligible influence of the penalty parameters on the 1D and 3D contact results. The highest penalty coefficients, for which convergence is reached, are chosen, so that the interpenetrations are minimized and can be neglected. These coefficients are denoted k n and k t , respectively for normal and tangential penalty parameters. For this case of study, k n = 10 1 4 N/m and k t = 10 9 N/m. Three contact cases are studied. In the first one, k ns = k n /10 and k ts = k t /10. In the second, the nominal parameters are considered, i.e., k n and k t and in the third one k nf = k n × 10 and k tf = k t × 10 are selected. If the three simulations lead to similar results, then the influence of the penalty coefficients is neglected. This is performed on both the 1D and the 3D contact models.

Fig. 8 shows a comparison of the orbits of point N c belonging to the 1D main rotor axis at the contact section location. It is obvious that the penalty coefficients have negligible influence on the rotor orbits. Thus, the chosen penalty parameters enable us to have a good physical representation of the contact phenomenon. The interpenetrations are negligible and do not exceed a maximum of 0.01 mm (note that the rotor-stator gap distance is about 8 mm).

The same parametric study is also performed on the 3D model, and for the same penalty coefficients. Results are presented in Fig. 9, and show a negligible sensitivity to the penalty parameters, a sensitivity that is even less than the one observed on the 1D contact problem. For the 3D model, each contact case requires a 30 h computational time on a 2.5 GHz QuadCore with 10 GB RAM machine, and for a time step equal to 5 × 10 -5 s. 3D contact problems are, therefore, CPU time consuming as illustrated in this study case, especially when an implicit integration technique is used.

In the following a comparison of the 3D rotor orbits (at point N c ) with the 1D rotor orbits is set. Fig. 10 shows an obvious difference between the 1D and the 3D contact models. It can be seen that before contact both models have similar responses, as it has been shown also earlier in a case without contact, see Fig. 5. Therefore the differences observed in Fig. 10 are only due to contact modeling and not due to the models themselves. For a better understanding of the difference sources, let us check the stator behaviors of the 1D and the 3D contact problems.

The stator of the 3D contact simulation at t = 0.875 s is amplified by a 1.5 × 10 5 factor and is shown in Fig. 11. The stator shows a local deformation at the contact with the rotor. However, if the stator of the 1D simulation is amplified by the same factor and shown at the same contact instance (see Fig. 12), no local deformation is found.

In fact, the rotor-stator contact in a 1D simulation engages a wider contact surface than a 3D contact simulation (α ≥ β in Fig. 13) due to the rigidity of the rotor contact cross-section in a beam model assumption. Thus, no local deformation of the stator is observed in a 1D rotor-stator contact model. This is, in a part, the source of difference between the 1D and the 3D contact models.

To numerically demonstrate this conclusion, we plot a 1.5 × 10 4 times amplified deformation of the 3D rotor crosssection in the contact zone at t = 0.875 s as shown in Fig. 13. The triangles represent the different cross-section nodes. This latter is no more a disk, and this can easily be seen if the cross section is compared to a disk having the same center as the cross-section and passing through the node that is the most far from the circle center. The cross-section deformation reduces the rotor to stator contact surface with respect to that of a rigid cross-section assumption in a beam model. This conclusion highlights the necessity of a 3D contact model for rotor-stator contact simulations, and is to be endorsed in the following study cases.

Rotor contacting a suspended stator

In this study case, the stator is suspended by blades and is presented in Fig. 14. Its dimensions are those of Fig. 2.

Fig. 15 shows simultaneously the 1D and the 3D rotor, as well as the 3D stator that is used in both 1D and 3D rotor-stator contact simulations.

In this study case the stator is meshed quadratically and counts 2280 nodes. 128 nodes belong to the contact surface 4 .

The rotor is spinning at 240 rpm. The time evolution of this rotation speed is smooth. It is described in Fig. 7 and presented in Fig. 16 with t m = 0.01 s. This allows us to minimize the transient stage due to initial conditions or when the rotation speed becomes constant since the accelerations are null at start up and at t = t m . The time step is 5 × 10 -5 s.

For a one second simulation period, the computational time is 60 h on a 8 GB RAM and a Quad2-Core 2.75 GHz CPU computer.

Since the penalty contact method is chosen, a parametric study is performed on the normal and tangential penalty coefficients to ensure that penetration is negligible and results are not dependent on these coefficients. The nominal values k n = 10 14 N/m (normal contact) and k t = 10 10 N/m (tangential contact) are thus taken, since it appears that multiplying or dividing these coefficients by 10 do not lead to solution changes. Now the results of the 1D and 3D contact models are compared. First, the rotor behaviors are examined. Fig. 17 shows a comparison between the orbits of point N c of both the beam and the 3D models.

Differences on the rotor orbits are small. In fact, a stator suspended by blades absorbs some of the impact energy and reduces local stator deformations.

We now examine the stator behavior. Fig. 18 shows a comparison between the displacements of point C 1 (see Fig. 2) for a beam and a 3D modeling of the rotor-stator contact.

Larger displacements are obtained (up to three times) for the beam modeling of the stator when compared to the 3D model, along the y-axis direction. This direction is indeed the tangential direction that is influenced by the friction force. This is also illustrated when checking the tangential displacement of point C 2 (in the x-axis direction) as shown in Fig. 19, but with a less extend in the normal contact direction, as in Fig. 20 with the displacements of point C 1 along the x-axis direction. In fact, there are two types of errors: one in amplitude, and the other one in phase. The first one is quite large. The second one is small as it can be seen at the beginning of the simulation, which results from the frequency consistency of both models. However, even though this phase error is small, it increases over the simulation and may lead to out of phase results between 1D and 3D models. Moreover, the tangential displacement of the 1D stator has a gradually varying amplitude. The change from one amplitude to another coincides with impacts, also seen on the rotor's orbit in Fig. 17.

Note that the 1D contact problem is solved via a non-damped Runge-Kutta integration scheme, while the 3D problem is solved via a damped Newmark integration scheme. It is not possible to use the same integration scheme and damping for both the 1D and the 3D contact problems. Therefore, some of the differences between the 1D and the 3D solutions is due to using different integration techniques. For example, we observe high frequency oscillations on the 1D solution in Figs. [START_REF] Edf R&d | Code Aster: a general code for structural dynamics simulation under gnu gpl licence[END_REF] and 18 for t ∈ [0.7, 0.8] s, and this is not the case for the 3D solution. These oscillations are not a source of error and do not lead to divergence. However, along the contact normal direction, the influence of the time integration scheme choice is clearer. In fact, the 3D solution is getting damped for t ∈ [0.6, 0.9] s and that is not the case of the 1D one as seen in Fig. 20 which shows the displacements along the x-axis of point C 1 for both the 1D and the 3D contact simulations along the contact normal direction. For a better clarity, the 1D and 3D solutions are also presented separately at the right hand side of Fig. 20. One can say that the difference between the two solutions is obvious and the use of different integration schemes makes this difference wider, but cannot be the main cause of the difference.

As a conclusion, the differences between the beam and 3D contact models arise from the contact force, both with its normal and tangential components, but with a predominance of the tangential one, so for the rotor-stator friction torque that differs between the two models. The main source of discrepancy arises from the rotor cross-section stiffness, leading to a wider contact surface for the beam model than for the 3D model, as discussed in Section 6.1 and illustrated in Fig. 13. To assess it, Fig. 21 depicts the deformed stators, with an amplification factor of 100, at t = 0.876 s. The stator rotation around its main axis is clearly larger for the 1D model than for the 3D one, as the blade deflections exemplified it. The rotational velocity being an input of the dynamic problem in our study case, the rotor-stator impact does not involve a rotor deceleration, but leads to a rotational deflection of the stator. In the converse case, when the rotational velocity is an unknown of the problem, its deceleration would have been larger for the 1D contact model than for the 3D one. This motivates the necessity of using a 3D model in the case of the analysis of a turbine slowing down due to an accidental condition. This 3D model improves the physical accuracy of the rotor-stator impact problem.

Rotor contacting a suspended stator in the vicinity of a natural frequency

In this example, the same rotor-stator system and parameters are used. The only difference is that the rotational velocity of the rotor is set to ω = 636 rpm, corresponding to 10.6 Hz which is the first natural frequency of the rotor. We first check that the behavior of both the 1D and the 3D rotors at this rotational velocity is consistent and that both the rotors are excited in the same manner. We then check the penalty parameters and ensure that the results are not dependent on them. The retained values are k n = 10 14 N/m and k t = 10 10 N/m. A 0.2 s simulation is performed. The computational cost of the 1D simulation does not exceed 3 min. However, the 3D computational cost is much higher. Due to the contact severity, the time step needed for the convergence dropped from 5 × 10 -5 s to 5 × 10 -7 s. The contact severity actually increases due to the vicinity of the natural frequency. This fact is highlighted in Fig. 22 that represents the orbits of point N c with and without a rotor-stator contact. The contact severity combined with an implicit solution via the Newmark scheme leads to a prohibitive computational cost: 170 computational hours on a Quad-Core 2.5 GHz machine with 12 GB RAM were needed for the 3D computations.

The friction torque is clearly overestimated as illustrated in Fig. 23, which represents the state of the 1D and 3D stators at t = 0.1567 s. Note that no amplification is needed to observe the deformation of the stator in this study case.

This application example highlights the increase in the computational cost of a 3D model when the rotational velocity approaches a natural frequency of the rotor. The 1D model presents non-negligible estimations that lead to an overestimation of the friction torque, and a stator behavior that is very different from the 3D one.

Conclusions

This paper aimed to compare 3D and 1D models of rotor to stator contact. The 1D model obviously relies on more restrictive assumptions than the 3D model, mainly because of the use of a beam model which is based on a rigid-body kinematics for the cross-section. Therefore, the 3D results can be considered as reference results which provide a mean to assess the validity of the 1D model. Simulation examples were carried out on industrial turbines inspired from the TAG of EDF nuclear parks and on the industrial finite element code Code Aster. The main conclusions are recalled in the following.

In a 3D model, the rotor cross-section deformation leads to a narrower rotor to stator contact interface, and thus leading to local deformations and a stator behavior that is explicitly different between a 1D and a 3D model. Moreover, the rigid rotor cross-section in a 1D model causes an overestimation of the friction torque. The friction torque is the main factor that decelerates a free rotating rotor, causes backwards whirl, etc. Thus, overestimating the friction torque leads to mistaken conclusions and raises the question of the validity of 1D rotor-stator contact models that dominate the literature and are basically based on rigid Jeffcott rotors and rigid stators. Such models suffer from limitations to describe real industrial problems.

Moreover, deciding whether a 3D model is imperative is not only a question of friction torque overestimation by a 1D model. The 3D model remains mandatory in spite of its high computational cost, if the stator is deformable and the contact is severe (in the vicinity of a natural frequency for example) in order to accurately account for the stator local deformations neglected by a 1D model due to the rigid rotor cross-section assumption.

The results showed also that, in spite of a better accuracy, the increase of the penalty coefficients makes the convergence less likely and increases the computational cost. In 1D modeling, it contributes to an increase in the memory storage needed for the solution.

This paper underlines the low performance of the 3D models, which in part is due to the time integration technique. One can use an explicit scheme. Another way to reduce the computational time consists in using a strategy which enables us to switch from the 1D model to the 3D model when contacts occur, as proposed by the authors in Tannous et al. [START_REF] Tannous | A beam to 3D model switch for rotor dynamics applications[END_REF]. This enables us to combine the advantages of 1D and 3D models. The 1D model is actually accurate enough without contact, while the 3D model is required in the nonlinear part of the simulation.

Figure 1 :

 1 Figure 1: Contact between two bodies.

Figure 2 :

 2 Figure 2: Dimensions of the rotor-stator system under study. (a) The dimensions of the stator (in meters) Ring Blades Inner radius R i 1.083 Length l 0.755 Thickness 0.207 Wideness L 0.038 Thickness along the z-axis (of Fig. 2) e = 0.15 (b) The materials of both the rotor and the stator Density (kg/m 3 ) 7800 Young Modulus E (Pa) 2.1 × 10 11 Poisson coefficient 0.3 Coulomb friction coefficient between the rotor and the stator µ 0.2

Figure 3 :

 3 Figure 3: The choc law basics (the dots are the different stator nodes).

Figure 4 :

 4 Figure 4: The eleventh natural mode of the rotor. (a) Beam rotor model. (b) 3D rotor model.

Figure 5 :

 5 Figure 5: Orbits of point N c of the rotor.

Figure 6 :

 6 Figure 6: 1D contact model, between a beam rotor model and a 3D stator, via the impact law.

Figure 7 :

 7 Figure 7: The 3D rotor-stator system of this study case.

Figure 8 :

 8 Figure 8: Influence of the penalty coefficients on the 3D contact problem.

Figure 9 :

 9 Figure 9: Influence of the penalty coefficients on the 3D contact problem.

Figure 10 :

 10 Figure 10: Orbits comparison of the 3D and the 1D rotor-stator contact formulation.

Figure 11 :

 11 Figure 11: 3D deformed stator amplified by a 1.5 × 10 5 scale at t = 0.875 s.

Figure 12 :

 12 Figure 12: 1D deformed stator amplified by a 1.5 × 10 5 scale at t = 0.875 s.

Figure 13 :

 13 Figure 13: 3D rotor's cross-section amplified deformation in the contact area at t = 0.875 s.

Figure 14 :

 14 Figure 14: 3D model of the rotor-stator system.

Figure 15 :

 15 Figure 15: 1D and 3D models of the rotor-stator system.

Figure 16 :

 16 Figure 16: The time evolution of the rotation speed ω(t).

Figure 17 :

 17 Figure 17: Orbits of the 3D rotor and the beam rotor.

Figure 18 :

 18 Figure 18: Tangential displacements of point C 1 (the y-axis direction).

Figure 19 :

 19 Figure 19: Tangential displacements of point C 2 (the x-axis direction).

Figure 20 :

 20 Figure 20: Normal displacements of points C 1 .

Figure 21 :

 21 Figure 21: Comparison of the deformation of 1D and 3D stator models, amplified by a factor 100.

Figure 22 :

 22 Figure 22: Orbits of points N c in a 1D model.

Figure 23 :

 23 Figure 23: State of the stator at t = 0.1576 s.

Table 2 :

 2 3D and beam rotors natural frequencies comparison.

	No	Natural frequencies (Hz) Mode nature
		Beam rotor 3D rotor
	1-2	23.77	23.75	Bending
	3-4	62.48	62.29	Bending
	5	79.85	79.87	Torsion
	6-7	115.33	114.69	Bending
	8	128.68	128.79	Traction/Compression
	9	159.71	159.86	Torsion
	10-11 178.1	176.66	Bending

The master-slave formation[START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] remains one of the most popular methods for solving the geometric problem. Solids S 1 and S 2 in Fig.1are susceptible of contact. The geometric problem consists in calculating the

The industrial stator has a larger number of blades. However, reducing the number of blades simplifies the model and reduces the computational cost while having minor influence on the objectives and conclusions of this research work.

The velocities are important since the friction torque depends on the contact velocities. The 1D and the

3D rotors velocities are very close.

A sufficiently refined 3D mesh with quadratic elements is suited to obtain reliable results. In our study case, refining the mesh leads to non-negligible computational cost for the same result precision.

Note that a more refined mesh does not have a remarkable effect on the results, while increasing significantly the computational cost.

Acknowledgments. The authors thank the French National Research Agency (ANR) in the frame of its Technological Research COSINUS program (IRINA, project ANR 09 COSI 008 01 IRINA).