Loïc Peletan
email: loic.peletan@insa-lyon.fr

Sebastien Baguet
email: sebastien.baguet@insa-lyon.fr

Georges Jacquet-Richardet

Mohamed Torkhani
email: mohamed.torkhani@edf.fr

USE AND LIMITATIONS OF THE HARMONIC BALANCE METHOD FOR RUB-IMPACT PHENOMENA IN ROTOR-STATOR DYNAMICS

In the present paper, a Harmonic Balance Method (HBM) coupled with a pseudo-arc length continuation algorithm is presented for the prediction of the steady state behaviour of a rotorstator contact problem. The ability of the HBM to reproduce the four most common phenomena encountered during rotor to stator contact situations (i.e. 'no-rub', 'full annular rub', 'partial rub' and 'backward whirl/whip') is investigated. A modified Jeffcott rotor model is used and results of the proposed algorithm are compared with traditional time marching solutions and analytical predictions. The advantages and limitations of the HBM for this kind of problem are analyzed. It is shown that the HBM is orders of magnitude faster than transient simulations, and provides very accurate results. However, in its current form it is unable to predict quasi-periodic behaviour. Detailed analysis of the transient solutions yields valuable information for the future extension of the HBM to efficient quasi-periodic simulations.

INTRODUCTION

One of the objectives of turbomachinery designers is to improve the machines efficiency. This can be achieved by reducing the operating gap between the rotating and stationary parts of the rotor, but leads to an increased risk of rotor to stator con-tact. For safety reasons, designers must ensure that such interactions cannot have serious consequences. When an accurate prediction of the system response is required, it is insufficient to consider the linear response of the system alone. When nonlinearities are considered, highly complex dynamic behaviours are revealed, even for very simple rotor models. The simplest rotor model encountered in the literature is the Jeffcott rotor. Although this model has only one node, with two degrees of translational freedom (one for each radial direction), it may exhibit a rich variety of phenomena, due to the presence of nonlinearities. Edwards [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF] studied a Jeffcott rotor with an additional degree of torsional freedom, when subjected to a mass imbalance and contact with a circular, rigid, fixed stator. The contact was modelled using a penalty method with a Coulomb friction law. The system of equations was solved by numerical integration (Runge-Kutta). Edwards showed that even when the rotor is excited by a periodic force (i.e. an imbalance force) only, this relatively simple system can exhibit periodic as well as chaotic forms of behaviour.

Many authors have studied nonlinear rotor models, using time marching techniques to reach the steady state response of the system [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF][START_REF] Feng | Rubbing phenomena in rotor-stator contact[END_REF][START_REF] Popprath | Nonlinear dynamics of a rotor contacting an elastically suspended stator[END_REF][START_REF] Han | Periodic motions of a dual-disc rotor system with rub-impact at fixed limiter[END_REF][START_REF] An | Dynamic response of a rub-impact rotor system under axial thrust[END_REF]. This approach, although essential when investigating the system's transient behaviour [START_REF] Feng | Rubbing phenomena in rotor-stator contact[END_REF][START_REF] Grapis | Overcritical high-speed rotor systems, full annular rub and accident[END_REF][START_REF] Braut | Dynamic analysis of the rotorstator contact due to blade loss[END_REF], is not optimal in terms of CPU time when the steady state response is investigated, because many cycles need to be simulated before the transient behaviour disappears completely. This can be problematic when parametric studies are carried out and/or when more complex models (such as finite element rotor models [START_REF] Han | Periodic motions of a dual-disc rotor system with rub-impact at fixed limiter[END_REF][START_REF] Qin | Nonlinear responses of a rub-impact overhung rotor[END_REF]) are used.

On the other hand, some authors such as Jiang [START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear jeffcott rotor with cross-coupling stiffness coefficients[END_REF][START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] or Bently et al. [START_REF] Bently | Full annular rub in mechanical seals, part ii: Analytical study[END_REF] have used mathematical tools to study a Jeffcott rotor, for which they analytically determined the limits of existence for different types of behaviour, according to a set of dimensionless parameters such as the dimensionless rotational speed of the rotor or the coefficient of friction. Periodic and quasi-periodic phenomena were analyzed, and the periodic phenomena were found to be: 'no-rub motion', in which the rotor eccentricity is smaller than the operating gap, and 'full annular rub', in which the rotor is in permanent contact with the stator in a forward whirl mode, with a constant eccentricity. The quasiperiodic phenomena are 'partial rub', in which the rotor comes into intermittent contact with the stator, and self-excited 'backward whip' or 'backward whirl', in which the rotor rolls (respectively with or without sliding) on the inner surface of the stator, in a reverse whirl motion, with large eccentricities and large contact forces. The latter phenomenon is known to be quite violent and destructive. The four aforementioned phenomena are the most commonly reported in rotor-to-stator contact scenarios. These analytical methods are appealing, because they prove that several stable solutions can coexist with a given set of parameters. For instance, in [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] Jiang shows that a Jeffcott rotor can, for a given set of parameters and a certain range of rotational speeds, be in periodic full annular rub, as well as quasi-periodic backward whip. Time integration methods would lead to only one of these two solutions, depending on the initial conditions. These analytical methods have been extended to multi-mode rotor models by Childs and Bhattacharya [START_REF] Childs | Prediction of dry-friction whirl and whip between a rotor and a stator[END_REF]. Their predictions appear to be in good agreement with numerical simulations and experiments. Nevertheless, these analytical methods require non-trivial mathematical developments, and further extensions are still needed before they can be used to analyse more complex models.

Alternatively, although widely used to solve other types of problems, numerical techniques such as shooting or harmonic balance methods are rarely used for rotor-to-stator contact problems. Although such techniques do not require transient behaviour simulations, and are in practice much faster, they are limited to the study of periodic movements and, to a certain extent, to quasi-periodic movements. This limitation to periodic movements is the reason for which such techniques are generally limited to the study of geometric nonlinearities or oil-film types of mounting nonlinearities (such as in [START_REF] Sundararajan | An algorithm for response and stability of large order non-linear systemsapplication to rotor systems[END_REF]), because it is known that with this kind of nonlinearity the system response remains periodic. Nevertheless, for rotor-stator contact problems coupled with a path-following continuation scheme and stability analysis, these methods prove to be helpful in parametric studies, thus minimising the use of long transient simulations.

Following a brief presentation of the harmonic balance method, a Jeffcott rotor model described in [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] is studied, and the numerical results provided by this technique are compared to Jiang's predictions, as well as to the solutions given by traditional time marching techniques. The ability of the HBM to reproduce the four aforementioned phenomena is studied, for several values of the dry friction coefficient µ. The advantages and limitations of the proposed method are analyzed. The simulated transient behaviour of the rotor is then analysed spectrally, and valuable results are derived for a future extension of this HBM to quasiperiodic situations.

HARMONIC BALANCE METHOD (HBM)

The Harmonic Balance Method (HBM) is a well-known technique for the computation of periodic solutions for dynamic systems. This method consists in solving the equations of motion in the frequency domain, rather than in the time domain. The general time-domain equation for a dynamic mechanical system is given by (1). R(q(t)) = Kq(t) +C q(t) + M q(t) + f (q(t), q(t))p(t) = 0 [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF] where q represents the displacement vector for all degrees of freedom (DOFs); K, C and M are the generalized stiffness, damping and mass matrices; f is the nonlinear force vector and p is the external excitation force vector. When the external excitations are periodic, it can be assumed that a steady state solution exists for [START_REF] Edwards | The influence of torsion on rotor/stator contact in rotating machinery[END_REF], and that this solution is also periodic. The displacements can thus be written as a truncated Fourier series:

q(t) = Q 0 + N ∑ j=1 Q 2 j-1 cos(m j ωt) + Q 2 j sin(m j ωt) (2)
where Q i , i ∈ [1 . . . N] denote the Fourier coefficients of q; ω is the fundamental frequency of the external excitation, and m i , for i ∈ [1 . . . N], are positive integers. Similarly, f (q) and p(t) can be rewritten using the same form:

p(t) = P 0 + N ∑ j=1 P 2 j-1 cos(m j ωt) + P 2 j sin(m j ωt) (3)
f (t) = F 0 + N ∑ j=1 F 2 j-1 cos(m j ωt) + F 2 j sin(m j ωt) (4)
As described in [START_REF] Petrov | Analytical formula-tion of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks[END_REF], equations (2) (3) and ( 4) are substituted into (1), and a Galerkin procedure is applied to obtain the following algebraic system:

R(Q, ω) = Z(ω)Q + F(Q) -P = 0 (5)
where:

Z = diag K, K -m 2 1 ω 2 M m 1 ωC -m 1 ωC K -m 2 1 ω 2 M , . . . , K -m 2 N ω 2 M m N ωC -m N ωC K -m 2 N ω 2 M (6) Q = Q T 0 , Q T 1 , . . . , Q T 2N T , F = F T 0 , F T 1 , . . . , F T

2N

T and P = P T 0 , P T 1 , . . . , P T

2N

T are the Fourier coefficient vectors for displacements, nonlinear forces and external excitations, respectively.

Equation ( 5) needs to be solved for Q. This equation is still nonlinear, such that an appropriate nonlinear solver should be used to determine the correct solution. A Newton-Raphson solver was used in the present case, which solves ( 5) in an iterative manner:

Q (k+1) = Q (k) + ∂ R(Q (k) ) ∂ Q -1 R(Q (k) ) (7)
where

R(Q (k) , ω) = Z(ω)Q (k) + F(Q (k)
) -P is the residual of equation ( 5). It has been shown that the combination of a HBM and a Newton-Raphson solver is equivalent to the Incremental Harmonic Balance Method (IHBM) [START_REF] Lau | A variable parameter incrementation method for dynamic instability of linear and nonlinear elastic systems[END_REF]. As the relationship between the displacements and the nonlinear forces is known in the time domain only (a priori unknown in the frequency domain ( 5)), an Alternating Frequency Time (AFT) algorithm [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF] was used for the computation of F(Q). Direct and inverse Fourier transforms must be used in order to determine the solution. Displacements (Q) in the frequency domain are transferred to the time domain using an inverse Fourier transform. Once the displacements have been determined in the time domain (q), the corresponding nonlinear forces are calculated, following which a direct Fourier transformation is used to obtain the nonlinear forces in the frequency domain (F). The next Newton-Raphson iteration can then be carried out. In practice, fast Fourier transform (FFT) algorithms are used due to their combined speed, accuracy and robustness. This procedure can be summarized as follows:

Q(ω) f f t -1 -→ q(t) -→ f (t) f f t -→ F(ω)
Similarly, for reasons of improved accuracy and faster convergence, the Jacobian term

∂ R(Q (k) )
∂ Q is calculated using analytical derivation rather than finite differences. The path following technique used for the parametric study was a pseudo-arc length continuation method. Contrary to transient time marching techniques, the combination of HBM and a continuation method converges to both stable and unstable solutions, and there is no way to make an a priori assessment of the solution's stability. Several definitions of the stability of a dynamic solution exist. In the present study, we investigate the local stability of the solutions. The local stability assessment consists in applying a small perturbation to the equilibrium solution, and then verifying that this perturbation subsides over time. The stability can be computed a posteriori, using different theories. One theory frequently cited in the literature is that of Floquet. In the following, the stability of the HBM solutions is calculated in the time domain, by means of a monodromy matrix computation. The monodromy matrix is a linear operator describing the evolution over one period of perturbations, applied to the periodic solution. Its eigenvalues are called the Floquet multipliers. If any of the Floquet multipliers has a module greater than one, the solution is unstable, otherwise it is stable.

STUDY OF A NONLINEAR JEFFCOTT ROTOR Description of the problem

In this section, a nonlinear Jeffcott rotor identical to that described by Jiang in [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] is studied (see Fig. 1). The mathematical model has two degrees of freedom and the dimensionless equations of motion are as follows:

X + 2ξ X + β X + Θ 1 - R 0 R (X -µsign(V rel )Y ) = Ω 2 cos Ωτ (8) Y + 2ξY + βY + Θ 1 - R 0 R (Y + µsign(V rel )X)
= Ω 2 sin Ωτ (9) with the dimensionless variables:

V rel = R disk Ω + Rω b (10)
X = x e , Y = y e , R 0 = r 0 e , R disk = r disk e , V rel = v rel e , 2ξ = c √ k b m , β = k s k b , ω 2 = k b m , Ω = ω ω 2 , ω b = ω w ω 2 , R = √ X 2 +Y 2 , Θ = 1 if R ≥ R 0 , 0 else.
, τ = ω 2 t where x and y represent the radial displacements along the x and y axes, respectively, e is the mass eccentricity, r 0 is the initial clearance between the rotor and the stator, r disk is the disk radius, v rel is the relative velocity between the rotor and the stator, c is the damping factor, m is the rotor mass, k s and k b are the rotor and stator stiffnesses, respectively, ω is the rotor's rotational speed, and ω w is the rotor's whirl frequency.

As seen in the equations of motion ( 8) and ( 9) above, the stator model is circular, rigid and fixed. There are no DOF associated with it. The contact is modelled by a penalty law using a dry friction Coulomb model.

Jiang analytically the existence of four possible types of steady-state rotor behaviour, resulting from imbalance and rotor-to-stator contact. These four types of behaviour are: 'no-rub motion' in which there is no contact between the rotor and the stator. The orbit is circular and the movement is periodic. 'full annular rub' in which the rotor remains in permanent contact with the stator. The orbit is also circular and the movement periodic. 'partial rub' in which the rotor intermittently comes into contact with the stator. The movement is quasi-periodic. 'backward whip' in which the rotor rolls (with slippage) on the inner surface of the stator. This movement is a selfexcited phenomenon and its frequency does not necessarily depend on the rotor's rotational speed. As a consequence, the movement is a priori quasi-periodic.

The stability of these four types of behaviour can be derived from the aforementioned dimensionless parameters. More specifically, Fig. 2 shows the rotor response characteristics in the Ω-µ parameter plane. According to Jiang, the straight lines Ω l and Ω u indicate the rotational speed limits when the linear rotor eccentricity is equal to the initial clearance. Thus, the 'norub motion' can appear only at rotational speeds below Ω l and above Ω u . The curves HP and SN correspond to the limits of 'full annular rub' motion. To the right of the line HP, 'full annular rub' becomes unstable and the movement becomes quasiperiodic (Hopf bifurcation). To the right of line SN, the 'full annular rub' no longer exists (saddle node bifurcation). The line DF denotes the presence of a backward whip response. The line DW indicates the limit beyond which backward whip is triggered by the imbalance force.

There are some regions in which several behaviours coexist. The regions labelled '0' are regions where no-rub motion and dry whip coexist. In the region labelled '1', both full annular rub and dry friction backward whirl are possible. The region labelled '2' shows the presence of partial rub and backward whirl. In '3', full annular rub and no-rub motions are simultaneously present. Finally, a small region '4' exhibits no-rub motion, full annular rub and backward whip.

The results shown in Fig. 2 have been confirmed via numerical simulations. Applying the same parameter values as those used in Fig. 2, the simulations were run with a fixed value for the friction coefficient µ and varying values of Ω. Only one harmonic was retained for the HBM simulations. As the no-rub and full annular rub motions have circular orbits, one harmonic is sufficient to provide an exact representation of these movements. The pseudo arc length continuation algorithm was used for the parametric study with respect to Ω, over the range from 0 to 1.2. In the case of the simulations based on the time marching technique, the solutions were calculated for discrete values of Ω. The initial conditions used for the next value of Ω correspond to the final steady state conditions of the current simulation. This corresponds to an incremental continuation technique. For each value of Ω, 100 rotor revolutions were simulated in order to attain the steady-state behaviour. Twenty-five additional revolutions were then simulated and saved for analysis. This procedure was performed twice for each value of µ. In the first pass, the simulations were made with progressively greater values of Ω, whereas in the second pass progressively smaller values of Ω were used. As numerical simulations (HBM or time integration) require the use of non dimensionless variables, the following parameters were set arbitrarily: e = 0.1 m, k s = 100 N.m -1 , m = 1 kg. The other parameters were derived. All of the simulations were performed using the Code Aster [START_REF]Code aster: A general code for structural dynamics simulation under gnu gpl licence[END_REF] and Scipy [START_REF] Jones | Scipy: Open source scientific tools for Python[END_REF] environments. Although the simulations were run for many values of µ, only the results obtained for µ = 0.1 and µ = 0.2, which are representative of all possible solutions, are presented in the following.

Numerical Results

Results for µ µ µ = = = 0 0 0. . .1 1 1. Fig. 3 shows the response curve for µ = 0.1, β = 0.04, ξ = 0.05, R 0 = 1.05 and R disk = 20R 0 calculated using HBM and time integration. The y-axis represents the rotor eccentricity normalized to the initial clearance, on a logarithmic scale, and the x-axis corresponds to the normalized rotational speed Ω = ω ω 2 . In this case, the time transient and HBM curves are very similar. As predicted by Jiang's analytical results (see Fig. 2), only no-rub motion and full annular rub are found for µ = 0.1. The transition between no-rub motion and full annular rub (which correspond to the Ω l line in Fig. 2) occurs between approximately Ω = 0.150 and Ω = 0.155, which is highly consistent with the analytically determined value (Ω l = 0.1536). One can indeed see that the normalized rotor eccentricity is less than unity for Ω ≤ 0.150 and greater than unity for Ω ≥ 0.155. The Ω u line is also correctly predicted by both HBM and the time integration methods (Ω u = 0.854). In addition, the region labelled '3' in Fig. 2 is also correctly predicted. Fig. 3 shows that no-rub motion and full annular rub are both possible for 0.850 < Ω < 0.995. The time integration method is characterised by a sudden 'drop' in eccentricity at Ω SN ≈ 0.995, corresponding to the 'run up' simulation, whereas a 'jump' occurs at Ω u ≈ 0.850 during the 'run down' simulation. The HBM coupled to arc length continuation and stability computation is characterised by three possible solutions within this same range of rotational speeds, although only two of these are stable (passage through a saddle node bifurcation). The Floquet multipliers of the unstable HBM solution at Ω ≈ 0.98 are plotted in Fig. 4. As the number of Floquet multipliers is equal to twice the number of degrees of freedom, four Floquet multipliers are visible. One of these has departed from the unit circle on the positive real axis, which corresponds here to a saddle node bifurcation. Although the HBM and time integration methods described here provide extremely similar results, the HBM technique is approximately 500 times faster than the time integration method. The circular orbit of the rotor during full annular rub (Ω = 0.5) is shown in Fig. 5. Of course, the no-rub motion is also characterised by a perfectly circular orbit.

Results for µ µ µ = = = 0 0 0. . .2 2 2 The response curve for µ = 0.2 is shown in Fig. 6. Here, the HBM curve is very similar to that obtained for µ = 0.1. However, stability analysis reveals a loss of stability beyond Ω HP ≈ 0.289. The Floquet multipliers of the unstable solution at Ω ≈ 0.3 can be seen Fig. 7. A pair of Floquet multipliers has departed from the left-hand side of the unit circle, with a non-zero imaginary component indicating a Hopf bifurcation and a loss of stability. A Hopf bifurcation corresponds to the appearance of a second frequency in the real rotor response, and the initiation of a quasi-periodic regime. This means that although, theoretically, the periodic full annular rub response still exists beyond this point, it is in practice no longer possible to ilar to those observed for µ = 0.1. The no-rub motion is the first to be observed, and is very close to the HBM prediction. Then, beyond Ω l ≈ 0.15, the full annular rub regime starts. Here again, the results are very close to those found with the HBM simulations. However, beyond Ω HP ≈ 0.289, the time transient line deviates from that predicted by the HBM, and the rotor motion becomes quasi-periodic. The rotor orbit computed using the time transient technique at Ω = 0.33 is shown Fig. 9. It can be seen that the rotor does not remain in permanent contact with the stator, which corresponds to the 'partial rub' motion predicted by Jiang. During this movement, the maximum rotor eccentricity is greater than that predicted during full annular rub. At Ω DW ≈ 0.460, it can be observed that the eccentricity suddenly increases by a factor of ten. Although slightly different from the value found by Jiang (≈ 0.41), the latter value of Ω is coherent with this author's prediction. In this regime, referred to as backward dry whip, the rotor's orbit is almost circular, as shown in Fig. 10. The rotor can be seen to remain in the backward whip regime until the end of the sweep up.

A substantially different behaviour is observed during the sweep down. The rotor starts in the no-rub regime, until Ω u ≈ 0.86, at which a sudden increase in eccentricity occurs, transferring the rotor motion to the backward whip regime. The rotor remains in this regime until a very small value of Ω DF ≈ 0.045 is reached, at which the eccentricity suddenly drops, leading to the no rub regime. This value of Ω DF is also in good agreement with Jiang's prediction (see Fig. 8).

From the preceding description, it can be seen that the HBM is able to predict only a portion of the time transient simulations. However, stability analysis can detect unstable solutions, such that the DF curve can be determined in addition to the Ω l and Ω u lines. The DW and DF lines cannot be predicted by the HBM in its current form. Nonetheless, the HBM technique could be extended to quasi-periodic analysis [START_REF] Lau | Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems[END_REF][START_REF] Roques | n-dimensional harmonic balance method extended to non-explicit nonlinearities[END_REF], theoretically allowing the partial rub and backward whip regimes to be predicted.

Transient simulations can provide valuable information concerning the nature of the quasi-periodic regimes. A fast Fourier transform of the nonlinear contact force, during partial rub at Ω = 0.33, is shown Fig. 11. By analyzing this diagram, it can be concluded that the peaks are present for linear combinations of two incommensurable frequencies. The lowest frequency (noted f 1 ) is naturally the rotational frequency of the rotor (1 × Ω). The second frequency (noted f 2 ) cannot be accurately determined, since it is not commensurable with f 1 . However, it can be approximated, for this value of Ω, to approximately 1.08 × Ω. Spectral analysis of the contact force during partial rub for several values of Ω show that, even though f 1 is naturally always equal to the rotating frequency ( f 1 = 1 × Ω), the second frequency f 2 is related to f 1 ( f 2 = f ( f 1 )). Although we noted that this second frequency f 2 decreases when Ω increases, no trivial relationship could be found relating them to each other. It was verified that, for every simulated partial rub movement, spectral peaks appeared in the contact force power spectrum, at linear combinations of the two incommensurable frequencies, as described by the following formulae:

f n 0 = (2n + 1) f 1 + (n + 0) f 2 for n ∈ N (11)
f n 1 = (2n + 1) f 1 + (n + 1) f 2 for n ∈ N (12)
The peaks indeed occur in pairs, and decrease in amplitude with n (see Fig. 11). The values of f n 0 and f n 1 are given in Table 1 for increasing values of n, (to be compared with Fig. 11). The spectral analysis provided in Fig. 12 shows that the rotor movement has one main frequency during backward whip. This predominant frequency does not correspond to the rotational speed of the rotor, but rather to the backward whip frequency. The rotational frequency of the rotor f 1 has a comparatively small amplitude, because the backward whip phenomenon is self-excited. Our computations show that the backward whip frequency has little to no dependence on Ω. For Ω ≤ 0.5 the backward whip frequency increases slightly when Ω increases (from 0.75 × ω 2 at Ω = 0.045 to 1 × ω 2 at Ω = 0.5). For Ω ≥ 0.5 the backward whip frequency is constant and equal to ω 2 . This type of backward whip frequency behaviour is frequently reported in the literature, as for example in [START_REF] Childs | Prediction of dry-friction whirl and whip between a rotor and a stator[END_REF].

Other simulations were made for a large set of values of µ. For each value tested, good agreement was found between the HBM, the time transient computations, and Jiang's prediction. Unfortunately, Jiang's predictions are limited to Jeffcott rotor models. Nonetheless, the HBM and time transient tools can be used for more complex models, such as finite element rotor models, with no need for any further developments.

CONCLUSION

The applicability of the Harmonic Balance Method (HBM), with pseudo arc-length continuation and stability assessment, to the computation of the steady-state behaviour of a rotor to stator contact problem, has been investigated. The HBM results have been compared with classical time marching solutions and analytical results [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF]. The rotor model was that of a simple Jeffcott rotor, with a circular, rigid, static stator. It has been demonstrated that although the rotor is submitted to imbalance only, there are four possible kinds of solution: two of these are periodic solutions (no-rub motion and full annular rub), whereas the remaining two are quasi-periodic (partial rub and backward whip). The HBM is shown to be capable (at speeds orders of magnitude faster than those given by time marching analysis) of accurately predicting the two periodic solutions. The HBM can also predict the initiation of the quasi-periodic partial rub regime, which corresponds to a loss of stability of the periodic solution. However, the partial rub and backward whip motions cannot be described by means of the HBM in the presented form. An extension of this method to the treatment of quasi-periodic solutions is mandatory [START_REF] Lau | Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems[END_REF][START_REF] Roques | n-dimensional harmonic balance method extended to non-explicit nonlinearities[END_REF], and is currently under study. Our detailed analysis of the partial rub and backward whip regimes, using time marching algorithms, has provided clues for the implementation of this extension.

FIGURE 1 .

 1 FIGURE 1. JIANG'S SCHEMATIC DIAGRAMS OF THE JEFF-COTT ROTOR WITH A CIRCULAR RIGID STATOR [10].

FIGURE 2 .

 2 FIGURE 2. ROTOR RESPONSE CHARACTERISTICS IN THE PARAMETER PLANE OF Ω-µ WITH β = 0.04, ξ = 0.05, R 0 = 1.05 AND R disk = 20R 0 , ACCORDING TO JIANG [10].

FIGURE 3 .

 3 FIGURE 3. RESPONSE CURVE FOR µ = 0.1. HBM STABLE( ), HBM UNSTABLE ( ), TIME INTEGRATION (×).

FIGURE 4 .

 4 FIGURE 4. REAL AND IMAGINARY PARTS OF THE FLOQUET MULTIPLIERS IN THE UNSTABLE HBM SOLUTION AT Ω ≈ 0.98 (SADDLE NODE BIFURCATION)

FIGURE 5 .

 5 FIGURE 5. PERFECTLY CIRCULAR ORBIT OF THE ROTOR DURING FULL ANNULAR RUB AT Ω = 0.5.

FIGURE 6 .

 6 FIGURE 6. RESPONSE CURVE FOR µ = 0.2. HBM STABLE ( ), HBM UNSTABLE ( ), TIME INTEGRATION (×).

FIGURE 7 .

 7 FIGURE 7. REAL AND IMAGINARY PARTS OF THE FLOQUET MULTIPLIERS IN THE UNSTABLE HBM SOLUTION AT Ω ≈ 0.33 (SECONDARY HOPF BIFURCATION).

FIGURE 8 .

 8 FIGURE 8. SAME AS FIG. 2 WITH NOTEWORTHY VALUES OF Ω HIGHLIGHTED FOR µ = 0.1 AND µ = 0.2 [10].

FIGURE 9 .

 9 FIGURE 9. ROTOR ORBIT DURING PARTIAL RUB AT Ω = 0.33 (TIME TRANSIENT RESULT).

FIGURE 10 .

 10 FIGURE 10. BACKWARD WHIP ORBIT AT Ω = 0.33. THE ROTOR ECCENTRICITY IS ONE ORDER OF MAGNITUDE GREATER THAN THE CLEARANCE (TIME TRANSIENT RE-SULT).

FIGURE 11 .

 11 FIGURE 11. FAST FOURIER TRANSFORM OF THE NONLIN-EAR CONTACT FORCE SIGNAL DURING PARTIAL RUB AT Ω = 0.33. PEAKS AT LINEAR COMBINATIONS OF THE TWO INCOM-MENSURABLE FREQUENCIES ARE VISIBLE.

FIGURE 12 .

 12 FIGURE 12. FAST FOURIER TRANSFORM OF THE NONLIN-EAR CONTACT FORCE SIGNAL DURING BACKWARD WHIP Ω = 0.33. ONLY TWO PEAKS CAN BE SEEN.

TABLE 1 .

 1 f n 0 AND f n 1 CALCULATED FOR THE FIRST FOUR VALUES OF n FOR THE PARTIAL RUB MOTION AT Ω = 0.33. HERE, f 1 = 1 × Ω AND f 2 ≈ 1.08 × Ω.

	n f n 0 (×Ω) f n 1 (×Ω)
	0	1	2.08
	1	4.08	5.16
	2	7.16	8.24
	3	10.2	11.3

ACKNOWLEDGMENT

This work was partially supported by the French National Agency (ANR) in the framework of its Technological Research COSINUS program (IRINA, project ANR 09 COSI 008 01 IRINA)


Mzaki Dakel
email: mzaki.dakel@insa-lyon.fr

Sébastien Baguet
email: sebastien.baguet@insa-lyon.fr

Régis Dufour
email: regis.dufour@insa-lyon.fr

INVESTIGATION ON THE DYNAMICS OF AN ON-BOARD ROTOR-BEARING SYSTEM

In ship and aircraft turbine rotors, the rotating mass unbalance and the different movements of the rotor base are among the main causes of vibrations in bending. The goal of this paper is to investigate the dynamic behavior of an onboard rotor under rigid base excitations. The modeling takes into consideration six types of base deterministic motions (rotations and translations) when the kinetic and strain energies in addition to the virtual work of the rotating flexible rotor components are computed. The finite element method is used in the rotor modeling by employing the Timoshenko beam theory. The proposed on-board rotor model takes into account the rotary inertia, the gyroscopic inertia, the shear deformation of shaft as well as the geometric asymmetry of shaft and/or rigid disk. The Lagrange's equations are applied to establish the differential equations of the rotor in bending with respect to the rigid base which represents a noninertial reference frame. The linear equations of motion display periodic parametric coefficients due to the asymmetry of the rotor and time-varying parametric coefficients due to the base rotational motions. In the proposed applications, the rotor mounted on rigid/elastic bearings is excited by a rotating mass unbalance associated with sinusoidal vibrations of the rigid base. The dynamic behavior of the rotor is analyzed by means of orbits of the rotor as well as fast Fourier transforms (FFTs).

INTRODUCTION

The rotating machines are among the vital parts of modern industrial domains and the rotor forms their key component. Base excitations of rotating machines can be caused by a machine carried on a base affected by ground movement or onboard moving systems. Rotors mounted on transportation systems are mainly excited by the rotating mass unbalance and the various motions of their base which can increase their lateral vibration and create a dynamic instability phenomenon.

The rotors are complex systems and their analysis is also complex. Many studies concentrated on making a clear understanding of the dynamics of symmetric/asymmetric rigid/flexible rotor systems supported by linear/nonlinear elastic bearings in the case of fixed base because these systems can be introduced in the modern industry [START_REF] Rao | Rotordynamics[END_REF][START_REF] Lalanne | Rotordynamics Prediction in Engineering[END_REF]. Berlioz et al. [START_REF] Berlioz | Dynamic Behavior of a Drill-String: Experimental Investigation of Lateral Instabilities[END_REF] observed theoretically and experimentally the lateral instabilities of a drill-string parametrically excited under rotating conditions. Nandi and Neogy [START_REF] Nandi | An Efficient Scheme for Stability Analysis of Finite Element Asymmetric Rotor Models in a Rotating Frame[END_REF] presented an efficient analysis of stability for finite element models of asymmetric rotors and investigated whether an unstable rotor could be stabilized using an isotropic viscous damper.

Samali et al. [START_REF] Samali | Random Vibration of Rotating Machines Under Earthquake Excitations[END_REF] used the Monte Carlo simulation to simulate the nonstationary earthquake ground motions and to determine the statistics of rotating machinery response. In [START_REF] Hori | Earthquake-Induced Instability of a Rotor Supported by Oil Film Bearings[END_REF], a seismic response of a Jeffcott rotor supported by oil film bearings to a real seismic wave was examined and its stability was investigated by calculating the loci of the disk and journal centers. Bachelet et al. [START_REF] Bachelet | Rotor Under Seismic Excitation: A Spectral Approach[END_REF] investigated the dynamic behavior of an asymmetric rotor excited by a base seismic translation and developed an original spectral approach in order to approximate the rotor response. Da Silva Tuckmantel et al. [START_REF] Silva Tuckmantel | An Analysis on the Supporting Structure Representative Model in Rotating Systems[END_REF] represented the supporting structure (foundation) of a rotating system by coupled as well as uncoupled modes. They found that there were limitations in identifying the vibratory modes and their corresponding damping factors in the system response. Kang et al. [START_REF] Kang | An Investigation in Stiffness Effects on Dynamics of Rotor-Bearing-Foundation Systems[END_REF] and Cavalca et al. [START_REF] Cavalca | An Investigation on the Influence of the Supporting Structure on the Dynamics of the Rotor System[END_REF] studied the effect of the flexible foundation on the dynamic characteristics of the rotor-bearing systems.

The number of works dealing with the investigation of dynamic behavior of a rotor during a harmonic motion of its base is quite low. Duchemin et al. [START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF] applied the method of multiple scales to a simple rotor model under a base sinusoidal rotation in order to observe its stability. They presented experimental results in order to validate the analytical study. Driot et al. [START_REF] Driot | Theoretical and Experimental Analysis of a Base-Excited Rotor[END_REF] described the orbits of a rotor induced by a base harmonic rotational movement using numerical methods. From their work, the comparison between numerical results and experimental investigation was quite satisfactory. El-Saeidy and Sticher [START_REF] El-Saeidy | Dynamics of a Rigid Rotor Linear/Nonlinear Bearings System Subject to Rotating Unbalance and Base Excitations[END_REF] obtained the responses of a rigid rotor-bearing system subjected to rotating mass unbalance plus base harmonic excitations along or around lateral directions using analytical solutions in the case of linear bearing and a time integration scheme in the case of bearing cubic nonlinearity. Das et al. [START_REF] Das | Active Vibration Control of Unbalanced Flexible Rotor-Shaft Systems Parametrically Excited Due to Base Motion[END_REF] performed a numerical simulation of a flexible rotor system excited by mass unbalance and rotational periodic motion of the base in order to investigate its active vibration control. This control was successful for avoiding the lateral parametric instability due to the sinusoidal rotation of the base. Among all the literature mentioned above, there are references studying base-excited rotor systems and whose few works deal with the harmonically excited on-board rotors [START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF][START_REF] Driot | Theoretical and Experimental Analysis of a Base-Excited Rotor[END_REF][START_REF] El-Saeidy | Dynamics of a Rigid Rotor Linear/Nonlinear Bearings System Subject to Rotating Unbalance and Base Excitations[END_REF][START_REF] Das | Active Vibration Control of Unbalanced Flexible Rotor-Shaft Systems Parametrically Excited Due to Base Motion[END_REF]. Moreover, these references concentrate on the investigation of dynamic behavior of either simple rotors (for example, Jeffcott rotor [START_REF] Hori | Earthquake-Induced Instability of a Rotor Supported by Oil Film Bearings[END_REF] and rotors modeled using the Rayleigh-Ritz method [START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF][START_REF] Driot | Theoretical and Experimental Analysis of a Base-Excited Rotor[END_REF]), rotors supported by elastic bearings with constant damping and stiffness coefficients (for example, [START_REF] El-Saeidy | Dynamics of a Rigid Rotor Linear/Nonlinear Bearings System Subject to Rotating Unbalance and Base Excitations[END_REF][START_REF] Das | Active Vibration Control of Unbalanced Flexible Rotor-Shaft Systems Parametrically Excited Due to Base Motion[END_REF]) or rotors excited by simple motions (rotation or translation) of the base. As a consequence, the applications proposed in these works are not suitable for realistic ones. In this paper, an improved model is presented. Namely, an asymmetric rotor is discretized using the finite element method based on the Timoshenko beam theory, mounted on hydrodynamic bearings linearized with damping and stiffness coefficients calculated using the Reynolds equation [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF] and excited by base combined motions (six types of deterministic rotational and translational motions). The rotary inertia, the gyroscopic inertia, the shear deformation of shaft as well as the geometric asymmetry of shaft and/or rigid disk are taken into consideration. By computing the Rayleigh damping coefficients, the effect of rotor internal damping is included in the study. The application of the Lagrange's equations using the finite element method gives the differential equations of the rotor in bending with respect to a noninertial reference frame connected to the rotor rigid base. The linear equations of motion point out periodic parametric coefficients due to the geometric asymmetry of the rotor and time-varying parametric coefficients due to the base rotational excitations. These parametric coefficients are considered as producers of internal excitation and can create lateral dynamic instability. In the presented applications, the rotor mounted on rigid or elastic linear bearings is subjected to rotating mass unbalance plus to base sinusoidal rotation without or with base sinusoidal translation. Numerical solutions for the linear equations of vibratory motion of the on-board rotor are computed and analyzed by means of orbits of the rotor as well as fast Fourier transforms (FFTs).

THEORTICAL ANALYSIS OF AN ON-BOARD ROTOR-BEARING SYSTEM

The rotor generally consists of the disk, shaft, mass unbalance, bearing and base. This section includes formulations for the rotor components and the equations of motion.

Preliminary Calculations

Three principal coordinate systems shown in Fig. 1 are introduced to take into account the movement of the rotor rigid base. They are attached to the ground R g , the rigid base R and the moving rotor R l .

Figure 1. REFERENCE FRAMES OF THE ON-BOARD ROTOR

The rotational motions of the rotor base are described by the angular velocity vector components ω x (t), ω y (t) and ω z (t) of the rigid base R with respect to the ground R g projected in the frame R. The translational motions of the rotor base are described by the coordinates x O (t), y O (t) and z O (t) of the position vector O g O projected in the frame attached to the base R. The Euler angles ψ(y,t), θ(y,t) and φ (t) (Fig. 2) permit describing the orientation of the rotor R l with respect to its base R. The angular velocity vector of the rotor R l with respect to the ground R g measured in the frame R l is given by , ,

g g l l l l l l T R R R x y z R R R R = + = ω ω ω ω ω ω ( 1
)
where T is the matrix transposition symbol. The components

( l x ω , l y ω , l z
ω ) are formulated as a function of (ψ,θ, φ ) and their time derivative as well as (ω x ,ω y ,ω z ). The rotor is supposed to rotate at a constant speed Ω. So the spinning angle φ is replaced by Ωt and its derivative φ ɺ by Ω (• refers to the differentiation with respect to time t). Let us consider a generic point C 0 on the elastic line of the nondeformed shaft. Its coordinates in the frame R are (0,y,0). The translational displacements u(y,t) and w(y,t) of the point C 0 due to bending of the flexible rotor are expressed respectively with respect to the Ox and Oz axes of the coordinate system R.

Energy and Virtual Work Calculations

The kinetic and strain energies are measured from the ground and their terms are written with respect to the frame linked with the rotor rigid base R.

Disk. The disk is assumed to be rigid. Therefore only its kinetic energy T d is calculated as follows [START_REF] Bachelet | Rotor Under Seismic Excitation: A Spectral Approach[END_REF][START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF] ( ) ( ) ( )

1 2 2 g g g g l l l l d T T R R R R d d m O O R R m T = + v v ω I ω (2)
( ) ( ) ( ) 2 2 2 2 2 2 2 2 2 1 2 l l l l l l l l d d d d d O O O mo x z y y di x z m m m m T u v w I I I = + + + + + + - ɺ ɺ ɺ ω ω ω ω ω (4)
Shaft. The shaft is supposed to be flexible. Therefore it is characterized by the kinetic and strain energies and modeled by beam elements. The kinetic energy of a shaft can be obtained by considering a shaft elementary volume which can be considered as a disk of thickness dy. Thus the kinetic energy of a shaft takes the following form [START_REF] Bachelet | Rotor Under Seismic Excitation: A Spectral Approach[END_REF][START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF] ( )

( ) ( ) 2
= + +  + + +     + -    ∫ ∫ ∫ ∫ ɺ ɺ ɺ ρ ρ ω ω ρ ω ρ ω ω (5)
where ρ sh , S sh and l sh are respectively the density, the crosssection and the length of the shaft.

      ∂ ∂   = +       ∂ ∂             ∂ ∂   + + + -       ∂ ∂              ∂ ∂    - -        ∂ ∂            ∂ ∂ - + - -     ∂ ∂     ∫ ∫ ∫ ( ) ( ) ( ) ( ) 0 0 0 2 2 2 2
ψ θ ψ θ ρ ω ω       Ω            ∂ ∂ ∂ ∂ - - + - Ω        ∂ ∂ ∂ ∂            ∂ ∂   + - + +       ∂ ∂       ∫ ∫ ∫ ∫ (6)
where E sh and G sh are respectively the Young's modulus and the shear modulus of the shaft (G sh =E sh /(2(1+ν sh )) hence ν sh is the Poisson's ratio). mo sh k and di sh k are respectively the mean shear correction factor of the cross-section and that relative to the section asymmetry of the shaft.

Mass Unbalance.

Let us consider a concentrated mass unbalance m mu placed at a point P mu on the disk (y mu =y d ) with a distance r mu from the shaft geometric center. Its initial angle with the Oz axis of the frame R at rest is η mu . The mass unbalance is characterized by the kinetic energy computed as follows [START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF] ( )

2 g g mu mu T R R mu mu P P m T = v v (7)
The components of the mass unbalance translational velocity vector

g mu R P
v are functions of (ω x ,ω y ,ω z ) and (x O ,y O ,z O ).

Bearing. Figure 3 shows a simple representation of a hydrodynamic bearing which is composed of a fixed journal containing a rotating shaft. The points O and O l =C sep are respectively the bearing center and the shaft geometric center. The radius, length and clearance of the bearing are respectively r be , l be and c be =r be -r sh where r sh is the shaft radius. At a constant speed of rotation Ω of the rotor and for a constant static load W r created by the rotor weight, the shaft geometric center C sep in the bearing holds a static equilibrium position defined by the displacement vector δ be,sep =<u be,sep ,w be,sep > R T expressed in the frame R or equivalently by the eccentricity , be be sep e = δ of the shaft center in the journal and the attitude angle φ be between the W r load direction and the line of centers OC sep .

Figure 3. SCHEMATIC VIEW OF A HYDRODYNAMIC BEARING

In the present study, the short bearing theory is considered (l be /d be ≤1/8 where d be =2r be ) and the static solution can be reached using the following relations obtained from the Reynolds equation with the Gümbel boundary conditions [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF] ( ) ( ) ( )

1 2 2 2 2 2 2 2 2
16 1 1 be be be be be r be be be be be

r l W r l c d + -     = Ω     -     ε ε π ε µ ε (8)
( )

2 1 4 be be be tg - = ε π ϕ ε (9)
where ε be is the relative eccentricity (ε be =e be /c be ) and µ is the fluid film dynamic viscosity. The nonlinear Eq. ( 8) is solved by an iterative Newton-Raphson method and provides the eccentricity ε be and then the static radial displacement e be of the shaft center. The components of the vector δ be,sep are obtained by a classical change of basis. The hydrodynamic fluid forces F be =<F be u ,F be w > R T produced by the bearings and projected in the frame R can be obtained by integration of the fluid film pressure (Reynolds equation in the dynamic regime) over the bearing. In order to apply the Lagrange's equations, the virtual work δW be of these forces has to be established ( )

, T be be be be be

W = F δ δ δ ɺ δ δ ( 10
)
If the transverse dynamic displacements δ be =<u be ,w be > R T of the shaft elastic line are supposed to be small in the vicinity of the static position δ be,sep , the linear analysis of bearings can be applied by constructing a first order Taylor expansion of the fluid film forces

    = =         c k (12)
, ;

be be be sep be be

= - = ∆δ δ δ ∆δ δ ɺ ɺ (13)
where c be and k be are the damping and stiffness matrices of the linearized hydrodynamic bearing (Fig. 4) whose analytical expressions can be found in [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF].

Figure 4. BEARING DAMPING AND STIFFNESS COEFFICIENTS

Equations of Motion

The finite element method is chosen for discretizing the rotor and describing its flexural motion as a function of the nodal displacement vector defined by δ n =<u n ,w n ,θ n ,ψ n > R T because the rotor has two lateral translations and two rotations at each node. The finite element used for the shaft modeling has two nodes and the shape functions are based on the Timoshenko beam theory.

The linear equations of motion of the finite element rotor model are obtained after applying the Lagrange's equations to the energies for the disk, the shaft finite elements and the mass unbalance as well as to the virtual work of the hydrodynamic bearings and assembling appropriately the produced matrices and vectors. They are written with respect to the noninertial reference frame connected to the rotor rigid base R as follows

( ) ( ) ( ) ( ) ( ) , , , r r r r r r r be r sep be r sep t t t t + + = + + M δ C δ K δ F F δ 0 K δ ɺɺ ɺ (14)
with M r (t), C r (t) and K r (t) the mass, damping and stiffness matrices with time-varying parameters due to the asymmetry of the rotating rotor and to its moving base. r δ ɺɺ , r δ ɺ and r δ are respectively the acceleration, velocity and displacement vectors of the on-board rotor-bearin system of dimension 4(n esh +1)×1 hence n esh is the number of shaft finite elements. F r (t) is the external force vector. In addition, the static hydrodynamic force vector F be (δ r,sep ,0) of the bearings is opposite and equal to the rotor weight vector , r W d sh F . Lastly, δ r,sep is the static solution vector of the rotor due to the hydrodynamic bearings.

The matrices of Eq. ( 14) are expressed in what follows t

( ) ( ) ( ) 2
= + + Ω + Ω Ω + Ω Ω + C C C C C C C ω ω (16) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2 2 2 2 , , , , , , 2 , , , , , , , , , 2 , , 2 , , , , , , , , , , , , , , cos 2 sin 2
+ + + + + + K K K K K K K K K K K K K K K ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ɺ ɺ ɺ ɺ ( ) ( ) ( 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 , , , , 2 , , , , , , , , , 2 2 , , , , , , , , , , , , 2 , , , , , , , , , , 2 , , , , cos 2 y
Ω Ω Ω + + + + Ω + + Ω + + + K K K K K K K K K ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ɺ ɺ ) ( ) sin 2 x z t Ω ω ω (17)
The subscripts "d", "sh", "be" and "b" refer to the disk, shaft, bearing as well as base respectively and express the contribution to the phenomenon represented by the corresponding matrix. The superscripts "c 2 " and "s 2 " denote the geometric asymmetry of the rotor expressed in terms of the time-varying trigonometric functions cos(2Ωt) and sin(2Ωt). The superscript "id" stands for the rotor internal damping introduced by estimating the Rayleigh damping coefficients, "g" for the rotor gyroscopic effect, "e" for the shaft elasticity corresponding to the bending and shear deformations, "re" for the rotational effects due to the base rotations (these effects come from the kinetic energies of the disk and the shaft) and "gse" for the geometric stiffening effects corresponding to the centrifugal stress due to the base rotations.

The vector F r (t) is defined as follows x

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 1 2 2 1 1 1 1 2 1 1 , , , , , , , , 2 2 , , , , , 2 2 , , , , , , 2 ,
= + + + + Ω + Ω = - + Ω Ω + Ω Ω + + Ω + + + + Ω F F F F F F F V V V V V V V V V ɺ ɺ ( ) ( ) ( ) ( ( ) ( )) ( ) ( 2 1 1 1 2 2 1 1 1 , , , 2 , , , , , , 2 2 , , , , , 2 2 ,
ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω Ω + + Ω + + + + Ω - + - + + - - - + - + - + + - V V V V V V V V ɺ ɺ ɺ ɺɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺɺ ɺ ( ) ( )) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2 , , , , , , , , , , , , , , , , , , , , , 2 cos 2 2
- Ω + + Ω + - -Ω - Ω + + Ω + Ω + - V V V V V V V V V ɺ ɺ ɺ ɺ ɺ ɺ ɺ ( ) ( ) ( ) ( ) 2 , , , sin 2 2 sin 2
z y z s z x x y d sh b t t θ ω ω ω ω ω ω ω Ω - Ω + + Ω + Ω V ɺ (18)
where

V mu (F mu (t)), V d,sh,b (F d,sh,b (t)) and V mu,b (F mu,b (t)
) are load vectors (force vectors) associated respectively with the mass unbalance, the inertia force due to base motions and that due to coupling between both phenomena. The superscripts "c 1 " and "s 1 " signify the components of the mass unbalance force expressed in terms of the time-varying trigonometric functions cos(Ωt) and sin(Ωt). The superscripts u, w, ψ and θ denote the direction of the action force components associated with the rotor base motions.

The transient dynamic motion of the rotor is then obtained by solving Eq. ( 14) by means of the implicit Newmark timestep integration algorithm based on the average acceleration. The static equilibrium position δ r,sep is used to initialize the transient dynamic problem. The final integration time is chosen such that the transient effects have disappeared and the steadystate regime has been reached.

NUMERICAL SIMULATIONS AND DISCUSSION

The symmetric rotor-rigid/elastic linearized bearing system presented in Fig. 5 is subjected to rotating mass unbalance as well as to base harmonic rotation without/with base harmonic translation: the rotation around the Ox axis is given by ω x =ω x,a cos(Ω x t) in rad/s, while the translation along the Oz axis is expressed as z O =Z O cos(Ω z t) in m.

Figure 5. ON-BOARD ROTOR-BEARING SYSTEM CONFIGURATION

The physical properties as well as the geometry of the rotor and the bearings are listed in Table 1. The shaft is discretized into eight identical 2-node Timoshenko beam finite elements. The disk is located at node 5 and the bearings # 1 and # 2 are respectively placed at nodes 1 and 9. The rotor runs at a speed of rotation Ω=1200 rpm (20 Hz=mass unbalance frequency).

The static equilibrium position of the shaft geometric center in the fluid film bearings is given by δ be,sep =<-5.71×10 -5 ,-1.76×10 -4 > R T m. The bearing damping and stiffness matrices are expressed in what follows

  × × =   × ×     × × =   × ×   c k (19)
The rotor system is supported either by rigid bearings or hydrodynamic ones leading in the case of fixed base either to symmetric damping and stiffness matrices with isotropic diagonal and cross-coupling terms or to nonsymmetric matrices respectively. ɺ in Eq. ( 18)). The two transverse displacements associated with these force components are different and this breaks the asymmetry of the rotor behavior. For the selected speed of rotation of the rotor, the z-amplitude is very higher than the x-amplitude which remains mostly the same. FFT exhibits two frequency components due to the mass unbalance excitation (20 Hz) and to the base harmonic motion (80 Hz). Figure 8 presents the disk orbit due to the mass unbalance effect in the presence of hydrodynamic bearings. The nonsymmetric damping and stiffness coefficients make the bearings anisotropic and the orbit elliptical with diagonal axes defining the phase between the mass unbalance excitation and the rotor response. The orbit center coincides with the point C sep (static position of the shaft center in the bearings). It is noted that the orbit is large compared to that corresponding to rigid bearings because of the combination between the bending modes of the rotor and the rigid body modes relative to the rotor motion in the hydrodynamic bearings. z-FFT indicates that the disk vibration has an amplitude of about 7×10 -6 m and the same mass unbalance frequency (20 Hz).

Table 1. ROTOR AND BEARING DATA

Figure 8. MASS UNBALANCE ORBIT AND ITS z-FFT AT THE DISK FOR FLEXIBLE BEARINGS AND FIXED BASE

Figure 9 displays the disk orbits for hydrodynamic bearings during base rotational motions. The orbit characteristics (shape and magnitude) change with amplitude and frequency of the base harmonic excitation. It should be noted that for very small amplitudes and different frequencies, the orbit shapes become more complicated with respect to those obtained when the base is fixed (see Fig. 9(a)). On the other hand, the base motion amplitudes change the x and z orbit magnitudes. z-FFTs exhibit two frequency components due to the mass unbalance excitation (20 Hz) and to the base harmonic motions (50 Hz or 80 Hz). They also show that the displacement amplitude due to the mass unbalance excitation (20 Hz) does not roughly change for base motions of frequency 50 Hz and 80 Hz, i.e., the influence of the base rotations on the mass unbalance excitation (see Eq. ( 18)) can be neglected at the lower amplitudes of the base excitation.

Figures 10 and11 display the orbits of the disk and their z-FFTs in the presence of mass unbalance forces and of a combination of base rotational and translational motions. In Fig. 10, the base rotation is kept constant (ω x,a =5×10 -2 rad/s and Ω x =50 Hz) while the base translation has two amplitudes and three frequencies. In Fig. 11, the base translation is kept constant (Z O =5×10 -5 and Ω z =50 Hz) while the base rotation has two amplitudes and three frequencies. Increasing the amplitudes and frequencies makes the orbits larger and more complicated especially when compared with the mass unbalance orbit for a fixed base; see Fig. 8. The assumption to have linear behavior of the bearing is questionable. FFTs in Figs. 10 and 11 exhibit three frequency components due to the mass unbalance excitation (20 Hz), to the base harmonic rotation (50 Hz in Fig. 10 and 80 Hz in Fig. 11) as well as to the base harmonic translation (80 Hz in Fig. 10 and 50 Hz in Fig. 11).

C sep

CONCLUSIONS

A finite element model is presented to analyze the dynamic behavior of a symmetric on-board rotor-rigid/hydrodynamic bearing system whose base is subjected to sinusoidal rotations without/with sinusoidal translations.

The rotational effects and the geometric stiffening effects relative to the centrifugal stressing due to the base rotations are taken into account.

The base rotations create time-varying parametric coefficients which can lead to lateral dynamic instability.

In the case of a rotor mounted on rigid bearings and running at the lower speeds of rotation, the base rotation around a transverse axis creates an orbit having its greatest magnitude in the perpendicular transverse direction. Unlike the previous case, the base rotation effects concern the two transverse directions when the rotor is supported by hydrodynamic bearings.

It is noted that the shape and the magnitude of the orbits can be significantly affected by the base motion frequency and amplitude respectively.

In the case of considerable rotation amplitude compared to the speed of rotation of the rotor, it is shown that the mass unbalance forces can depend on the base rotation around a transverse axis.

The frequency components due to the mass unbalance excitation and to the base harmonic motions appear in FFTs of the rotor flexural vibrations.

In the case of large orbits, the assumption of a linearized hydrodynamic bearing model is questionable. A nonlinear model is to be considered and the hydrodynamic bearings have to be treated as external nonlinear forces acting on the shaft within the bearings.

Figure 2 .

 2 Figure 2. EULER ANGLES

I

  are the inertias of the disk mass about the Ox, Oy and Oz axes respectively. In addition, explain the effects due to the mean inertia of the disk mass and those due to the inertia modeling the disk geometric asymmetry. Thus when the disk is asymmetric, the disk center are functions of the components (ω x ,ω y ,ω z ) and the coordinates (x O ,y O ,z O ). The final expression of the kinetic energy of the disk having a mass center positioned at the arbitrary abscissa y d relative to the frame R is then of the following form[START_REF] Duchemin | Dynamic Behavior and Stability of a Rotor Under Base Excitation[END_REF]

  the mean inertia of the cross-section and the inertia characterizing the asymmetry of the shaft. The rigid base motion relative to the ground has no influence on the strain energy of the shaft because this energy depends only on the stresses and therefore on the transverse deflection of the shaft with respect to the rotor base R. In addition to the bending deformation, the shear effects highlighted by Timoshenko and the geometric stiffening effects corresponding to the centrifugal stressing due to the base rotational motions are taken into account. The strain energy of a shaft is defined by

Figure 6

 6 Figure6shows the classical disk orbit due to the mass unbalance for rigid bearings and a fixed base. Since the matrices of the rotor-bearing system are symmetric and skewsymmetric with isotropic diagonal and cross-coupling components, the dynamic behavior is symmetric and the orbit is circular. Its center coincides with the point O (bearing center). FFT of the orbit shows that the rotor displacement has an amplitude almost of order 1×10 -7 m and the same mass unbalance frequency (20 Hz).

Figure 6 .Figure 7

 67 Figure 6. MASS UNBALANCE ORBIT AND ITS z-FFT AT THE DISK FOR RIGID BEARINGS AND FIXED BASE

Figure 7 .

 7 Figure 7. DISK ORBIT AND ITS z-FFT FOR RIGID BEARINGS AND BASE HARMONIC ROTATION: ω x,a =5×10 -2 rad/s AND Ω x =80 Hz

Figure 9 .Figure 10 .Figure 11 .

 91011 Figure 9. DISK ORBITS AND THEIR z-FFTs FOR FLEXIBLE BEARINGS AND TWO AMPLITUDES OF BASE HARMONIC ROTATIONS: ω x,a = (a) 10 -2 rad/s, (b) 5×10 -2 rad/s, (Ω x =20 Hz, 50 Hz AND 80 Hz)

  [START_REF] Lalanne | Rotordynamics Prediction in Engineering[END_REF]

	r	t	be	id d sh	g d sh	g c d sh	t
			g s d sh		t	y d sh b re	y


Marc Berveiller
email: marc.berveiller@edf.fr

Bruno Sudret
email: bruno.sudret@edf.fr

Maurice Lemaire
email: maurice.lemaire@ifma.fr

Stochastic finite element : a non intrusive approach by regression

Keywords: stochastic finite element, non intrusive method, polynomial chaos, stochastic response surface method MOTS-CLÉS : Eléments finis stochastiques, méthode non intrusive, chaos polynomial, surfaces de réponse stochastiques

come

Introduction

The Stochastic Finite Element Method (SFEM) developed by [START_REF] Ghanem | Stochastic finite elements -A spectral approach[END_REF] allows to solve stochastic boundary value problems involving spatially randomly varying material properties usually described as Gaussian or lognormal random fields. The method is based on the discretization of the input random fields and the expansion of the mechanical response onto the so-called polynomial chaos. A similar procedure allowing to model random material properties and loads by means of an arbitrary number of random variables of any type has been recently proposed [START_REF] Sudret | Eléments finis stochastiques en élasticité linéaire[END_REF], Sudret et al., 2005[START_REF] Berveiller | Eléments finis stochastiques : approches intrusive et non intrusive pour des analyses de fiabilité[END_REF]. In both cases, the coefficients of the response expansion are computed using a Galerkin procedure, which leads to a linear system whose size is equal to the number of degrees of freedom of the finite element model multiplied by the number of coefficients retained in the response expansion. The output of this method is the expansion of response quantities S(θ) (such as displacement, strain or stress components) onto the polynomial chaos.

S(θ) =

P -1 j=0 s j Ψ j ({ξ k (θ)} M k=1 ) [1]
where s j are the coefficients to be determined, {Ψ 0 , • • • , Ψ P -1 } are multidimensionnal Hermite polynomials of order less or equal than p, {ξ k (θ)} M k=1 are standard Gaussian random variables, M is the number of random input parameters and P denotes the size of the polynomial chaos basis, which is given by:

P = (M + p)! M ! p! [2]
In this paper, a non intrusive method (called regression method) is proposed for computing these coefficients using a series of deterministic analysis.

Non intrusive regression method

The non intrusive regression method presented in this section is based on a leastsquares minimization between the exact solution and the approximate solution based on the polynomial chaos [START_REF] Isukapalli | Uncertainty Analysis of Transport-Transformation Models[END_REF][START_REF] Mahadevan | A stochastic response surface method for random field problems[END_REF], Berveiller et al., 2004[START_REF] Berveiller | Eléments finis stochastiques : approches intrusive et non intrusive pour des analyses de fiabilité[END_REF]. First the input random variables (gathered in a random vector X whose joint PDF (Probability Density Function) is prescribed) are transformed into a standard Gaussian vector ξ. If these M variables are independent, the one-to-one mapping reads :

ξ i = Φ -1 (F i (X i )) [3]
where Φ is the standard Gaussian CDF (Cumulative Density Function) and

{F i (X i ), i = 1, • • • , M } are the marginal CDF of the X i 's. Let us denote by {ξ (k) , k = 1, • • •
, n} n outcomes of the standard Gaussian random vector ξ. For each outcome ξ (k) , the isoprobabilistic transform yields a vector of input random variables X (k) (Eq.( 3)). Using a classical finite element code, the response vector S (k) can be computed. According to Eq.( 1), each scalar component of the response S (displacement, strain, stress, etc.) may be approximated as follows:

S ≈ S(ξ) = P -1 j=0 s j Ψ j (ξ) [4]
where {s j } P -1 j=0 are coefficients to be computed. The regression method consists in finding the set of coefficients that minimizes the difference:

ΔS = n k=1 S (k) -S(ξ (k) ) 2 [5]
These coefficients are solution of the following linear system:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ n k=1 Ψ 0 (ξ (k) )Ψ 0 (ξ (k) ) • • • n k=1 Ψ 0 (ξ (k) )Ψ P -1 (ξ (k) ) . . . . . . . . . n k=1 Ψ P -1 (ξ (k) )Ψ 0 (ξ (k) ) • • • n k=1 Ψ P -1 (ξ (k) )Ψ P -1 (ξ (k) ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ • ⎛ ⎜ ⎝ s 0 . . . s P -1 ⎞ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ n k=1 S (k) Ψ 0 (ξ (k) ) . . . n k=1 S (k) Ψ P -1 (ξ (k) ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ [6]
Note that the P × P matrix on the left hand size is evaluated once and for all. Moreover it is independent on the mechanical problem under consideration. Then the coefficients of the expansion of each response quantity s i are obtained by the resolution of the linear system Eq.( 6). The crucial point in this approach is to properly select the collocation points, i.e. the outcomes {ξ (k) , k = 1, • • • , n}. Note that n ≥ P is required so that a solution of (6) exist. [START_REF] Webster | Application of the Probabilistic Collocation Method for an Uncertainty Analysis of a Simple Ocean Model[END_REF] and [START_REF] Isukapalli | Uncertainty Analysis of Transport-Transformation Models[END_REF] choose for each input variable the (p + 1) roots of the (p + 1)-th order Hermite polynomial H p+1 (x), and then built (p + 1) M vectors of length M using all possible combinations between theses different roots (Figure 1). Then they select n outcomes {ξ (k) , k = 1, • • • , n} out of these (p + 1) M possible combinations as follows :

- [START_REF] Webster | Application of the Probabilistic Collocation Method for an Uncertainty Analysis of a Simple Ocean Model[END_REF] selected n = P + 1 and the (P + 1) outcomes which

maximize ϕ n (ξ (k) ) = (2π) -M/2 exp -1/2 ξ (k) 2
. This corresponds to selecting the points that minimize ξ 2 .

- [START_REF] Isukapalli | Uncertainty Analysis of Transport-Transformation Models[END_REF] selected by the same method 2(P + 1) outcomes, adding the null vector ξ = 0 if it is not already included in the experimental design.

Figure 1. Scheme of the selection of the collocation points (

M = 2, p = 4, n = 9)
No indication of the reasons of these choices for n could be found in the literature. In the sequel, a parametric study is carried out to investigate the influence of n on the accuracy of the results. Indeed, using the selection method mentioned above, it is easy to increase the accuracy by taking more collocation points, without losing the benefit of previous computations. It is also possible to define an error estimator for any response quantity S. Suppose that we have n = kP collocation points, where k = 1, 2, • • • , and that the minimal k for a given accuracy is looked after. For each k ≥ 1, an estimation of the response quantity S (assumed not vanishing) is computed together with a convergence rate ε k :

ε k = S k+1 -S k S k [7]
If this rate is smaller than a tolerance, say 1%, the regression scheme is said to be accurate. Otherwise P new collocation points are added to the scheme (leading to P additional finite element runs). The procedure is followed up until convergence.

Application example : settlement of a foundation

Description

Let us consider an elastic soil mass made of two layers of different isotropic linear elastic materials lying on a rigid substratum (Figure 2(a)). A foundation on this soil mass is modeled as two uniform pressure loads P 1 and P 2 applied respectively over a length 2B 1 and 2B 2 of the free surface (Figure 2(b)). Due to the symmetry, half of the structure is modeled by finite element. The mesh comprises 99 nodes and 80 4-node elements which allows a 1%-accurate evaluation of the maximal settlement compared to a reference solution. In the sequel, three different analysEs of this problem are performed. First, the statistical moments of the maximum settlement are computed. Then a reliability analysis with respect to admissible settlements is carried out. Finally, the approach is used in a purely deterministic context, where the parameters are given interval values.

6 ? 6 ? A - 2B 1 = 10m t 1 = 7, 75 m t 2 = 22, 25 m E 1 , ν 1 E 2 , ν 2 P 2 ???? P 1 - 2B 2 = 5m

Statistical moments of the response

In this section, the six parameters (namely Young modulus and Poisson ration in each layer and two loads) are modelled as random variables described in Table1. We consider a polynomial chaos of order three. Thus P = (6+3)! 6!3! = 84 coefficients are to be determined using the non intrusive method. We are interested in the evaluation of the maximum settlement of the foundation obtained at point A and denoted by u A .

The four statistical moments of the maximal settlement of the foundation obtained by the non intrusive method are compared to reference results. From Eq.( 1), the statistical moments of any response quantity can be computed analytically. The mean and the variance of S are given by:

E[S] = s 0 [8]
Var[S] = σ 2 S = P -1 i=1 E[Ψ 2 i ]s 2 i [9]
The skewness and the kurtosis coefficients are given by:

δ S = 1 σ 3 S P -1 i=1 P -1 j=1 P -1 k=1 E[Ψ i Ψ j Ψ k ]s i s j s k [10] κ S = 1 σ 4 S P -1 i=1 P -1 j=1 P -1 k=1 P -1 l=1 E[Ψ i Ψ j Ψ k Ψ l ]s i s j s k s l [11]
The analytical formulae to compute the expectation of products of Ψ j polynomials are given in [START_REF] Berveiller | Eléments finis stochastiques : approches intrusive et non intrusive pour des analyses de fiabilité[END_REF]. Reference results are obtained by Monte Carlo simulation (10000 samples) using a direct coupling between the finite element code Code_Aster1 and the probabilistic code PROBAN. Figure 3 shows the evolution of the statistical moments versus the number of samples in the Monte Carlo simulation. Table 2 gathers the reference results and those obtained by the non intrusive method.

With the non intrusive regression method, good results are obtained for the four first moments if more than 420 points (i.e. (M -1)P points) are taken. Figure 4 presents the evolution of the four first statistical moments (divided by the converged value in the regression method) versus the number of collocation points. We can note that from 420 points (i.e. 5 × 84 ≡ (M -1)P points), results have converged. Note that the skewness and kurtosis coefficients obtained by Monte Carlo simulation slightly differ from the regression results.

Reliability analysis

In reliability analysis, the failure criterion of a structure is defined in terms of a limit state function g(X, S(X)) which may depend both on basic random variables X and response quantities S(X). The domain defined by g(X, S(X)) < 0 is the failure domain and g(X, S(X)) > 0 defines the safe domain. The boundary between these two domains is the limit state surface. When performing a stochastic finite element analysis, the true limit state function is replaced by its polynomial approximation onto the polynomial chaos. Thus the approximation is defined by definition in the so-called standard normal space, which is the space used for the First Order Reliability Method (FORM) in reliability analysis:

g(X, S(X)) ≡ g ⎛ ⎝ {ξ k } M k=1 , P -1 j=0 S j Ψ j ({ξ k } M k=1 ) ⎞ ⎠ [12]

Figure 4. Evolution of moments vs. the number of collocation points in the non intrusive regression method

Thus the reliability problem, which is already formulated in the standard Gaussian space, may be solved by any available method including Monte-Carlo Simulation, FORM/SORM, Importance Sampling (IS), etc. [START_REF] Ditlevsen | Structural reliability methods[END_REF].

The reliability of the foundation with respect to the maximum admissible settlement ū is investigated. Hence the limit state function can be written:

g(X) = ū -u A (E 1 , E 2 , ν 1 , ν 2 , P 1 , P 2 ) [13]
where u A (E 1 , E 2 , ν 1 , ν 2 , P 1 , P 2 ) is the maximal displacement. The failure probability is computed by FORM and importance sampling for different values of ū using the two strategies mentioned above, namely :

-a direct coupling between the finite element code Code_Aster and the probabilistic code PROBAN by using FORM and importance sampling. Eq.( 13) is used for this purpose, when u A (•) is obtained with the finite element code. 1,000 samples are used in IS allowing a coefficient of variation of the simulation of 1%.

-the non intrusive regression method followed by a FORM reliability analysis and importance sampling (1,000 samples and a coefficient of variation of 1% for the simulation). Note that FORM as well as IS are performed using analytical limite state function in this case :

g(ξ) = ū - 83 j=0 u Aj Ψ j ({ξ k } 6 k=1 ) [14]
In this approach, the impact of the number of collocation points used to compute the u Aj 's onto the accuracy is studied.

Table 3 presents the probability of failure obtained by direct coupling and by regression using various numbers of collocation points. Figure 5 shows the evolution of the ratio between the logarithm of the probability of failure (divided by the logarithm of the converged probability of failure) versus the number of collocation points for several values of the maximum admissible settlement ū. Again, accurate results are obtained when using 420 collocation points or more for different values of the failure probability (from 10 -1 to 10 -4 ). When taking less than 420 points, results are inaccurate. When taking more than 420 points, the accuracy is not improved. It seems that the best compromise between accuracy and efficiency is (M -1)P collocation points in this case, as it was also observed in the previous section. 3.09.10 -1 1.62.10 -1 2.71.10 -1 3.31.10 -1 3.23.10 -1 3.32.10 -1 15 6.83.10 -2 6.77.10 -2 6.90.10 -2 8.43.10 -2 6.73.10 -2 6.93.10 -2 20 2.13.10 -3 -9.95.10 -5 8.22.10 -4 2.01.10 -3 1.98.10 -3 22 4.61.10 -4 -7.47.10 -7 1.31.10 -4 3.80.10 -4 4.24.10 -4 0.9

Deterministic parametric analysis

In this section, the formalism of stochastic finite element is used to build a deterministic response surface of the 6-parameter model. In this sense, the input parameters are given a range of variation [a; b]. To be able to apply the proposed method, these intervals are then considered as support of uniform random variables (Table 4). (a) shows the relative error between a deterministic finite element analysis and the approximation of the settlement u A by a third order polynomial chaos using 420 collocation points. We can note that the maximum error is about 4%. Figure 6(b) shows the Henry's line [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF] of the residual which allows to determine whether the residual follows a Gaussian distribution. The straight line (dotted) draw on this figure is the fitted line to the Henry's line. We can also determine the mean and the standard deviation of the residual fitted to a Gaussian distribution and compared them with those obtained empirically. The mean of the residual is equal to -3.96.10 -6 m (the mean settlement is -0.082m). The standard deviation is equal to 7.94.10 -4 . The skewness coefficient and the kurtosis coefficient are respectively equal to -0.56 and 0.64. For a Gaussian distribution these values are respectively 0 and 3. Values of skewness and kurtosis coefficients allow to conclude that the residual is not fitted to a Gaussian distribution. The Table 5 gathers the R 2 coefficient [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF] of the regression versus the number of collocation points. Again this coefficient shows an excellent quality of the regression whatever the number of points. It is concluded that P collocation points are sufficient to build a deterministic response surface of the foundation. Note that the deterministic response surface built using the polynomial chaos is now straight forward to use. For any values of the input parameters (in the range given in Table 4), the corresponding standard normal values are computed using Eq.(3). Then the quantities are used to evaluate the polynomial expansion.

Conclusion

The paper presents a non intrusive method to compute coefficients of the polynomial series expansion of the response in stochastic finite element analysis. This method is based on a least-squares minimization of the distance between the exact solution and the polynomial expansion. The result is interpreted as a stochastic polynomial response surface. The collocation points are selected among M -uplets of roots of Hermite polynomials. The method is illustrated by the analysis of the settlement of a foundation. Various analysis are performed such as computation of the statistical moments of the response and reliability analysis. Based on the selection scheme, it appears that a number of points equal to (M -1) times the size of the polynomial chaos (i.e. 420 points in the application example) provides excellent accuracy both in the mean region (moments of the response) and in the tail (probability of exceedance of a threshold). As far as a deterministic parametric study is concerned, P collocation points are sufficient to get an accurate deterministic response surface. The great advantage of the non intrusive approach compared to the classical Galerkin approach is that only deterministic finite element models are used. Hence the full non linear capabilities of the code may be used without additional implementation (see for instance [START_REF] Berveiller | Non linear non intrusive stochastic finite element method -Application to a fracture mechanics problem[END_REF]). Note that if several output quantities are of interest, the marginal cost to estimate the response coefficients is low. Indeed the deterministic finite element analysis are usually the expensive part of the evaluation.

Figure 2 .

 2 Figure 2. (a) : Scheme of the foundation -(b) : Mesh of the foundation

Figure 5 .

 5 Figure 5. Evolution of the logarithm of the failure probability divided by the converged value versus the number of collocation points

Figure 6 .

 6 Figure 6. (a) : Relative error (%) -(b) : Henry's line (420 collocation points)

Table 1 .

 1 Definition of the input random variables for moments and reliability analysis

	Parameter	Notation	Distribution	Mean	Coeff. Var.
	Young modulus upper layer	E 1	Lognormal	50 MPa	20%
	Young modulus lower layer	E 2	Lognormal	70 MPa	20%
	Poisson ratio upper layer	ν 1	Uniform [0.28; 0.32]	0.3	4%
	Poisson ratio lower layer	ν 2	Uniform [0.28; 0.32]	0.3	4%
	Load #1	P 1	Weibull	0.2 MPa	30%
	Load #2	P 2	Lognormal	0.4 MPa	20%

Table 2 .

 2 Statistical moments of the maximal settlement

		MC Simulation		Non intrusive method
	#FE runs	10000	84	168	336	420	4096
	Mean	-1.097	-0.566	-0.735	-1.040	-1.108	-1.108
	Std. Dev.	0.0256	0.0802 0.0747	0.0439	0.0250 0.0250
	Skewness	-0.6019	1.6664 1.9638	2.3967 -0.5182 -0.5303
	Kurtosis	3.7214	8.7149 10.9272 23.7676 3.4826 3.5095

Table 3 .

 3 Failure probability P f

	ū	Direct			Non intrusive
	(cm) Coupling	84	168	336	420	4096
	12

Table 4 .

 4 Definition of input parameters for the parametric analysis

	Parameter	Notation	Interval
	Young's modulus upper layer	E 1	[50 ; 70] MPa
	Young's modulus lower layer	E 2	[70 ; 90] MPa
	Poisson's ratio upper layer	ν 1	[0.1 ; 0.4]
	Poisson's ratio lower layer	ν 2	[0.1 ; 0.4]
	Load #1	P 1	[0.1 ; 0.3] MPa
	Load #2	P 2	[0.1 ; 0.3] MPa
	Figure 6

Table 5 .

 5 R 2 coefficient

	Number of collocation points	84	168	252	336	420
	R 2	0.999 0.992 0.993 0.993 0.993

http://www.code-aster.org