Jérôme Didier

Jean-Jacques Sinou

Béatrice Faverjon

Study of the non-linear dynamic response of a rotor system with faults and uncertainties

Keywords: rotor system, non-linear dynamic, uncertainties

In this paper, the quantification of uncertainty effects on the variability of the non-linear response in rotor systems with multi-faults (such as unbalance, asymmetric shaft, bow, parallel and angular misalignments) is investigated. To take account of uncertainties in this kind of non linear problem, it is proposed to use the Harmonic Balance Method (HBM) with a Polynomial Chaos Expansion (PCE). The efficiency and robustness of the proposed methodology is demonstrated by comparison with Monte Carlo Simulations (MCS) for different kinds and levels of uncertainties.

Introduction

During a phase of preliminary design of dynamic systems such as rotors, many physical parameters can have a significant influence on the vibration response of the system. Also, some features (such as unbalance, misalignment, bow, asymmetric shaft, etc) that can generate linear or non-linear responses need to be taken into account to avoid unexpected vibrational amplitudes and worst design for a rotor system. A state of the art of faults in rotating machinery can be found in [START_REF] Edwards | Fault diagnosis of rotating machinery[END_REF]:

• Rotor unbalance is one of the most current faults in a rotordynamics problem. This phenomenon is due to the fact that the center of mass of a disc or a shaft element is not on the axis of rotation of the system. It is generally modeled by an eccentric mass which generates a synchronous force.

• Misalignment in coupling is a common fault in rotating machinery, nevertheless this phenomenon remains incompletely understood. It can be noted that there are two cases of misalignment, firstly the parallel misalignment occurs when two shafts have a relative displacement, secondly the angular misalignment occurs when it exists an angular fault between two 1 shafts of a rotor system. The presence of a misalignment generates reactions forces and moments in the coupling which are source of vibrations. Consequently the dynamic behavior of the system depends of the characteristics of the coupling. An other distinction have to be done between flexible coupling (such as universal joint) and rigid coupling. Xu and Marangoni [START_REF] Xu | Vibration analysis of a Motor-Flexible Coupling-Rotor system subject to misalignment and unbalance, part i: Theoretical model and analysis[END_REF][START_REF] Xu | Vibration analysis of a motor-flexible Coupling-Rotor system subject to misalignment and unbalance, part II: experimental validation[END_REF] present a theoretical model and experimental validation for a flexible coupling. Theoretical development for a rigid coupling can also be found in [START_REF] Al-Hussain | Dynamic response of two rotors connected by rigid mechanical coupling with parallel misalignment[END_REF]. An other modeling has been developed in [START_REF] Lees | Some studies on misalignment in rigidly coupled flexible rotors[END_REF][START_REF] Lees | Misalignment in rigidly coupled rotors[END_REF]: Lees presents a simple model of a rigid coupling system with a parallel fault and investigates the influence of misalignment on the dynamic behavior of the system.

• The asymmetry is an other common fault in a rotor problem. Usually this fault can be due to the geometry of the shaft or due to the presence of a crack. Numerous studies have been done to take account of the presence of a crack, and the breathing behavior of the cracked rotor [START_REF] Penny | Simplified modelling of rotor cracks[END_REF][START_REF] Sinou | A non-linear study of a cracked rotor[END_REF]. Mathematical development for an asymmetric shaft cross-section has been developed by Oncescu et al. in rotating frame [START_REF] Oncescu | Investigation of the stability and steady state response of asymmetric rotors, using finite element formulation[END_REF] or by Lalanne and Ferraris in fixed frame [START_REF] Lalanne | Rotordynamics-Prediction in Engineering[END_REF].

• The bow is an initial deformation of the rotor which can be caused by different factors. In [START_REF] Pennacchi | Use of modal representation for the supporting structure in model-based fault identification of large rotating machinery: part 1-theoretical remarks[END_REF],

Pennacchi presents different sources of bow such as gravity effects on off-line machines or thermal effects. An asymmetric thermal distribution can be localized on an element (local bow) or all parts of the rotor (extended bow). Darpe [START_REF] Darpe | Dynamics of a bowed rotor with a transverse surface crack[END_REF] proposes to investigate the influence of a bow on the dynamic response of the rotor system with other kinds of faults such as unbalance and crack.

• Other kinds of faults exist, such as rubs or ovalization errors [START_REF] Richard | Model based fault identification in rotor systems by least squares fitting[END_REF], but they will not be studied in this paper.

As previously mentioned, all these faults can significantly reduce the life time of a rotating machine. Even if usually the main aim of research studies is to locate individual faults in machine operation, this study proposes to demonstrate the capability to include the treatment of uncertainties in multi-faults in the dynamic analysis of a rotor system, believing that this goal should be useful in the design stage of a machine. For example, some physical parameters may vary in an uncertain way during the manufacturing monitoring or the use process (like flight or others) and so the response may change also in one uncertain way. So, the novelty of the present paper is to propose a formulation of dynamical systems such as rotor with multi-faults, which includes uncertainties in physical parameters. Thus, the goal is to allow an estimate of the overall nonlinear dynamic responses generated by these considerations. Considering this, there are different methods to take account of uncertainties in this kind of problem such as perturbation method, Monte Carlo simulations or Polynomial Chaos Expansion. The perturbation and Neumann methods based respectively on the expansion of random quantities into Taylor series [START_REF] Nayfeh | Perturbation Methods[END_REF] and Neumann series [START_REF] Benaroya | Finite element methods in probabilistic structural analysis: A selective review[END_REF][START_REF] Yamazaki | Neumann expansion for stochastic finite element analysis[END_REF], provide acceptable results for small random fluctuations but are not adapted for solving a dynamic problem close to resonance frequency. The Monte Carlo simulations (MCS) is well adapted to include uncertainties in a deterministic finite element model by generating samples of the random parameter. Nevertheless MCS have a high computational cost due to the fact that a high number of samples is necessary to have convergence of this method. The Polynomial Chaos Expansion (PCE) [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] is a method based on the representation of the stochastic processes and variables in a set of orthogonal bases of random variables. This method has shown to be a successful approach to solve uncertainty quantification problems that is why it is used in the present study.

So the objective of the present study is to propose a methodology to undertake the effects of uncertainties on the non-linear dynamic response of a rotating machinery with multi-faults. The present paper is organized as follows: in a first part, a brief description of the rotor system is given. Then the different faults in rotating machine studied in this paper are presented. As previously explained, the novelty in the attempted work appears in the formulation of such rotor system model that includes the uncertainties in physical parameters. So the second part of the paper presents the methodology of the Polynomial Chaos Expansion (PCE) with the Harmonic Balance Method in order to estimate evolutions of the n× amplitudes and the global non-linear responses of a rotor system with multi-faults and uncertainties. Then, the expansions of the random quantities for all the models of faults of the rotor system are presented. Finally the efficiency of the method is demonstrated through various numerical simulations and a comparison of the results obtained via Polynomial Chaos Expansion (PCE) and via Monte Carlo Simulation (MCS).

Rotor system 2.1 General equation of a rotor system

The system under study is illustrated in Figure 1. The rotor consists of a shaft with two discs. The shaft is discretized into 10 Euler beam finite elements with four degrees of freedom at each node, with an asymmetrical section on the elements 3, 4 and 5 and a constant circular section on the other elements. Two bearing supports are added at each end of the rotor. All the values of the parameters are given in Table 1.

Figure 1: Rotor system

After assembling the shaft elements, the discs and the bearing supports, the equation of motion of the complete rotor system is given by:

Mẍ + (C + ωG) ẋ + (K s + K b )x = f (1)
where M and G are the mass and gyroscopic matrices of the system. K s and K b are respectively the stiffness matrices of the asymmetric shaft and the bearings. C is the damping matrix of the shaft. f are the forces applied on the system including gravitational forces, unbalance forces and forces due to faults. ẍ, ẋ and x are the acceleration, the velocity and the displacement of the degrees of freedom Table 1: Model parameters of the rotor system. ω is the rotating speed of the shaft. It should be noted that all these matrices are expressed in the fixed frame, their construction will be detailed in the following parts of the paper.

Faults in rotating machinery

In this paper it has been chosen to investigate the most current faults in a rotor system. Four kinds of faults are studied here: an unbalance rotor, an asymmetric shaft cross section, a bended shaft and a misalignment fault.

In this section the rotating and fixed frames are respectively noted R(X, Y, z) and R 0 (x, y, z), and are linked by the following basis change matrix defined as

P R→R 0 =   cos(ωt) -sin(ωt) 0 sin(ωt) cos(ωt) 0 0 0 1   (2)
Figure 2: Unbalance disc

Unbalance

The unbalance is one of the most common fault in a rotor problem. It appears when the center of mass of a shaft element or a disc is not on the axis of rotation of the system. So this fault is usually modeled by an eccentric mass on a disc or on the shaft, as illustrated in Figure 2.

For the degrees of freedom [v w θ ψ] T of a disc, the eccentric mass can be modeled by the following forces system

where m e and d e are respectively the mass unbalance and the eccentricity. ψ u defines the initial angular position in relation to z-axis.

Rotor asymmetry

The shaft cross-section of the elements 3, 4 and 5 is asymmetric. The geometry of the section is given in Figure 3.

Consequently the asymmetric section has different principal moments of inertia in rotating frame (X a , Y a , z) which are given by the following expressions

I X = Y 2 dS I Y = X 2 dS (4) (5)
After calculations, we obtain where α, R and h define the angle, the radius and the depth of the asymmetric section as illustrated in Figure 3.

I X = R 4 2 (α - sin(2α) 2 ) + √ 2hR -h 2 (R -h) 3 I Y = R 4 2 (α + sin(2α) 2 ) + (2hR -h 2 )
In the following part of this section, the stiffness and mass matrices due to this rotor asymmetry will be calculated by using expressions of the strain and kinetic energies.

Strain energy

The nodal displacement vector is defined as

δ f ix = [v 1 w 1 θ 1 ψ 1 v 2 w 2 θ 2 ψ 2 ] T in fixed frame or δ rot = [V 1 W 1 Θ 1 Ψ 1 V 2 W 2 Θ 2 Ψ 2 ]
T in rotating frame (see Figure 4).

The strain energy of the shaft element is defined as

U = 1 2 l 0 EI X (V ′′ ) 2 + EI Y (W ′′ ) 2 dz (8)
where V and W are the displacements in rotating frame. l is the length of the element. E is the Young modulus. Shear effects are neglected. By using the following relations V = cos(ωt + ψ a )v + sin(ωt + ψ a )w W = -sin(ωt + ψ a )v + cos(ωt + ψ a )w [START_REF] Oncescu | Investigation of the stability and steady state response of asymmetric rotors, using finite element formulation[END_REF] with v and w the displacements in fixed frame, the strain energy can be rewritten as follows

U = 1 2 l 0 EI m v ′′ 2 + w ′′ 2 dz + 1 2 l 0 EI d v ′′ 2 -w ′′ 2 cos(2ωt + 2ψ a )dz + l 0 EI d v ′′ w ′′ sin(2ωt + 2ψ a )dz (10)
where I m and I d are the mean and deviatoric moments of the area and are equal to

I m = I X +I Y 2 and I d = I X -I Y 2 .
This energy can be calculated by using classical shape functions [START_REF] Lalanne | Rotordynamics-Prediction in Engineering[END_REF] defined as

v w = N 1 v 1 + N 2 v 2 + N 3 ψ 1 + N 4 ψ 2 N 1 w 1 + N 2 w 2 -N 3 θ 1 -N 4 θ 2 (11)
where

N 1 = 1 -3z 2 l 2 + 2z 3 l 3 , N 2 = 3z 2 l 2 -2z 3 l 3 , N 3 = -z + 2z 2 l -z 3 l 2 and N 4 = z 2 l -z 3 l 2 .
After calculations, the strain energy can be expressed as follows

U = 1 2 δ T f ix K e δ f ix (12)
where the stiffness matrix K e is defined in fixed frame as

K e = K e 0 + K e c cos(2ωt + 2ψ a ) + K e s sin(2ωt + 2ψ a ) (13)
The expressions of K e 0 , K e c and K e s are given by

K e 0 = EI m l 3             12 0 0 -6l -12 0 0 -6l 12 6l 0 0 -12 6l 0 4l 2 0 0 -6l 2l 2 0 4l 2 6l 0 0 2l 2 12 0 0 6l 12 -6l 0 Sym. 4l 2 0 4l 2             (14)
K e c = EI d l 3             12 0 0 -6l -12 0 0 -6l -12 -6l 0 0 12 -6l 0 -4l 2 0 0 6l -2l 2 0 4l 2 6l 0 0 2l 2 12 0 0 6l -12 6l 0 Sym. -4l 2 0 4l 2             (15)
K e s = EI d l 3             0 12 6l 0 0 -12 6l 0 0 0 -6l -12 0 0 -6l 0 -4l 2 -6l 0 0 -2l 2 0 0 6l -2l 2 0 0 12 -6l 0 0 0 6l Sym. 0 -4l 2 0             ( 16
)
where l is the length of the shaft element.

Kinetic energy

The kinetic energy of the shaft element is defined as

T = 1 2 ρS l 0 v2 + ẇ2 dz + 1 2 ρ l 0 I X ω 2 X + I Y ω 2 Y + 2I m ω 2 Z dz (17)
where ω X , ω Y and ω are the angular velocities of the shaft element and can be approximated as [START_REF] Lalanne | Rotordynamics-Prediction in Engineering[END_REF]

       ω X = θ sin(ωt + ψ a ) -ψ sin(ωt + ψ a ) ω Y = ψ cos(ωt + ψ a ) + θ sin(ωt + ψ a ) ω Z = ω + ψθ (18)
After calculations, the kinetic energy can be rewritten as follows

T = 1 2 ρS l 0 v2 + ẇ2 dz + 1 2 ρI m l 0 θ2 + ψ2 dz + 1 2 ρI d l 0 θ2 + ψ2 cos(2ωt + 2ψ a )dz + 1 2 ρI d l 0 ψθ sin(2ωt + 2ψ a )dz + ρI m ω l 0 ψθdz (19)
By using the shape functions and after integration of the previous equation, it leads to express the gyroscopic matrix and the mass matrices.

The gyroscopic matrix is defined as

G e = ρI m 15l             0 -36 -3l 0 0 36 -3l 0 0 0 -3l -36 0 0 -3l 0 -4l 2 -3l 0 0 l 2 0 0 -3l -l 2 0 0 -36 3l 0 0 0 3l skew-sym. 0 -4l 2 0             (20)
The elementary mass matrix is defined in fixed frame as

M e = Mt e 0 + Mr e 0 + M e c cos(2ωt + 2ψ a ) + M e s sin(2ωt + 2ψ a ) (21)
The expression Mt e 0 , Mr e 0 , M e c and M e s are given by

Mt e 0 = ρSl 420             156 0 0 -22l 54 0 0 13l 156 22l 0 0 54 -13l 0 4l 2 0 0 13l -3l 2 0 4l 2 -13l 0 0 -3l 2 156 0 0 22l 156 -22l 0 Sym. 4l 2 0 4l 2             (22)
Mr e 0 = ρI m 30l             36 0 0 -3l -36 0 0 -3l 36 3l 0 0 -36 3l 0 4l 2 0 0 -3l -l 2 0 4l 2 3l 0 0 -l 2 36 0 0 3l 36 -3l 0 Sym. 4l 2 0 4l 2             (23)
M e c = ρI d 30l             36 0 0 -3l -36 0 0 -3l -36 -3l 0 0 36 -3l 0 -4l 2 0 0 3l l 2 0 4l 2 3l 0 0 -l 2 36 0 0 3l -36 3l 0 Sym. -4l 2 0 4l 2             (24)
M e s = ρI d 30l             0 -36 -3l 0 0 36 -3l 0 0 0 3l 36 0 0 3l 0 4l 2 3l 0 0 -l 2 0 0 -3l -l 2 0 0 -36 3l 0 0 0 -3l Sym. 0 4l 2 0             ( 25
)
where ρ is the density.

Rotor bow

In this section, the study of the rotor in the presence of an initial deformation of the shaft is proposed. This deformation that is illustrated in Figure 5, is defined in fixed frame as x 0r , so the equation of the dynamics problem is given by

M r ẍr + D r ẋr + K r (x r -x 0r ) = f r (26)
where M r , D r and K r are respectively the mass, the damping and the stiffness matrices expressed in the rotating frame. The subscript r denotes the fact that the quantity is expressed in the rotating frame. For the reader comprehension, the deformation

x 0r in rotating frame R b (X b , Y b , z) is defined as x 0r = [x 0 (1) r ... x 0 (k) r ... x 0 (n) r ] T (27)
where x 0

(k)
r are the displacements of the k th node in R b . Consequently the initial deformation can be modeled as a force defined by

f b r = K r x 0r (28)
In the present paper, it has been chosen to study a bowed rotor which has the geometric properties

x 0 (k) r =[V (k) 0 0 Φ (k) ] T =[R 0 sin( k -1 n -1 π) 0 0 R 0 π L cos( k -1 n -1 π)] T (29)
where R 0 is the maximum of amplitude of the bow. Finally the force f b r is expressed in the fixed frame by using the basis change matrices between the rotating frames R b (X b , Y b , z) and R 0 (x, y, z). So after calculations this force f b r can be rewritten as synchronous rotating forces defined in the fixed frame by

f b = f b c cos(ωt + Ψ b ) + f b s sin(ωt + Ψ b ) (30)

Misalignment

Parallel misalignment In this part, the modeling proposed by Lees in [START_REF] Lees | Misalignment in rigidly coupled rotors[END_REF] is used. The coupling of the two parallel rotors of relative displacement δ p , is composed of N bolts. On the first shaft the bolts are distributed on a circle of radius r, on the second rotor the bolts are distributed on a same circle with an offset δ p from the axis of rotation. This configuration is illustrated on Figure 6.

In the present study it will be considered that the two rotors are acting at the same rotation speed. So the torsional deflection is not taking into account. The relative displacement is only considered in the vertical direction. The initial configuration of the coupling composed of N = 3 bolts is defined in Figure 7.

Due to the geometry of the coupling the extremities of each bolt have different trajectories. Consequently during a rotation of the coupling system, the energy stored in each bolt will evolve. This energy can be calculated with the trajectories of the bolt extremities and depends on the stiffness k b of bolt in x-direction and the y-direction.

For the first disc, the trajectory ---→ O 1 M j 1 of the j th bolt (for j = 1, 2 or 3) is defined in rotating frame as where φ j = (j -1) 2π N is the angle between two consecutive bolts (N defines the number of bolts). By introducing the basis change and the degrees of freedom [v c w c θ c ψ c ] T of the shaft extremity, these trajectories can be rewritten as

---→ O 1 M j 1 =   -r sin(φ j ) r cos(φ j ) 0   R (31)
---→ OM j 1 =   v c w c 0   R 0 +   -r sin(ωt + φ j ) r cos(ωt + φ j ) rθ cos(ωt + φ j ) + rψ sin(ωt + φ j )   R 0 (32)
For the second disc, the trajectory of the j th bolt (for j = 1, 2 or 3) is defined in rotating frame as

---→ O 2 M j 2 =   -r sin(φ j ) δ p + r cos(φ j ) 0   R (33)
or in fixed frame as

---→ O 2 M j 2 = P R→R 0   -r sin(φ j ) δ p + r cos(φ j ) 0   R =   -r sin(ωt + φ j ) -δ p sin(ωt) r cos(ωt + φ j ) + δ p cos(ωt) 0   R 0 (34) and ---→ OM j 2 =   -r sin(ωt + φ j ) -δ p sin(ωt) r cos(ωt + φ j ) -δ p (1 -cos(ωt)) 0   R 0 (35)
Figure 8: Model of angular misalignment By using the previous equations, the strain energy can be expressed as

U m1 = 1 2 k b N j=1 [( ---→ OM j 1 - ---→ OM j 2 ). -→ x ] 2 + [( ---→ OM j 1 - ---→ OM j 2 ). -→ y ] 2 = 1 2 k b N j=1 (v c + δ p sin(ωt)) 2 + (w c + δ p (1 -cos(ωt))) 2 (36)
(37)

After summation on each bolt of the same properties, the total strain energy due to the parallel misalignment can be rewritten as

U m1 = 1 2 Nk b [(v c + δ p sin(ωt)) 2 + (w c + δ p (1 -cos(ωt))) 2 ] (38)
Finally the Lagrange's equation is applied to the degree of freedom [v c w c θ c ψ c ] T of the shaft extremity. Consequently after calculations, it results that the presence of a coupling with a parallel misalignment leads to add to the previous system: a stiffness component due to the coupling which is defined as

K m 1 = Nk b     1 1 0 0     (39)
and a synchronous force due to the misalignment

f m 1 = Nk b δ p [sin(ωt) 1 -cos(ωt) 0 0] T (40)
Angular misalignment The proposed model of angular misalignment is illustrated in Figure 8. It is composed of two discs with an angular fault α m and N = 3 bolts. The first bolt has a stiffness k + k ′ in z-direction and the others bolts have a stiffness k.

In the following part, an approach similar to these developed in the case of a parallel misalignment is proposed. The position ---→ OM j 1 of the left extremity of the j th bolt (with j = 1, 2 or 3) is identical to 32). Then, the second shaft have an angular fault α m in x-axis as illustrated in Figure 8. Consequently the position ---→ O 2 M j 2 can be defined as

j k j φ j 1 k + k ′ 0 2 k 2π/3 3 k 4π/3
---→ O 2 M j 2 = r   -sin(ωt + φ j ) cos(ωt + φ j ) cos(α m ) cos(ωt + φ j ) sin(α m )   R 0 (41)
Assuming that the angular α m is small, the following approximations cos(α m ) ≃ 1 and sin(α m ) ≃ α m can be done. So the position ---→ O 2 M j 2 can be rewritten as

---→ O 2 M j 2 = r   -sin(ωt + φ j ) cos(ωt + φ j ) α m cos(ωt + φ j )   R 0 (42)
Consequently the strain energy of one bolt can be rewritten as

U j m2 = 1 2 k j [( ---→ OM j 1 - ---→ OM j 2 ). -→ z ] 2 = 1 2 k j r 2 [(θ -α m ) cos(ωt + φ j ) + ψ sin(ωt + φ j )] 2 (43)
(44)

So the strain energy for all bolts is defined as

U m2 = N i=j U j (45)
which after calculations gives the following expression

U m2 = 1 4 (3k + k ′ )r 2 [(α m -θ) 2 + ψ 2 ] + 1 4 k ′ r 2 [((α m -θ) 2 + ψ 2 ) cos(2ωt) -2ψ(α m -θ) sin(2ωt)] (46)
Finally by using the Lagrange's equation, it is possible to show that an angular misalignment fault can be modeled with a parametric stiffness and a force component. The expression of the stiffness matrix is similar to the development for an asymmetric shaft. The final expression of this stiffness matrix K m 2 is defined as follows

K m 2 = 1 2 (3k + k ′ )r 2     0 0 1 1     + 1 2 k ′ r 2     0 0 cos(2ωt) sin(2ωt) sin(2ωt) -cos(2ωt)     (47)
The angular misalignment force is given by

f m 2 = - 1 2 α m r 2 [0 0 3k + k ′ + k ′ cos(2ωt) k ′ sin(2ωt)] T (48)

Final expression

Finally the equation of motion of the complete system, that has been previously defined in Equation [START_REF] Edwards | Fault diagnosis of rotating machinery[END_REF], is rewritten as

Mẍ + D ẋ + Kx = f (49) with M = M 0 + M c cos(2ωt) + M s sin(2ωt) K = K 0 + K c cos(2ωt) + K s sin(2ωt) D = C + ωG (50) (51) (52)
M, K and D are the global mass, stiffness and damping matrices. They are determinated by assembling the shaft elementary matrices M e , C e , G e , K e , K m 1 and K m 2 . f is defined by assembling the different sources of external forces such as unbalance, misalignment reaction force, bow and gravity. In this study, the damping is given as follows

C = α d K 0 + β d M 0 (53)
where α d and β d are the damping coefficients and C is the Rayleigh damping matrix. For the interested reader, it may be remained that the damping is overwhelmingly in the bearings in large machines on oil film bearings. In this case, more complex representations of damping need to be taken into account [START_REF] Bently | Fundamentals of Rotating Machinery Diagnostics[END_REF].

Non-linear dynamic response of the rotor system with uncertainties

The fault parameters can be partly unknown, consequently, to model them correctly, the choice has been to take these uncertainties in the model of each one developed in the previous section. Then, the response is uncertain and we chose to model it via the Polynomial Chaos Expansion. So, the aim of this section is to mix the Harmonic Balance Method and the Polynomial Chaos Expansion to obtain the non-linear response of a rotor with uncertain parameters. In a first part the Harmonic Balance Method which is one of the mathematical approaches to solve equations with nonlinear terms is presented. In a second part, the random character of the fault parameters will be investigated through stochastic model of each fault. Then a methodology will be proposed to take into account the uncertainties in this kind of problem by using a Polynomial Chaos Expansion with the Harmonic Balance Method.

Deterministic problem

Due to the presence of faults, the rotor system is subjected to periodic excitation. So the response x(t) and the excitation f(t) can be expressed as finite Fourier series of order m [19]

x(t) = m n=0 (A n cos(nωt) + B n sin(nωt)) (54)
and

f(t) = m n=0 (S n cos(nωt) + T n sin(nωt)) (55)
Consequently the amplitude of the n th response is defined as |A n + iB n | with i the complex number. By substituting the previous equations in the rotor system equations, it leads to solve a linear system of size (2m + 1)n dof where n dof is the number of degrees of freedom

LX = f (56)
with L, X and f are defined as

L = I 1 + I 2 + I 3 X = [A 0 A 1 B 1 ... A m B m ] T f = [S 0 S 1 T 1 ... S m T m ] T (57) (58) (59)
I 1 , I 2 and I 3 are respectively the contributions due to the constant terms, the parametric terms of the stiffness and the parametric terms of the mass.

In our study, some coupling terms appear due to the parametric stiffness (see Section 3.1.2). So an excitation of n× induces an evolution of the n -2 and n + 2 super harmonic responses. Consequently the truncation order of the non-linear response (54) has to be chosen in order to show this phenomenon. In this paper the m th order is equal to 4.

In the following sections, definitions of the terms I 1 , I 2 and I 3 are investigated.

Constant terms -Definition of I 1

Considering the constant terms of the left side of the equation of motion of the complete rotor system defined in Equation (49), it leads to the following ormulation

M 0 ẍ + D ẋ + K 0 x (60)
By using Equations (54-55) in Equation (60) and after calculations, the formulation can be written as the linear expression I 1 X with The non-constant terms of the left side of the equation of motion defined in Equation ( 49) can be rewritten such as follows

I 1 =      K 0 K 0 -ω 2 M 0 ωD -ωD K 0 -ω 2 M 0 . . . K 0 -m 2 ω 2 M 0 mωD -mωD K 0 -m 2 ω 2 M 0      (
[K c cos(2ωt) + K s sin(2ωt)]x [M c cos(2ωt) + M s sin(2ωt)]ẍ (62) (63)
The displacement x and the acceleration ẍ can be expressed as truncated Fourier series. By using the definition of x(t) given in Equation (54) in Equation (62), the parametric stiffness system can be rewritten as

m n=0 (K c A n + K s B n ) cos((n -2)ωt) + (K c A n -K s B n ) cos((n + 2)ωt) + (K c B n -K s A n ) sin((n -2)ωt) + (K c B n + K s A n ) sin((n + 2)ωt) = 0 (64)
Considering a truncation of order m, the previous system can be written as the linear system I 2 X = 0 where

I 2 = 1 2              Kc Ks Kc Ks Kc Ks Ks -Kc -Ks Kc 2Kc 2Ks Kc -Ks Ks Kc . . . Kc Ks -Ks Kc . . . Kc -Ks Ks Kc              (65)
Similarly, using the following expression of the acceleration,

ẍ(t) = -ω 2 m n=0 n 2 A n cos(nωt) + n 2 B n sin(nωt) (66)
one development on the parametric mass system can be done to Equation (63) which leads to linear system I 3 X = 0 where

I 3 = - ω 2 2               4Mc 4Ms Mc Ms 9Mc 9Ms Ms -Mc -9Ms 9Mc 0 0 Mc -Ms Ms Mc . . . m 2 Mc m 2 Ms -m 2 Ms m 2 Mc . . . (m-2) 2 Mc -(m-2) 2 Ms (m-2) 2 Ms (m-2) 2 Mc               (67)

Uncertainties on n× Harmonic response via the Polynomial Chaos Expansion

Since the material and geometrical properties of the rotor and the excitation force are taken as uncertain, we need to model the response in a stochastic way. Material and geometrical properties are called M(τ ), D(τ ) and K(τ ), the excitation force is written as f (τ ), then Equation (56) can then be rewritten

L(τ ) X(τ ) = f (τ ) (68)
where argument τ denotes the random character and

L(τ ) = Ĩ1 (τ ) + Ĩ2 (τ ) + Ĩ3 (τ ) (69)
and Ĩi (τ ), i = 1 to 3, depend on the material and geometrical properties of the rotor. Each uncertain quantity (material, geometrical or fault parameter of the rotor) is modeled by implementing the Karhunen-Loeve expansion in the Galerkin formulation of the finite element method [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF], it gives :

ã = a + L l=1 ξ l a l (70)
for one uncertain quantity ã (scalar or matrix), where {ξ 1 , ...ξ L } is a set of orthonormal random variables, the notation a refers to the mean of the quantity ã and a l is its lth Karhunen-Loeve expansion term. Integrating these models into Ĩi (τ ), quantities L(τ ) and f(τ ) can then be expressed on a basis of orthogonal random variables, the Polynomial Chaos basis (details will be given in the following subsection) such that

L(τ ) = ∞ i=0 L i Ψ i (ξ(τ )), f (τ ) = ∞ j=0 f j Ψ j (ξ(τ )) (71) (72)
where Ψ j (ξ(τ )) are multidimensional Hermite polynomials and where L i (i = 0 to ∞) are obtained by identification with I 1 (τ ), I 2 (τ ) and I 3 (τ ).

Finally, the response has to be also expanded on the Polynomial Chaos basis such as

X(τ ) = ∞ j=0 X j Ψ j (ξ(τ )) (73)
where X j is the unknown deterministic j th vector associated with Ψ j (ξ(τ )) and ξ = {ξ r } [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Loève | Probability Theory[END_REF]. So the system to be solved, when expanded on the polynomial chaos basis, becomes

∞ j=0 L(τ )X j Ψ j (ξ(τ )) = f(τ ) (74)
that can be rewritten as

∞ i=0 L i Ψ i (ξ(τ )) ∞ j=0 X j Ψ j (ξ(τ )) = ∞ j=0 f j Ψ j (ξ(τ )) (75)
Projecting the system of equations on the subspace spanned by

{Ψ k } ∞ k=0 yields ∞ i=0 ∞ j=0 E{Ψ i Ψ j Ψ k }L i X j = f k E{Ψ 2 k } k = 0, 1, ..., ∞ (76)
where E{} denotes the operation of mathematical expectation.

It should be noted that, in practice, the expansion can be truncated after the P th term. P is the total number of polynomial chaoses used in the expansion excluding the 0 th order term and can be determined by P + 1 = (p+r)! p!r! in which p is the order of homogeneous chaos used and r is the number of random variables. Linear system Equation (76) to be solved in practice is then given by

L X = f (77)
with components

Ljk = P i=0 E{Ψ i Ψ j Ψ k }L i j, k = 0, ..., P (78)
and where X = [X 0 ...

X P ] T , f = [f 0 E{Ψ 2 0 } ... f P E{Ψ 2 P }] T (79)
It should be noted that coefficients E{Ψ i Ψ j Ψ k } and E{Ψ 2 k } only have to be calculated once.

Uncertainties in multi-faults rotor system

It is difficult to model accurately faults in a deterministic way. We then chose to model them in a stochastic way. Here, we present one model for each treated fault.

Unbalance

We consider the parameters m e and ψ u as stochastic. From Equation 70, they can be developed as following

me (τ ) = m e (1 + δ m ξ 1 ) ψu (τ ) = σ ψu ξ 2 (80)
(81) where ξ 1 and ξ 2 are Gaussian random variables. Unbalance mass me is defined by its mean m e and its variation coefficient δ m ; angular position of the force ψu has mean 0 and standard deviation σ ψu . By substituting Equations (80-81) in Equation [START_REF] Xu | Vibration analysis of a motor-flexible Coupling-Rotor system subject to misalignment and unbalance, part II: experimental validation[END_REF]

cos(σ ψu ξ 2 ) = m l=0 (-1) l (σ ψu ξ 2 ) 2l (2l)! sin(σ ψu ξ 2 ) = m l=0 (-1) l (σ ψu ξ 2 ) 2l+1 (2l + 1)! (85) (86)
the unbalance force can be rewritten such as follows

fd = 1 k=0 m l=0 fd kl ξ k 1 ξ l 2 ( 87
)
where fd kl are the unbalance force components in a basis (ξ k 1 , ξ l 2 ) and are defined as

fd kl = (-1) l (δ m ξ 1 ) k (σ ψu ξ 2 ) 2l (2l!) ([cos(ωt) sin(ωt) 0 0] T + σ ψu ξ 2 2l + 1 [-sin(ωt) cos(ωt) 0 0] T ) (88)
Finally, the uncertain loading can be expanded on the polynomial chaos basis as

fd = R j=0 fd j Ψ j (ξ 1 , ξ 2 ) ( 89
)
where R is the number of the polynomial chaoses and parameters fd j are deduced from fd kl after one identification procedure between Ψ j (ξ 1 , ξ 2 ) basis and (ξ k 1 , ξ l 2 ) basis (see Table 3 for the case of Two-Dimensional Polynomial Chaoses). i j

ξ i 1 ξ j 2 0 0 Ψ 0 1 0 Ψ 1 0 1 Ψ 2 2 0 Ψ 3 + Ψ 0 1 1 Ψ 4 0 2 Ψ 5 + Ψ 0 3 0 Ψ 6 + 3Ψ 1 2 1 Ψ 7 + Ψ 2 1 2 Ψ 8 + Ψ 1 0 3 Ψ 9 + 3Ψ 2
Table 3: Expression of ξ i 1 ξ j 2 in polynomial chaos basis for order 2 of the expansion

Rotor asymmetry

Material parameters and geometrical properties can be often subject to uncertainties. Then, we can model Young modulus and the stochastic depth of the asymmetry as random.

From Equation (70), we can write :

Ẽ(τ ) = E(1 + δ E ξ 3 ) (90)
where ξ 3 is one Gaussian random variable, E and δ E is the mean and the variation coefficient of the Young modulus respectively. The uncertain geometrical property considered here is one characteristic of the asymmetric shaft cross section : the depth of the asymmetry, as illustrated in Figure 9.

Then, following again Equation (70), we have :

h(τ ) = h(1 + δ h ξ 4 ) (91)
where ξ 4 is one random Gaussian variable, h and δ h the mean and the variation coefficient of the depth. Let us denote

h ′ (τ ) = hδ h ξ 4 (92)
since we can write Equation (91) as h(τ ) = h + h ′ . h and h ′ are presented in Figure 9 for one ξ 4 given. Consequently, the moments of inertia are stochastic quantities and are defined by

ĨX (τ ) = I X + bh ′ 3 6 + 2bh ′ (R -h - h ′ 2 ), ĨY (τ ) = I Y + h ′ b 3 6 (93) with b = 2 √ 2hR -h 2 (94)
Figure 9: Uncertainty on the shaft section I X and I Y are deterministic moments defined by Equations (6-7). Substituting Equation (92) in Equation ( 93) and after developing these equations, we have

ĨX (τ ) = I X + 2bh(R -h)δ h ξ 4 -bh 2 δ h 2 ξ 2 4 + bh 3 6 δ h 3 ξ 3 4 , ĨY (τ ) = I Y + hδ h b 3 6 ξ 4 (95)
which can also be expanded as

ĨX (τ ) = 3 i=0 ĨX i ξ i 4 , ĨY (τ ) = 1 i=0 ĨY i ξ i 4 ( 96
)
After one identification procedure between Equation (95) and Equation (96), we obtain components ĨX i and I Y i given by

ĨX 0 = I X , ĨX 1 = 2bh(R -h)δ h , ĨX 2 = -bh 2 δ h 2 , ĨX 3 = bh 3 6 δ h 3 , ĨY 0 = I Y , ĨY 1 = hδ h b 3 6 . (97)
(98)

It yields the stochastic expressions of the mean and deviatoric moments of inertia

Ĩm (τ ) = ĨX (τ ) + ĨY (τ ) = 3 i=0 ( ĨX i + ĨY i )ξ i 4 = 3 i=0 Ĩm i ξ i 4 Ĩd (τ ) = ĨX (τ ) -ĨY (τ ) = 3 i=0 ( ĨX i -ĨY i )ξ i 4 = 3 i=0 Ĩd i ξ i 4 (99) (100)
Introducing stochastic expressions Ẽ, Ĩm and Ĩd (Equations ( 90), ( 99), (100)) in the deterministic expressions of the elementary stiffness matrices given by Equations (14-16) and also valid for stochastic quantities, the stochastic elementary stiffness matrices are given by

Ke 0 = 3 i=0 1 j=0 Ke 0 ij ξ j 3 ξ i 4 , Ke c = 3 i=0 1 j=0 Ke c ij ξ j 3 ξ i 4 , Ke s = 3 i=0 1 j=0 Ke s ij ξ j 3 ξ i 4 (101) with Ke 0 ij = Eδ j E Ĩm i K e 0 EI m , Ke c ij = Eδ j E Ĩd i K e c EI d , Ke s ij = Eδ j E Ĩd i K e s EI d (102)
Then, after assembling the elementary matrices, the stiffness matrices of the structure are given by

K0 = 3 i=0 1 j=0 K0 ij ξ j 3 ξ i 4 , Kc = 3 i=0 1 j=0 Kc ij ξ j 3 ξ i 4 , Ks = 3 i=0 1 j=0 Ks ij ξ j 3 ξ i 4 ( 103
)
where the components K0 ij , Kc ij and Ks ij are constructed by assembling the matrices Ke 0 ij , Ke c ij and Ke s ij . Rewriting the developments of the stiffness matrices (Equation ( 103)) onto a polynomial chaos basis yields

K0 = n j=0 K0 j Ψ j (ξ 3 , ξ 4 ), Kc = n j=0 Kc j Ψ j (ξ 3 , ξ 4 ), Ks = n j=0 Ks j Ψ j (ξ 3 , ξ 4 ) (104)
where n is the number of Hermite polynomials. K0 j , Kc j and Ks j are the rearrangement of K0 ij , Kc ij and Ks ij on the polynomial chaos basis and are obtain after one identification procedure between the two bases. Similarly, we obtain from the deterministic expressions of the elementary mass matrices given by Equations (23-25) and also valid in the stochastic domain, the following expressions of the stochastic mass matrices

Mr e 0 = 3 i=0 Mr e 0 i ξ i 4 , Me c = 3 i=0 Me c i ξ i 4 , Me s = 3 i=0 Me s i ξ i 4 ( 105
)
where coefficients of the expansions are given by

Mr e 0 i = Ĩm i I m Mr e 0 , Me c i = Ĩd i I d M e c , Me s i = Ĩd i I d M e s ( 106
)
After assembly and identification procedures, the mass matrices can completely be defined by

M0 = n j=0 M0 j Ψ j (ξ 3 , ξ 4 ) , Mc = n j=0 Mc j Ψ j (ξ 3 , ξ 4 ) , Ms = n j=0 Ms j Ψ j (ξ 3 , ξ 4 ) (107)
Finally, the linear system to solve is given by Equation (68) with Equation (69). However, Equation (69) can be rewritten as

L(τ ) = n j=0 Lj Ψ j (τ ) (108)
where Lj = Ĩ1 j + Ĩ2 j + Ĩ3 j . Components Ĩi j , i = 1 to 3, can be constructed by one identification procedure between Equation (69) with Equations ( 61), ( 65), (67) also valid in the stochastic domain and the following expansions of Ĩi , i = 1 to 3, on the Polynomial Chaos basis :

Ĩ1 (τ ) = n j=0 Ĩ1 j Ψ j (ξ 3 , ξ 4 ) , Ĩ2 (τ ) = n j=0 Ĩ2 j Ψ j (ξ 3 , ξ 4 ) , Ĩ3 (τ ) = n j=0 Ĩ3 j Ψ j (ξ 3 , ξ 4 ) (109)

Rotor bow

In this model of bowed shaft, the parameters which seem the most subjected to uncertainties are the amplitude and the phase of the bow. Considering the uncertainties on the phase of the bow, a mathematical development similar to that exposed for the unbalance force can be used. Then, we chose to investigate the random character of the amplitude parameter R0 alone defined by

R0 (τ ) = σ R ξ 5 (110)
where ξ 5 is one random Gaussian variable and σ R is the standard deviation of the amplitude of the bow. We know that the deterministic equivalent force f b is expressed as one linear function of R 0 from Equation (30) with Equations ( 28) and (29). Therefore, fb (τ ) verifies the same relation with R0 (τ ) in the stochastic domain :

f b = α R 0 , fb (τ ) = α R0 (τ ) (111) (112)
Substituting Equation (110) into Equation (112) and using α from Equation (111) gives

fb (τ ) = fb 1 ξ 5 with fb 1 = σ R R 0 f b (113)

Misalignment

Parallel misalignment

The uncertain parameter in this modeling is the relative displacement δp which can be defined as

δp (τ ) = σ δp ξ 6 ( 114
)
where ξ 6 is a random Gaussian variable and σ δp is the standard deviation of the amplitude of parallel displacement. Substituting Equation (114) in Equation (40), the synchronous force due to the parallel misalignment can be defined by

fm 1 (τ ) = fm 1 1 ξ 6 (115)
where

fm 1 1 = Nk b σ δp [sin(ωt) 1 -cos(ωt) 0 0] T (116)
Angular misalignment Uncertainties can also occur in the determination of the angular fault α m which needs to be modeled in a stochastic way. Let ξ 7 be one random Gaussian variable, the angular fault αm is defined by

αm (τ ) = σ αm ξ 7 ( 117
)
where σ αm is the standard deviation of the angular fault. Finally, substituting Equation (117) in Equation ( 48), the moment due to the coupling misalignment can be expressed as

fm 2 (τ ) = fm 2 1 ξ 7 ( 118
)
where

fm 2 1 = - 1 2 σ αm r 2 [0 0 3k + k ′ + k ′ cos(2ωt) k ′ sin(2ωt)] T (119)

Numerical results

In this section, the effects of uncertainties on the dynamic response of a rotor system with faults are quantified via the Polynomial Chaos Expansion method. The random dynamic response of the rotor system is defined by its mean value and its variance, besides the envelope of the random response is presented. It is constructed by calculating the maximum and the minimum of all the responses computed by the PCE approach from samples generated following the random law of each input parameter.

The basic deterministic parameters used in all cases are given in Table 1. For the reader comprehension, the initial deterministic rotor system includes multi-faults such as unbalance, asymmetric shaft, bow, parallel and angular misalignments. In addition to this deterministic formulation, uncertainties on the different sources of fault are investigated (see Table 4 for the values of parameters and the order of the PCE taken). For the reference solution computed from the Monte Carlo Simulations the number of samples is taken equal to 1000. The PCE envelope is constructed with 20000 samples. All the results detailed in this section are presented for a rotor frequency evolution from 5 Hz to 275 Hz with 400 points of discretization.

Stiffness uncertainties

In this part, we consider the uncertainty present on the stiffness properties of the rotor. The variation coefficient on the Young modulus is taken as δ E = 5%. Computations are done with the MCS and the PCE approaches. Figure 10 presents the mean and variance of the Frequency Response Functions (FRF) of the n× amplitudes (for n=1,2,3 and 4) and the global non-linear response (at node 6 in x-direction) for case 1. It is observed that the mean and variance of the n-FRF, obtained by applying the PCE, are very similar to those of the MCS. Figure 11 presents all the FRF samples obtained by MCS compared to the envelope built by PCE. It can be noted that the contributions of 3× and 4× super-harmonic components are very small (respectively 10 -7 m and 10 -9 m of maximum

δ E δ m σ ψu δ h σ R σ δp σ αm Chaos order Case 1 5% - - - - - - 2 Case 2 - 5% 0.06 rad - - - - 2 Case 3 - - - 2% - - - 2 Case 4 - - - - 10 -5 m - - 2 Case 5 - - - - - 10 -4 m - 2 Case 6 - - - - - - 5
× 10 -4 rad 2 Case 7 2% 2% 0.03 rad 1% 10 -5 m --2 Case 8 5% 5% 0.06 rad 2% 10 -5 m 10 -4 m 5 × 10 -4 rad 2 Table 4: Sets of parameters amplitude), so the global response is mainly due to the 1× harmonic and the 2× super-harmonic components (respectively 10 -2 m and 10 -4 m of maximum amplitude) as illustrated in Figure 11(f). For this example the computational cost of the two methods have been compared and it shows that the Polynomial Chaos Expansion is 40 times faster that the Monte Carlo Simulations.

In this study, the PCE calculations are done with a chaos order equal to 2. It is observed that the difference between both PCE and MCS methods for this order is not significant either for the mean or the variance. However, a better accuracy on the model especially around critical speeds, will be obtained with a higher order than 2, but it would increase significantly the CPU time. Therefore, we consider the PCE approach with order 2 sufficient and it will be used for all computations.

Uncertainties in the unbalance force

The effects of uncertainties in the unbalance force are now studied. Variation of 5% on the mass unbalance (δ m ) and variation of 0.06 rad on the angular position (ψ u ) are considered. Figure 12 shows the mean and the variance of the MCS and the PCE results for the n× amplitudes of the rotor system (with n=1,2,3 and 4) and the global non-linear response. As previously observed for case 1, the results obtained by the PCE formulation are similar to those of the MCS. Due to the fact that comparisons of the mean and the variance for the MCS and the PCE processes will give the same results for the cases 3, 4, 5 and 6 (i.e. the mean and the variance for the MCS and PCE results will be similar), these representations will not be given in the following part of the paper in order to avoid repetitive informations.

Then, Figure 13 shows the Frequency Response Functions for the global non-linear response and the n× super-harmonic component (for n = 1, . . . , 4) for both MCS and PCE methods. A very good agreement is obtained between the results from the MCS and the PCE: all the MCS results are bounded by the lower and upper envelopes that have been obtained by applying the PCE process. As expected, these random parameters influence mainly the 1× response. However, due to the coupling terms of asymmetry, evolution of the 3× super harmonic components is also observed. Indeed, coupling terms between the 1× and 3× contributions are presents in the expressions of I 2 and I 3 (see Equations ( 65) and ( 67)). So a variation on the terms of unbalance causes a variation on the terms of order 3 (see Figure 13 where it is shown that variations of the n× amplitudes are only present for the first and third orders of the response. For the static component, the second and fourth orders, the variations of the response amplitudes are negligible (equal to zero if we consider the numerical errors inherent in calculations). Figure 13(f) indicates also that small variations of the amplitudes around the critical speeds of the 1× harmonic component are present in the global non-linear response of the rotor system. The range of amplitude is 2.2 × 10 -3 m to 2.7 × 10 -3 m, 5.4 × 10 -4 m to 1.6 × 10 -3 m and 3.2 × 10 -3 m to 4.0 × 10 -3 m for the first critical speed (28.7 Hz), the second critical speed (99.0 Hz) and the third critical speed (241.8 Hz) respectively.

Uncertainties in the asymmetric section properties

In this subsection, we are interested in the effect of variations on the properties of the shaft crosssection that is through h the depth of the asymmetry and δ h equal to 3% is taken. Figure 14 shows the n-FRF for both the Monte Carlo simulations and the lower and upper envelopes that have been calculated by the PCE method. Figures 14(a), 14(c), 14(d) and 14(e) show variations of the 0×, 2×, 3× and 4× harmonic components. Figure 14(b) indicates that variations of the 1× harmonic component are very small. Consequently, Figure 14(e) shows that variations on the global non-linear response are small too due to the fact that 1× harmonic component represents the predominant vibration of the global non-linear response. Indeed, it is interesting to note here that a variation in the depth of the asymmetry causes significant changes to the contributions of the matrix I 2 (see Equation ( 65)), while it has little impact on the responses of order 1 because its contribution in the expression of the matrix I 1 (see Equation (61)) is not predominant. We recall here that the presence and the variations observed for the 0×, 2×, 3× and 4× harmonics are due to the relations between the different contributions (i.e. static terms and the terms of orders 2, 3 and 4 given in the analytical expressions of Equation (65)).

Uncertainties in bow characteristics

Now, we present the impact of the presence of uncertainties on the parameters characterizing the extended bow. A variation of the bow amplitude R 0 is considered with a standard deviation σ R of 10 -5 m.

Figure 15 shows the n-FRF for the MCS and the PCE method (with n=1,2,3 and 4). Figures 15(a), 15(c) and 15(e) indicate no variability on the 0×, 2× and 4× super harmonic components whereas Figures 15(b) and 15(d) show a high variability on the 1× and 3× super harmonic components due to uncertainties in bow characteristics. For the reader comprehension, it is recalled that deterministic multi-faults such as unbalance, asymmetric shaft, bow, parallel and angular misalignments are present in the initial rotor system. Indeed, coupling between the 1× and 3× harmonic components exists, as indicated in the expressions of I 2 and I 3 (see Equations ( 65) and (67)). So a variation of the bow characteristics that generates a variation of the 1× (as previously indicated in Equation ( 40)) causes also a variation on the terms of order 3. For example, at the first 1× critical speed (28.7 Hz), it is observed variation from 9.3 × 10 -5 m to 4.5 × 10 -3 m for the 1× harmonic component and variation from 2.4 × 10 -11 m to 1.2 × 10 -9 m for the 3× harmonic component. Interpretations of these results are similar to the explanations given for unbalance force due to the fact that the forces that model the bow are acting as 1× component and consequently affect the 1× and 3× non-linear responses via the coupling terms between odd terms of Equation (65). Moreover, a high variation of the antiresonance position is observed (as illustrated in Figure 15(b)). It illustrates the fact that the addition of an initial deformation modifies the dynamic response and particularly the position of the minimum of amplitude on the shaft, especially when the initial bow is opposed to the deformation due to unbalance. This phenomenon is illustrated on the first critical speed on the Monte Carlo Simulations : for higher amplitudes than deterministic response, there is no antiresonance (same direction for bow and dynamic deformations) while for lower amplitudes, an antiresonance appears (opposite direction for bow and dynamic deformations).

Uncertainties in misalignment

In this part two cases are studied: the first case deals with an uncertainty on the parallel misalignment fault δ p for a standard deviation σ δp of 10 -4 m and the second case deals with an angular misalignment α m for a standard deviation σ αm of 5 × 10 -4 rad.

Figure 16 gives the n× super harmonic components (with n=1,2,3 and 4) and the global non-linear response by using the Monte-Carlo simulations. The upper and lower envelopes that have been calculated by applying the PCE method are also indicated. Firstly, a small variability of the 1× harmonic component is observed, due to the excitation of the first order given by the parallel misalignment (see Equation 40). Secondly, a very small variation of the Frequency Response Functions of order 3 appears. As explained in the previous paragraphs, the origin of this phenomenon is due to the coupling terms between the first and third orders in the matrices I 2 and I 3 , as seen in Equations ( 65) and (67).

Then, Figure 17 shows the n× harmonic components (n=1,2,3 and 4) and the global non-linear response for a parallel misalignment via Monte-Carlo simulations for an angular misalignment. The lower and upper envelopes via the PCE method are also given for all the Frequency Response Functions of the n× super-harmonic components (n = 1, . . . , 4) and the global non-linear response. Here again, the PCE approach is validated by the Monte-Carlo simulations. It is clearly shown that the 2× super harmonic component has a high variability due to the fact that the angular misalignment forces are acting as 2× component (see Equation 48). Moreover, variabilities of the 0× and 4× harmonic components can be observed due to the coupling terms in Equations ( 65) and (67) between the static components, the second and fourth orders.

Finally, differences between the results obtained for the variability of the parallel misalignment and angular misalignment fault are clearly observed (see Figures 16(a),(c) and (e) and Figures 17(a),(c) and (e)). It can be concluded that the variability of an angular misalignment has important effects on the contributions of the 0×, 2×, and 4× super-harmonic components, whereas an parallel misalignment affects only slightly the 1× responses of the rotor system.

Uncertainties in multi-faults

The main objective of this part of the paper is to demonstrate the efficiency of the new procedure (i.e. the combination of the Harmonic Balance Method and the Stochastic Finite Element Method using the Polynomial Chaos Expansion) when uncertain quantities come from all the parameters simultaneously. Case 7 deals with variations of mechanical properties of the shaft (with a variation coefficient of the Young modulus δ E = 2%), an unbalance excitation (with a variation coefficient mass unbalance δ m = 2% and a standard deviation of the angular position σ ψu = 0.03 rad) and a bow amplitude (with a standard deviation of the bow σ R = 10 -5 m). For the case 8, parallel and angular misalignments of 10 -4 m and 5.10 -4 rad respectively are added to the previous case. All the uncertainties on the different sources of fault are given in Table 4.

The mean and the variance of the n× harmonic components (n=1,2,3 and 4) and the global nonlinear response for both methods MCS and PCE are presented on Figure 18 for case 7 and on Figure 20 for case 8. Figures 19 and21 show the n-FRF, the global non-linear response and the upper and lower envelopes via the PCE method, respectively for case 7 and case 8. The results of the global non-linear amplitudes and the n× harmonic components via the Monte-Carlo simulations are also given in gray lines in Figures 19 and21. These results given for variations on material and geometrical parameters, excitations and faults characteristics show that the mean and the variance of all the responses as well as the upper and lower PCE envelopes agree very well with the Monte Carlo samples. So it proves the efficiency of the PCE method to evaluate the non-linear response of a rotor system in presence of multi-faults such as unbalance, misalignment, bow and asymmetric shaft. For these two examples, the computational cost of the Polynomial Chaos Expansion increases significantly due to the fact there is several random variables. So it increases the size of the Polynomial Chaos basis and consequently the size of the linear system to solve Equation (77). For cases 7 and 8, the chaos order is taken equal to 2 for all random quantities, consequently the Polynomial Chaos basis size is equal to 21 and 36 respectively.

Figures 19(b,f) and 21(b,f) show that for these two cases the global non-linear displacement is mainly due to the first harmonic contribution. We recall that the appearances of the various components (i.e. evolutions of the mean and variance for the 1×, 2×, 3× and 4× components) are due to the presence of all the faults in the system, which creates excitations of orders 1 and 2 and the coupling terms in the expressions given in Equations 65 and 67.

By comparing Figures 19 and21, it appears that differences between these two cases are the more significant for 0×, 2× and 4× components, whereas the 1× and 3× responses seem to be very similar. It is due to the fact that uncertainties on the angular misalignment are added for case 8: as previously seen, an angular misalignment can affect the 0×, 2× and 4× harmonic components via the expressions of matrices I 2 and I 3 . For the reader comprehension, it can be noted that, a variability of n th order induces a variability of the (n -2) th th and (n + 2) th orders due to expressions of I 2 and I 3 in Equations ( 65) and (67).

Orbits

In this part of the paper, results are investigated through the orbits of the shaft in order to understand the dynamic behavior of the shaft with the observation of simple, double or complex loops for different rotor frequencies. Orbits are constructed by using displacements in x-direction and y-direction at the node under study.

Figure 22 shows the evolution of orbits at node 6 for various rotation speeds with a set of random parameters defined in case 7. It can be noted that the orbits are presented without the static deflection of the rotor in order to show the complexity of the non-linear response that generally occurs at low amplitudes. Figures 22(b), 22(d) and 22(f) show the orbit of the shaft for the first, second and third critical speeds respectively (at 28 Hz, 99 Hz and 243 Hz). Figures 22(a) and 22(c) show orbits at the first and second 2× critical speeds respectively (at 14 Hz and 45 Hz). For the 1× critical speeds simple loops are observed. For the 2× critical speeds double loops are observed, so it highlights the presence of 2× harmonic component. These results are in perfect agreement with the previous evolutions of the n-FRF illustrate on Figure 21.

Conclusion

This paper describes the use of the Polynomial Chaos Expansion approach coupled with the Harmonic Balance Method in a rotordynamics problem. The efficiency and robustness of the proposed method has been tested through various numerical results. The stochastic results are validated by comparison with the results obtained by Monte Carlo simulations. These results show that variations of faults such as unbalance, bow and parallel misalignment may affect all the harmonic component of the rotor system, and variations of faults such as asymmetric section or angular misalignment can affect all the harmonic components.

2 . 1 × 3

 213 10 11 Nm 2 Shear modulus G 8.0 × 10 10 Nm 2 Poisson ratio ν 0.Density ρ 7800 kg m -3 Mass unbalance 0.05 g Eccentricity of the mass unbalance 0.01 m Damping coefficient α d 2.5 × 10 -5 Damping coefficient β d 2.8 Depth of the asymmetry h 0.002 m Phase of the asymmetry ψ a 0 • Phase of the bow ψ b 0 • Axial stiffness of bolt k 5.2 × 10 3 Nm -1 Transversal stiffness of bolt k b 13 × 10 3 Nm -1

  f d = [m e d e ω 2 cos(ωt + ψ u ) m e d e ω 2 sin(ωt + ψ u ) 0 0] T

Figure 3 :Figure 4 :

 34 Figure 3: Asymmetric shaft cross-section

Figure 5 :

 5 Figure 5: Coordinate system of the bow

Figure 6 :

 6 Figure 6: Model of parallel misalignment

Figure 7 :

 7 Figure 7: Initial positions of bolts

61) 3 . 1 . 2

 312 Parametric terms -Definition of I 2 and I 3

Figure 10 :

 10 Figure 10: Frequency Response Functions (Case 1); Mean and Variance of the FRF with PCE (red solid line) and MCS (black dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 11 :

 11 Figure 11: Frequency Response Functions (Case 1); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 12 :

 12 Figure 12: Frequency Response Functions (Case 2); Mean and Variance of the FRF with PCE (red solid line) and MCS (black dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 13 :

 13 Figure 13: Frequency Response Functions (Case 2); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 14 :

 14 Figure 14: Frequency Response Functions (Case 3); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 15 :

 15 Figure 15: Frequency Response Functions (Case 4); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 16 :

 16 Figure 16: Frequency Response Functions (Case 5); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 17 :

 17 Figure 17: Frequency Response Functions (Case 6); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 18 :

 18 Figure 18: Frequency Response Functions (Case 7); Mean and Variance of the FRF with PCE (red solid line) and MCS (black dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 19 :

 19 Figure 19: Frequency Response Functions (Case 7); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 20 :

 20 Figure 20: Frequency Response Functions (Case 8); Mean and Variance of the FRF with PCE (red solid line) and MCS (black dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 21 :

 21 Figure 21: Frequency Response Functions (Case 8); Lower and upper envelopes (red solid line); MCS (gray dotted line); Deterministic response (black dotted dashed line): (a) static component; (b) 1× harmonic component; (c) 2× super-harmonic component; (d) 3× super-harmonic component; (e) 4× super-harmonic component; (f) global non-linear response.

Figure 22 :

 22 Figure 22: Orbits; PCE (green); MCS (black dashed line); Deterministic response (blue solid line): (a) f = 14 Hz ; (b) f = 28 Hz ; (c) f = 45 Hz ; (d) f = 99 Hz ; (e) f = 135 Hz ; (f) f = 243 Hz .

Table 2 :

 2 Coupling parameters the development given in Equation (

  , also valid in the stochastic domain, unbalance force fd can be rewritten as fd = m e d e ω 2 (1 + δ m ξ 1 )[cos(ωt + σ ψu ξ 2 ) sin(ωt + σ ψu ξ 2 ) 0 0] T

		(82)
	Using
	cos(ωt + σ ψu ξ 2 ) = cos(ωt) cos(σ ψu ξ 2 ) -sin(ωt) sin(σ ψu ξ 2 ) sin(ωt + σ ψu ξ 2 ) = sin(ωt) cos(σ ψu ξ 2 ) + cos(ωt) sin(σ ψu ξ 2 )	(83) (84)
	and assuming that

(R -h)

(6)[START_REF] Penny | Simplified modelling of rotor cracks[END_REF]

Acknowledgments

The first author gratefully acknowledges the French Education Ministry which supports this research. The authors gratefully acknowledge the financial support of the French National Research Agency through the Young Researcher program ANR-07-JCJC-0059-01-CSD2.


Jean Krivine

A verification technique for reversible process algebra

A verification method for distributed systems based on decoupling forward and backward behaviour is proposed. This method uses an event structure based algorithm that, given a CCS process, constructs its causal compression relative to a choice of observable actions. Verifying the original process equipped with distributed backtracking on non-observable actions, is equivalent to verifying its relative compression which in general is much smaller. The method compares well with direct bisimulation based methods. Benchmarks for the classic dining philosophers problem show that causal compression is rather efficient both time-and space-wise. State of the art verification tools can successfully handle more than 15 agents, whereas they can handle no more than 5 following the traditional direct method; an altogether spectacular improvement, since in this example the specification size is exponential in the number of agents.

Introduction

Backtracking is commonplace in transactional systems where different components, such as processes accessing a distributed database, need to acquire a resource simultaneously. To ensure unconditional correctness of the overall execution of the transaction, one usually provides a code that incorporates explicit escapes from those cases where a global consensus cannot be met. Such an upfront method generates a large and unstructured state space, which often means verification based on proving that the code is bisimilar to a reference specification becomes unfeasible. Based on earlier work, we propose here an indirect verification method, and show on an example that it can handle larger specifications. The idea is to break down the distributed implementation of a given reference specification in two steps. First, one writes down a code which is only required to meet a weaker condition of causal or forward correctness relative to the specification. This condition is parameterized by a choice of observable actions corresponding to the actions of the specification. Second, the obtained code is equipped with a generic form of distributed backtracking on non-observable actions. A general theorem reduces the correctness of the latter partially reversible code to the causal correctness of the former [START_REF] Danos | Transactions in RCCS[END_REF]. In many transactional examples, this structured programming method works well, and obtains codes which are smaller, and simpler to understand [START_REF] Vincent Danos | Self-assembling trees[END_REF]. It also seems interesting from a correctness perspective, since one never has to deal with the full state space, and it is enough to consider the much smaller state space of the forward code causal compression relative to observable actions. Thus it obtains codes which are also easier to prove correct. It is only natural then to ask whether and to which extent such indirect correctness proofs can be automated. This is the question we address in this paper. Specifically we propose an algorithm, which, under certain rather mild assumptions about the system of interest, will compute its causal compression relative to a choice of observables. The true concurrency semantics tradition of using event structures as an intrinsic process representation comes to the rescue here. Besides event structures are uniquely suited to the handling of causal relationships between various events triggered by a process [START_REF] Winskel | Event structure semantics for CCS and related languages[END_REF]. For these reasons our procedure includes a translation of the process as a recursive flow event structure, and computes the relative causal compression on this intermediate representation. Benchmarks given for the classical example of the dining philosophers show a significant state compression, and a relatively low cost incurred by compression. Direct programming generates a state space that is already too big for being constructed by bisimulation verifiers for 6 agents, whereas our method can go well beyond 15. The language we use to formalize concurrent systems is the Calculus of Communicating Systems (CCS) [START_REF] Milner | Communication and Concurrency[END_REF]. This is a slightly more expressive language than basic models of communicating automata, in that processes can dynamically fork. On the other hand, this communication model includes no name-passing, which is a severe limitation in some applications.

Section 2 starts with a quick recall of CCS [START_REF] Milner | Communication and Concurrency[END_REF]. Section 3 develops its reversible variant RCCS, together with the central notion of causal correctness, and the fundamental result connecting causal correctness of a CCS process and full correctness of its lifting as a partially reversible process in RCCS [START_REF] Danos | Transactions in RCCS[END_REF]. The relative causal compression algorithm, and the accompanying verification method are explained in Section 4. Section 5 compares this method with the traditional direct method, using the dining philosphers problem as a benchmark. The conclusion discusses related work and further directions.

CCS

Syntax

CCS processes interact through binary communications on named channels: an output on channel x is written x, an input on the same channel is simply written x.

Processes p, q ::= a.p

| (p | q) | p + q | D(x) := p | (x)p | 0
We write P for the set of processes, A for the set of actions, and A * for the free monoid of action words. Restriction (x)p binds x in p and the set of free names of p is defined accordingly. In a recursive definition D(x) := p free names of p have to be x.

Operational semantics

A labelled transition system (LTS) is a tuple S, s, L, → where S is called the state space, s the initial state, L the set of labels, and → ⊆ S × L × S the transition relation. One uses the common notation s → a t, and for m = a 1 . . . a n ∈ A * , s → * m t means s → a1 s 1 , . . . , s n-1 → an t for some states s 1 , . . . , s n-1 . The operational semantics of a CCS term p is given by means of such an LTS (P, p, A, →), written TS(p), where → is given inductively by the rules:

a.p + q → a p (act) p → a p q → ā q p | q → τ p | q (synch) p → a p p | q → a p | q (par) p → a p a ∈ {x, x} (x)p → a (x)p (res) p ≡ p → a q ≡ q p → a q (equiv)
The equivalence relation ≡ is the classical structural congruence for choice and parallel composition, together with the recursion unfolding rule D(ỹ) ≡ p {ỹ/x} if D(x) := p.

Process equivalence

Given a set of observable actions, a basic requirement is to decide whether two processes have an equivalent interaction capacity with their environment. Several variants of observational equivalence for CCS processes have been considered. We use here a variant of weak bisimulation based on the choice of a countable distinguished subset K of the set of actions A, which we fix here once and for all. Actions in K are called observable actions. The complement A \ K of nonobservable actions is denoted by K c and also taken to be countable. Let S 1 = (S 1 , s 1 , A, →) and S 2 = (S 2 , s 2 , A, →) be LTSs both with labels in A, a relation R over S 1 × S 2 is said to be a weak simulation between S 1 , S 2 , if s 1 R s 2 and whenever p 1 R p 2 : -if p 1 → a q 1 , a ∈ K c , then p 2 → * m q 2 with m ∈ (K c ) * , and q 1 R q 2 ; -if p 1 → a q 1 , a ∈ K, then p 2 → * m q 2 with m ∈ (K c ) * a(K c ) * , and q 1 R q 2 . The idea is that S 2 has to simulate the behaviour of S 1 regarding observable actions, but is free to use any sequence of non observable ones in so doing. Such a relation R is said to be a weak bisimulation if both R and its inverse R -1 are weak simulations. When there is such a relation, S 1 and S 2 are said to be bisimilar, and one writes S 1 ∼ S 2 . A CCS process p is said to be a correct implementation of a specification LTS S, if TS(p) ∼ S. When the specification is clear from the context, we may simply say p is correct. One thing to keep in mind is that all these definitions are relative to a choice of K. Usually, K is taken to be A \ {τ }, but this more flexible definition will prove convenient.

We turn now to a quick intuitive introduction to RCCS. Consider the following CCS process:

(x) x | x | x.x.a.p | x.x.b.q (1)
Both subprocesses a.p and b.q require two communications on x to execute, so the whole process may reach a deadlocked state (x) x.a.p | x.b.q where neither a nor b may be triggered. If the intention is that the system implements the mutual exclusion process a.p + b.q, a possible fix is to give both subprocesses the possibility to release x:

(x) x | x | R p (x, a) | R q (x, a) (2) with R p (x, a) := x. τ.(R p (x, a) | x) + x.(τ.(R p (x, a) | x | x) + a.p) .
This example helps in realising two key things: first the original code (1) although not correct, is partially correct in the sense that any successful action a or b leads to a correct state p or q; second the proposed fix can be made an instance of a generic distributed backtracking mechanism. The idea of RCCS is to provide such a mechanism, in a way that partial or causal correctness (yet to be defined formally) in CCS, can be proved to be equivalent to full correctness of the same process once lifted to RCCS [START_REF] Danos | Reversible communicating systems[END_REF].

Syntax

RCCS forward actions are the same actions as CCS, namely A. Recall these are split into K and its complement K c . In the RCCS context actions in K are also called irreversible, or sometimes commit actions (following the transaction terminology); actions in K c are also called reversible, since these are the ones one wants to backtrack. RCCS therefore also has backward actions written a -, with a ∈ K c .

RCCS processes are composed of threads of the form m p, where m is a memory, and p is a plain CCS process:

r ::= m p | (r | r) | (x)r
Memories are stacks used to record past interactions:

m ::= θ, a, p • m | θ • m | ↑ • m |
where θ is a thread identifier drawn from a countable set. Open memory elements θ, a, p are used for reversible actions and contain a thread identifier θ, the action last taken, and the alternative process that was left over by a choice if any.

Closed memory elements θ are used for irreversible actions, and only contain an identifier. Eventually the memory element ↑ keeps track of the forking structure of processes (see congruence rules below)1 . The prefix relation on memories is defined as m m if there is an m such that m • m = m .

Processes are considered up to the usual congruence for parallel composition together with the following specific rules:

m D(ỹ) ≡ m p {ỹ/x} if D(x) := p m (p | q) ≡ ( ↑ • m p) | ( ↑ • m q) m (x)p ≡ (x)(m p) if x ∈ m
Any CCS process p can be lifted to RCCS with an empty memory (p) := p, and conversely, there is a natural forgetful map ϕ erasing memories and mapping back RCCS to CCS. Clearly ϕ( (p)) = p. When we want to insist that the lift operation is parameterised by the set K, we write K (p).

Operational semantics

The operational semantics of RCCS is also given as an LTS with transitions given inductively by the rules:

a ∈ K c θ ∈ m m a.p + q → θ a θ, a, q • m p (act) a ∈ K c θ, a, q • m p → θ - a m a.p + q (act * ) k ∈ K θ ∈ m m k.p + q → θ k θ • m p (commit) r → Θ a r Θ ∈ s r | s → Θ a r | s (par) r → Θ a r s → Θ ā s r | s → Θ τ r | s (synch) r → Θ a r a = x, x (x)r → Θ a r (res) r ≡ r → Θ a s ≡ s r → Θ a s (equiv)
In the contextual rules Θ stands either for θ or θ -. The freshness of the thread identifier θ is guaranteed by the side conditions θ ∈ m in the (act) and (commit) rules, and θ ∈ s in the (par) rule. The use of such identifiers corresponds to the notation introduced in Ref. [START_REF] Phillips | Reversing algebraic process calculi[END_REF] and equivalent to the one introduced originally for RCCS [START_REF] Danos | Transactions in RCCS[END_REF], as shown in Ref. [START_REF] Krivine | Algèbres de Processus Réversibles[END_REF]. Note that backtracking as defined in the operational semantics is a binary communication mechanism of exactly the same nature as usual forward communication. However, since threads are required to backtrack with the exact same thread with which they communicated earlier, backtrack can be shown to be confluent, at least for those processes that are reachable from the lifting of a CCS process. The (commit) rule uses a closed memory element θ • m indicating that the information contained in m is no longer needed, since by definition actions in K are not backtrackable. Supposing r is a process where any recursive process definition is guarded by a commit, an assumption to which we will return later on, this bounds the total size of open memory elements in any process reachable from r.

The fundamental property

The question is now to determine what are the possible (definitive) interactions of K (p) with the context. A first approach would be to find, for each p, a specification that would be bisimilar to the LTS engendered by K (p) (in which we would not observe θ on transitions). But then RCCS would be mere progress over CCS with explicit backtracking since one would need to consider every transitions of K (p) to check for bisimulation.

The question is now to see whether it is possible to obtain a characterisation of the behaviour of a lifted process K (p) solely in terms of p. Intuitively, K (p) being p enriched with a mechanism for escaping computations not leading to any observable actions, one might think that K (p) is bisimilar to the transition system generated by those traces of p which lead to an observable action. This is almost true.

To give a precise statement, we need first a few notations and definitions. An RCCS transition as defined above is fully described by a tuple t = r, a, Θ, r where r is the source of t, r its target, a its label and Θ its identifier. If a ∈ K we say that t is a commit transition, otherwise it is a reversible transition. If Θ = θ (Θ = θ -) we say t is forward (backward ). A trace is a sequence of composable transitions, and we write r → * σ s (p → * σ q) whenever σ is an RCCS (CCS) trace with source r (p) and target s (q). A trace is said to be forward if it contains only forward transitions.

A final and key ingredient is the notion of causality between transitions in a given forward trace. For CCS this is usually defined using the so-called proof terms [START_REF] Boudol | Permutation of transitions: An event structure semantics for CCS and SCCS. In Linear Time, Branching Time and Partial Order in Logics and Models for Concurrency[END_REF], but one can also use RCCS memories.

The set of memories involved in a forward transition t = r, a, θ, r is defined as µ(t) := {m ∈ r | ∃a, q : θ, a, q .m ∈ r }; this is either a singleton, if no communication happened, or a two elements set, if some did.

Definition 1 (Causality). Let σ : t 1 ; . . . ; t n be a forward RCCS trace: -t i and t j with i < j, are in direct causality relation, written t i < 1 t j if there is m ∈ µ(t i ), m ∈ µ(t j ) such that m m ; one says that t i causes t j , written t i < t j , if t i < * 1 t j . -σ is said to be causal if for all transitions t i with i < n, t i < t n ; it is said to be k-causal if it is causal, its last transition t n is labelled with k ∈ K, and all preceding transitions are labelled in K c .

One extends this terminology to CCS traces by saying a CCS trace p → * σ p is causal, if it lifts to a causal trace K (p) → * σ r with ϕ(r ) = p . We are now ready to state the property that explicits the effects of a commit transitions. Say a process r is initial if there is no possible backward transition with source r.

Proposition 1 (Committed trace.). Let σ be a RCCS trace with an initial source. Then the target of σ is also initial if and only if σ is in k-causal form.

This property says that whenever a commit is taken, then any other process can no longer cancel an action that played a role in the transaction; this expresses the durability of the transaction.

With the notion of causal trace in place, we can define the causal compression of a process p relative to K.

Definition 2 (Relative causal compression). Let p be a CCS process, its causal compression relative to K, written CTS K (p), is the LTS P, p, K, where k is defined as q k q if q → * σ q for some k-causal trace σ. We are now ready to state the theorem that characterizes the behaviour of K (p) in terms of the simpler process p.

Theorem 1 ([DK05]). Let TS K (p) := R, K (p), A, → be the LTS associated to the lift K (p), TS K (p) ∼ CTS K (p).
As said above, it is not true that TS K (p) is bisimilar to the transition system of traces of p leading to observable actions, one has to be careful to restrict to causal traces. A trivial but useful rephrasing of this result is:

Corollary 1. Let p be a CCS process, and S be its specification, if

CTS K (p) ∼ S then K (p) ∼ S.
In words, this says that to check the correctness of K (p) with respect to S, it is enough to check the correctness of CTS K (p).

If one goes back to the example at the beginning of this section, this says that {a,b} ((x) x | x | x.x.a.p | x.x.b.q ) is equivalent to a.p + b.q, as long as the causal compression of p = (x) x | x | x.x.a.p | x.x.b.q relative to {a, b} is. This is easily seen in this example, and in fact, as often in practice, CTS K (p) and S turn out to be equal.

The interest of this fundamental property lies in the fact that the causal compression relative to K, CTS K (p), is significantly smaller than the partially reversible process K (p). A natural question is therefore, given a process p, to compute CTS K (p). By finding an efficient way to do this, one would obtain an efficient verification procedure. This is the object of the next section.

Causal compression

A first idea to extract the causal transition system of a process p is to use the LTS generated by (p) and screen off non causal traces. One cannot know however whether a trace can be extended into a k-causal form until a commit is effectively taken, and such an approach would likely lead to both superfluous (because lots of traces will not be causal) and redundant (because of trace equivalence) computations. A more astute approach is to look only at traces that will eventually be in a k-causal form. This requires a bottom up view of traces where one starts from commits inside a term, and then reconstructs causal traces triggering this commit by consuming its predecessors in every possible way.

However, there is no need to work directly in the syntax, and event structures [START_REF] Winskel | Event structure semantics for CCS and related languages[END_REF] provide exactly what is needed here: a truly concurrent semantics that abstracts from the interleaving of concurrent transitions, and more importantly an explicit notion of causality. Among the various types of event structures the most often considered are prime ones, because consistent runs can be simply characterized. Yet they lead to quite large data structures. 2 Our algorithm uses instead flow event structures (FES) [START_REF] Boudol | Permutation of transitions: An event structure semantics for CCS and SCCS. In Linear Time, Branching Time and Partial Order in Logics and Models for Concurrency[END_REF][START_REF] Boudol | Flow event structures and flow nets[END_REF][START_REF] Van Glabeek | Well-behaved flow event structures for parallel composition and action refinement[END_REF]]. On the one hand, there is a simple inductive translation of CCS terms into FESs that incurs no computational cost; on the other hand, FES are algorithmically convenient compact forms of event structures.

We first explain how to extract the causal compression CTS K (p) from the translation of p into an FES. Then we discuss computational issues such as how to make this an algorithm, and how some of the apparent computational costs can be circumvented at the level of the implementation.

Flow event structures

A (labelled) flow event structure is a tuple E = E, ≺, #, λ where -E is a set of events, -≺ ⊆ E × E is the flow relation which has to be irreflexive, -# ⊆ E × E is the conflict relation which is symmetric, -and λ : E → A a labelling function. The idea is that the flow relation gives all immediate possible causes of an event, while the conflict relation indicates a conflicting choice between two events.

Definition 3. Let E = E, ≺, #, λ be an FES, a set X ⊆ E is a configuration of E, written X ∈ C(E), if it is: -conflict free: # ∩ (X × X) = ∅,
cycle free: ≺ * /X is a partial order, -and left-closed up to conflicts: if e ∈ X and there is d ∈ E such that d ≺ e then either d ∈ X or there exists f ∈ X such that f ≺ e and f #d.

The last two conditions are the price to pay for working with FESs, and are not needed for prime ones. The first one will require some optimised structuring of the conflict relation, we'll return to this point soon.

A configuration X in E with e ∈ X is e-minimal if ∀e ∈ X : e ≺ * e. The set of e-minimal configurations is denoted by C E, e .

There is an easy inductive translation u unfolding any CCS process into a FES, where events correspond to communications, and configurations are those subsets of events that a trace can trigger. We recall now this translation, defined by induction on the structure of the CCS process, from [START_REF] Boudol | Permutation of transitions: An event structure semantics for CCS and SCCS. In Linear Time, Branching Time and Partial Order in Logics and Models for Concurrency[END_REF]. def = E q , ≺ q , # q , λ q with E p ∩ E q = ∅, we have u(p + q) def = E p E q , ≺, #, λ where:

-(Prefix) Let u(p) def = E, ≺, #,
-e#e if (e, e ) ∈ E p × E q or if either e# p e or e# q e def = E q E q E, ≺, #, λ where:

-≺ / E p × E p def = ≺ p and ≺ / E q × E q def = ≺ q -(
-

E def = e (e ,e ) | (e , e ) ∈ E p × E q & λ p (e) = λ q (e ) -∀e ∈ E p E q E, λ(e) def = λ p (e) if e ∈ E p , λ(e) def = λ q (e) if e ∈ E q and λ(e) def = τ else.
-∀(e, e ) ∈ (E p E q E) 2 , we have e ≺ e if either:

-π 0 (e) ≺ p π 0 (e ) -π 1 (e) ≺ q π 1 (e ) -e#e if:

-π i (e) = π i (e ) for some i ∈ {1, 2} -π 0 (e)# p π 0 (e ) -π 1 (e)# q π 1 (e ) The correctness of the unfolding function u is given by the following representation theorem:

Theorem 2 ( [START_REF] Boudol | Flow event structures and flow nets[END_REF]). Let p be a CCS process, and T (p) stand for the traces of p quotiented by trace equivalence, then (T (p), ≤) and (C(u(p)), ⊆) are isomorphic.

One can define a transition system out of an FES. To do this, we define E|X, the residual of E by a configuration X in C(E).

Definition 4 (Residual)

. Let E = E, ≺, #, λ be an FES, X be a configuration of E, and define X # := {e ∈ E | ∃e ∈ X : e #e}. The residual of E by X is E|X := E , ≺ , # where:

E := E \ (X ∪ X # ) ≺ :=≺ ∩ (E × E ) # := # ∩ (E × E )
The LTS associated to E = E, ≺, #, λ has initial state E, and transition relation given by E

→ X E if X ∈ C(E ) and E = E |X.
It is here that our reframing of the compression question in terms of event structures pays off, since to obtain the causal compression of the transition system above, all one has to do is to restrict labels to e-minimal configurations such that λ(e) ∈ K. The causal LTS associated to E, written CTS K (E), has initial state E, and transition relation given by E k E if there is an event e ∈ E such that E → X E with X ∈ C E , e and λ(e) ∈ K. As a consequence of the representation theorem one gets: Lemma 1. Let p be a CCS process, then CTS K (p) and CTS K (u(p)) are isomorphic.

At this point, we have an equivalent definition of CTS K (p) in terms of the FES u(p), and it remains to see how one can turn this definition into an algorithm. This is what we discuss now.

Algorithmic discussion

First, the unfolding u(p) is in general an infinite object even if we restrict to finite state processes. To keep with finite internal data structures, we require each recursive process definition to be guarded by a commit action. This seems a reasonable constraint, in that there is a priori no reason to model a transactional mechanism with a process that allows infinite forward inconclusive traces.

To compute CTS K (u(p)), we use instead of u, a partial unfolding u fin that coincides with u except it does not unfold any recursive definition. The constraint above ensures that every commit k that is reachable by a single causal transition can be seen by this partial unfolding. Only after triggering the event corresponding to k, are the recursive calls guarded by k (if any) unfolded, and their translations by u fin added to the residual of the obtained event structure. One then checks whether the obtained residual event structure is isomorphic with some obtained previously, and adds it to the state space if not. Given a process p, the algorithm to compute CTS K (u(p)) proceeds as follows:

0. E = E, ≺, #, λ := u fin (p) 1.
For all e ∈ E such that λ(e) ∈ K, compute the e-minimal configurations X e ∈ C E, e . 2. For each such X e build the residual E|X e , with recursive definitions guarded by e unfolded using u fin .

Add the transitions E

k E|X e , where k = λ(e), to the CTS under construction. 4. For each residual E|X e not isomorphic to any previous one, set E := E|X e and goto step 1.

By the representation theorem, this algorithm will terminate as soon as CTS K (p) is finite. In practice most of the isomorphism tests can be avoided by using a quite discriminative equality test between FES signatures which is linear in the number of events. Another efficiency problem one has to deal with is the internal representation of the conflict relation (which is involved in step 1 because of the conflict-free condition on configurations). In prime event structures conflict is inherited by causality, that is to say if e#e and e ≺ e , then e#e . Hence a rather compact way to represent conflict is to keep only (e, e ) ∈ # and deduce when needed that e#e by heredity.

We have found that a similar compact structure, which we call a conflict tree can be used for FESs. Conflict trees are built during process partial unfoldings, and result in a typically logarithmically compact representation of conflict, for a low computational cost. An example of a conflict tree is given on the right: conflicts are predicated of intervals, and [nm]#[nm ] means that any pair of events indexed within {n, . . . , m} × {n , . . . , m } is in conflict.

Benchmark

To get a sense of how well our verification technique performs compared with straight bisimulation based verification, we ran several tests 2 using encodings the dining philosophers problem. This timeless example of distributed consensu involves n philosophers eating together around a table. Each of them needs tw chopsticks to start eating, and has to share them with his neighbours. Whe a philosopher has eaten, he releases his chopsticks after a while and goes bac to the initial state. In the partial implementation, say p part , once a philosophe takes a chopstick he never puts it back unless he has successfully eaten. In th fully correct one, say p full , he may release chopsticks at any time (thus avoidin deadlocks). The CCS processes p part and p full for n = 2 correspond roughly t

Benchmark

The relative compression algorithm was implemented as a prototype in Ocaml in order to get a sense of how well our verification technique performs compared with a straight bisimulation based verification. To do so we ran several tests 4 using encodings of the dining philosophers problem. This timeless example of distributed consensus involves n philosophers eating together around a table. Each of them needs two chopsticks to start eating, and has to share them with his neighbours. When a philosopher has eaten, he releases his chopsticks after a while and goes back to the initial state. In the partial implementation, say p part , once a philosopher takes a chopstick he never puts it back unless he has successfully eaten. In the fully correct one, say p full , he may release chopsticks at any time (thus avoiding deadlocks). The CCS processes p part and p full for n = 2 correspond roughly to the earlier examples (1) and ( 2). (See [START_REF] Danos | Transactions in RCCS[END_REF] for a general definition and detailed study.) There are two main reasons for taking the dining philosophers example. First it is a paradigmatic example of distributed consensus, so the way to solve it without access to the scheduler (by adding additional semaphores for instance) has to involve backtracking. Second, it turns out that the number of possible states of the specification is given by a Fibonacci sequence.

S(1) = 1 S(2) = 3 S(n + 1) = S(n) + S(n -1)
This is convenient in that it gives a simple means to compare the time of computation with the size of the specification state space. Verifying correctness of p full using the Mobility Workbench (MWB) [START_REF] Victor | The Mobility Workbench -a tool for the π-calculus[END_REF] (see top curve, Fig. 1) proved to be impossible beyond 5 philosophers (around 160 specification states) because of memory limitations. By using first the our prototype (see bottom curve, Fig. 1) to extract the causal transition system of p part , we could verify up to 19 philosophers (around 15, 000 specification states) within a time which stayed roughly proportional to the number of states. Since CTS(p part ) is in this case equal to the specification, the remaining part of the correctness proof takes negligible time (MWB needs 0.4s for 10 philosophers).

Conclusion

We have proposed a method for the verification of distributed systems which uses an algorithm of relative causal compression. The method does not always apply: the process one wants to verify must use a generic backtracking mechanism. This may seem a limitation, but it often obtains a much simpler code, and many examples of distributed transactions lend themselves naturally to this constraint. When the method does apply, however, it proves very effective as we have shown in the dining philosophers example.

State space explosion in automated bisimulation proofs is a well known phenomenon, and trace compression techniques have been proposed to avoid the redundancy created by the interleaving of transitions [START_REF] Boudol | Permutation of transitions: An event structure semantics for CCS and SCCS. In Linear Time, Branching Time and Partial Order in Logics and Models for Concurrency[END_REF][START_REF] Godefroid | Using partial orders for the efficient verification of deadlock freedom and safety properties[END_REF], and used in model-checking applications [BCDP95,AQR + 04]. These compressions preserve bisimilarity, whereas our does not, and is of a completely different nature. Besides, and because our algorithm uses event structures, we also benefit from this classical kind of compression.

Conclusion

We have proposed a method for the verification of distributed systems which uses an algorithm of relative causal compression. The method does not always apply: the process one wants to verify must use a generic backtracking mechanism. This may seem a limitation, but it often obtains a much simpler code, and many examples of distributed transactions lend themselves naturally to this constraint. When the method does apply, however, it proves very effective as we have shown in the dining philosophers example.

State space explosion in automated bisimulation proofs is a well known phenomenon, and trace compression techniques have been proposed to avoid the redundancy created by the interleaving of transitions [6,12], and used in modelchecking applications [13,14]. These compressions preserve bisimilarity, whereas our does not, and is of a completely different nature. Besides, and because our algorithm uses event structures, we also cash in on this classical kind of compression.

There is no reason why this verification method should be limited to CCS. Other concurrent models can be equipped with backtracking, and forward and backward aspects of correctness can be split there as well. Recent work extends the concept of partially reversible computations to various process algebras [15],

Relative causal compression using our algorithm:

Conclusion

We have proposed a method for the verification of distributed systems which uses an algorithm of relative causal compression. The method does not always apply: the process one wants to verify must use a generic backtracking mechanism. This may seem a limitation, but it often obtains a much simpler code, and many examples of distributed transactions lend themselves naturally to this constraint. When the method does apply, however, it proves very effective as we have shown in the dining philosophers example.

State space explosion in automated bisimulation proofs is a well known phenomenon, and trace compression techniques have been proposed to avoid the redundancy created by the interleaving of transitions [6,12], and used in modelchecking applications [13,14]. These compressions preserve bisimilarity, whereas our does not, and is of a completely different nature. Besides, and because our algorithm uses event structures, we also cash in on this classical kind of compression.

There is no reason why this verification method should be limited to CCS. Other concurrent models can be equipped with backtracking, and forward and backward aspects of correctness can be split there as well. Recent work extends the concept of partially reversible computations to various process algebras [15], There is no reason why this verification method should be limited to CCS. Other concurrent models can be equipped with backtracking, and forward and backward aspects of correctness can be split there as well. Recent work extends the concept of partially reversible computations to various process algebras [PU06,PU07,LMS10,LMSS11], and it is possible to define an analogue of RCCS for the π-calculus. New advances in event structure semantics for πcalculus [VY10,CVY12] might allow to extend the causal compression algorithm, so as to cover the important case of name-passing calculi.

2

  λ be the FES corresponding to P and let e ∈ E, then for any prefix action α we have u(α.p) def = E {e} , ≺ , #, λ where:-≺ / E × E def = ≺ and ∀e ∈ E, e ≺ e Specifically in prime event structure causes of an event must be uniquely determined, and this forces duplication of the future of an event each time it is engaged in a synchronization.-λ / E × E def = λ and λ(e) def = α -(Choice) Let u(p) def = E p , ≺ p , # p , λ p and u(q)

  Parallel product) with use the notation e (e0,e1) to denote the event resulting from the synchronization of events e 0 and e 1 . We use the traditional projection π 0 (e (e0,e1) ) def = e 0 and π 1 (e (e0,e1) ) def = e 1 . In addition for every event e that is not a synchronization we consider π 0 (e) def = π 1 (e) def = e. With these conventions, we define the unfolding u(p | q)

0 def=

 0 Consider for instance the process p = a.b.0 | b.ā.0. Using the above unfolding function one obtains the following FES, where arrows represent causality and dotted edges represent conflict: With the convention that λ(e α ) = α and λ(e α, ᾱ) = τ , for all α ∈ a, ā, b, b . It follows from Definition 3, that maximal configurations 3 of u(p) are: X {e a , e b , eb, e ā} X 1 def = e a , e b, b, e ā X 2 def = {eb, e a,ā , e b }

Fig. 5 .

 5 Fig. 5. Conflict tree of a 3 .(b 0 | c 2 + d 1 ) + e 4

Direct bisimulation test for p full :

Fig. 6 .Fig. 7 .

 67 Fig. 6. Direct bisimulation test for p full .

Fig. 6 .Fig. 7 .

 67 Fig. 6. Direct bisimulation test for p full .

Fig. 1 .

 1 Fig. 1. Benchmark results for the dining philosophers

The convention used here for keeping track of forking processes differs slightly from Ref.[START_REF] Danos | Reversible communicating systems[END_REF].

A configuration X is maximal if there is no additional event e such that e X is a valid configuration.

Tests were made with an Intel Pentium 4 CPU 3.20GHz with 1GB of RAM.

This work has been partially supported by the French National Research Agency (ANR), project REVER n • ANR 11 INSE 007