Nicolas Davister

Fabrice Frebel

Christophe Geuzaine

Comparison of Finite Element Formulations for the Modelling of High-Frequency Effects in Electromagnetic Coils

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The accurate prediction of coupled inductive and capacitive effects in electromagnetic coils is becoming of crucial importance in several industrial applications, either due to the increase of operating frequency or to the increase of voltage levels. In this paper we review and compare 2D and 3D finite element formulations of electromagnetic coils able to predict both inductive and capacitive effects, in order to design approximate 2D models whose computational complexity is compatible with current industrial design practice.

Index Terms-Finite element method, electromagnetic coils, high-frequency effects.

I. INTRODUCTION

H IGH frequency effects, i.e. the coupling of inductive and capacitive effects, could long be ignored in the design of electromagnetic coils for electrotechnical applications. With the advent of extremely fast switching electronic components based on silicon carbide or gallium arsenide, this coupling can however not be neglected anymore when designing the magnetic components (inductors and transformers) of either extremely fast-switching, or very high-voltage, power electronic circuits. Resonances due to parasitic capacitances cannot be neglected at the considered frequencies, in particular for electromagnetic compatibility issues.

While the design of such magnetic components has traditionally been (and for the most part still is) carried out analytically, the last decade has increasingly seen the adoption in engineering offices of finite element software tools. To keep the computational cost acceptable, only 2D models are commonly used. Moreover, these models usually consider decoupled magnetic and electric phenomena [START_REF] Biela | Using transformer parasitics for resonant converters -a review of the calculation of the stray capacitances of the transformers[END_REF][START_REF] Asensi | Modeling high-frequency multiwinding magnetic components using finite-element analysis[END_REF].

This paper follows up on the comparison carried out in [START_REF] De Grève | Finite element models for studying the capacitive behaviour of wound components[END_REF] for a multi-turn coreless inductor, which in addition to fully decoupled magnetic and electric formulations considered a weak (one-way) 3D coupling [START_REF] Dular | Dual finite element formulations for the 3D modeling of both inductive and capacitive effects in massive inductors[END_REF][START_REF] Greve | Numerical modeling of capacitive effects in HF multiwinding transformers -Part I: A rigorous formalism based on the electrostatic equations[END_REF]. Here we compare both decoupled and fully-coupled models, and propose 2D approximations of the fully coupled model that aim at bringing down the computational cost to a level acceptable for current industrial design practice.

II. FINITE ELEMENT FORMULATIONS FOR ELECTROMAGNETIC COILS

The coupling level between inductive and capacitive effects occurring in an electromagnetic coil depends on the ratio between its length and the excitation frequency. For this reason, three main approaches are possible, which consider decoupled, weakly coupled or fully coupled formulations [START_REF] De Grève | Finite element models for studying the capacitive behaviour of wound components[END_REF]. Fig. 1 and 2 illustrate the main characteristics of these methods and the results they provide on a simple transformer made of two circular coils, respectively. In the simplest approach, in which the electric and magnetic fields are considered to be totally decoupled, a magnetodynamic formulation is used to evaluate the resistive and inductive effects whereas the capacitive effects are dealt using an electrostatic formulation [START_REF] Biela | Using transformer parasitics for resonant converters -a review of the calculation of the stray capacitances of the transformers[END_REF][START_REF] Asensi | Modeling high-frequency multiwinding magnetic components using finite-element analysis[END_REF]. Note that the source of the latter can be either a linear distribution of the scalar electric potential along the winding (2D case) [START_REF] Greve | Numerical modeling of capacitive effects in HF multiwinding transformers -Part I: A rigorous formalism based on the electrostatic equations[END_REF] or come from a previous solution of the considered problem using an electrokinetic formulation (3D case).

A weak coupling between the magnetodynamic and electrostatic formulations can also be considered, where the source term of the electrostatic formulation originates from the electric field obtained by the magnetodynamic computation inside the conductors [START_REF] Dular | Dual finite element formulations for the 3D modeling of both inductive and capacitive effects in massive inductors[END_REF]. This has the advantage to better capture the first resonance, but is only applicable in 3D [START_REF] De Grève | Finite element models for studying the capacitive behaviour of wound components[END_REF].

Fully coupled formulations lead to the classical full-wave formulation of Maxwell's equations.

Those three approaches are represented in terms of coupling level and computing cost in Fig. 1. Decoupled and full-wave formulations are compared in Fig. 2. These preliminary results illustrate the numerous high frequency resonances not captured by the decoupled model.

III. TOWARDS CHEAPER FULLY-COUPLED FORMULATIONS

While the full-wave formulation takes all physical effects into account, it has two main drawbacks. First, as is well known, the classical full-wave formulation is unstable in the low frequency regime, which is however crucial for power electronics applications where the whole spectrum (from DC to high frequencies) is of practical interest. This loss of stability can be overcome thanks to various stabilization strategies proposed over the last decade [START_REF] Hiptmair | A Robust Maxwell Formulation for All Frequencies[END_REF][START_REF] Eller | A Symmetric Low-Frequency Stable Broadband Maxwell Formulation for Industrial Applications[END_REF].

Second, the full-wave formulation (like the weakly-coupled formulation) is intrinsically 3D, and thus leads to a computational cost that is incompatible with industrial R&D. To overcome this second drawback, we propose a family of 2D fullwave formulations, which do not assume any weak coupling, and can be derived as follows. Assuming for simplicity a Cartesian coordinate system (x, y, z) and a geometrical invariance along z (axisymmetric configurations could be considered as well), we suppose that the electric and magnetic fields are of the form E(x, y, z, t) = (e(x, y)e -i(ωt-iγz) ),

H(x, y, z, t) = (h(x, y)e -i(ωt-iγz) ),
where ω is the pulsation and γ a propagation constant along the invariance direction. Faraday's and Maxwell-Ampère's equations can thus be written respectively as curl γ e = iωµh and curl γ h = (σ -iω )e, where curl γ • := curl ( • e iγz )e -iγz and where µ, σ and are the magnetic permeability, electric conductivity and permittivity, respectively. Developing e(x, y) in its transverse and longitudinal components e t (x, y) and e z (x, y), such that e := e t + e z 1 z , we get curl γ e = curl e t + (grad e z -iγe t ) × 1 z , which leads to the following weak formulation [START_REF] Guenneau | A finite element formulation for spectral problems in optical fibers[END_REF]: find e t and e z satisfying appropriate boundary conditions such that (curl e t , curl ē t ) + (grad e z , grad ē z ) -(iγe t , grad ē z ) + (iγgrad e z , ē t ) + (γ 2 e t , ē t ) = ((iωµσ + w 2 µ )e t , ē t ) + ((iωµσ + w 2 µ )e z , ē z ) holds for appropriate test functions e t and e z , where (•, •) represents the integral over the domain of study of the scalar product of its arguments. This 2D formulation involves both components of the electric field in and perpendicular to the plane of study and is fully coupled, and can be readily discretized using 2D edge finite elements for e t and nodal elements for e z . For small coils compared to the wavelength, it decouples and leads to two uncoupled wave systems.

The full paper will detail this formulation and analyse its performance on test-cases (like the one from Figure 2) for which it does not introduce any approximation compared to the 3D formulation. We will then investigate how, in the case of multi-turn windings, several slices can be coupled through circuit equations, leading to an alternative to more classical transmission line matrix coil modelling techniques [START_REF] De Gersem | Transmission-line modelling of wave propagation effects in machine windings[END_REF].

Fig. 1 .

 1 Fig. 1. Coupling level of the three approaches vs computing cost

Fig. 2 .

 2 Fig. 2. Bode plots of the transformer impedance with full-wave and decoupled formulations

  Comparison of Finite Element Formulations for the Modelling of High-Frequency Effects in Electromagnetic Coils N. Davister 1 , F. Frebel 1 , C. Geuzaine 1 and G. Parent 2 1 University of Liège, Institut Montefiore B28, Dept. of Electrical Engineering and Computer Science, B-4000 Liège, Belgium 2 Univ. Artois, EA 4025, Laboratoire Systèmes Électrotechniques et Environnement (LSEE), F-62400 Béthune, France


Jean-Jacques Sinou
email: jean-jacques.sinou@ec-lyon.fr

Non-linear dynamics and contacts of an unbalanced flexible rotor supported on ball bearings

This study deals with the non-linear dynamic response of a flexible rotor supported by ball bearings. The excitation is due to unbalance force. The finite element rotor system is composed of a shaft with one disk, two flexible bearing supports and a ball-bearing element where the non-linearities are due to both the radial clearance and the Herztian contact between races and rolling elements. A numerical analysis is performed to analyze the non-linear behavior of this bearing rotor by using the Harmonic Balance Method with appropriate condensation located only on the non-linear coordinates of the system in order to minimize computer time. The condensation process reduces the original non-linear rotor system by focusing only on the solution of the non-linear equations of the Fourier coefficients associated with the system's non-linear components. In this study, the procedure is developed for the estimation of the harmonic and super-harmonic responses of the complex rotor system. Consequently, the non-linear unbalance responses and the associated orbits of the bearing rotor will be investigated. Moreover, the transition from contact to no-contact states between rolling elements and races, and the associated restoring contact forces are calculated for different speeds of the unbalanced rotor. Finally, hardening-type nonlinearity or softening-type nonlinearity due to the effects of radial clearance and unbalance mass are examined.

Introduction

One of the most important non-linear mechanical components of rotating machines to be taken into account are bearings due to their considerable influence on the dynamic behavior of rotor systems. For example, various works have studied non-linear behavior due to active magnetic bearings [START_REF] Chinta | Stability and bifurcation of rotor motion in a magnetic bearing[END_REF][START_REF] Virgin | Nonlinear behavior of a magnetic bearing system[END_REF][START_REF] Ji | Non-linear oscillations of a rotor-magnetic bearing system under superharmonic resonance conditions[END_REF] and fluid-film bearing [START_REF] Erich | Handbook of Rotordynamics[END_REF][START_REF] Childs | Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis[END_REF][START_REF] Yamamoto | Linear and Nonlinear Rotordynamics: a Modern Treatment with Applications[END_REF][START_REF] Ding | Non-stationary processes of rotor/bearing system in bifurcations[END_REF][START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF]. The third major type of bearings are the non-linear rolling-bearing elements currently used to support gas turbine engine rotors due to their durability and low power requirements [START_REF] Erich | Handbook of Rotordynamics[END_REF][START_REF] Yamamoto | Linear and Nonlinear Rotordynamics: a Modern Treatment with Applications[END_REF][START_REF] Harris | Rolling Bearing Analysis[END_REF]. It is well known that ball bearings are a source of nonlinearities for rotor systems. The nonlinear behavior of rotor systems involving bearing clearances and Hertzian ball-race contact have been studied by several investigators. For example, the dynamic analysis of rotors with bearing clearance was studied by Ehrich [START_REF] Ehrich | stator whirl rotors in bearing clearance[END_REF][START_REF] Ehrich | High order subharmonic response if high speed rotors in bearing clearance[END_REF]. Sunnersjo [START_REF] Sunnersjo | Ball passage vibrations of rolling bearings[END_REF] studied ball passage vibrations theo-1 retically, taking into account inertia and damping forces. Comparisons with experimentally obtained results were made. Choi and Noah [START_REF] Choi | Nonlinear steady state response of a rotor-support system[END_REF] used the harmonic balance method with the discrete Fourier transform procedure to analyze the non-linear behavior of a rotor bearing system. Tiwari et al. [START_REF] Tiwari | Effect of a ball bearing on the dynamics of a balanced horizontal rotor[END_REF] studied the non-linear dynamics of a horizontal rigid rotor: the appearance of instability and chaos in the dynamic response are observed both numerically and experimentally. They also studied the effect of radial internal clearance of a ball bearing and the appearance of regions of periodic, subharmonic and chaotic behavior [START_REF] Tiwari | Dynamics response of an unbalanced rotor supported on ball bearings[END_REF]. Later, Harsha et al. [START_REF] Harsha | Non-linear dynamic behaviors of rolling element bearings due to surface waviness[END_REF] proposed to take into account various sources of nonlinearity such as Hertzian contact force, surface waviness and internal radial clearance resulting in the transition from contact to no contact state between races and rolling elements. Periodic, quasiperiodic and chaotic behavior were analyzed in detail. Moreover, Harsha [START_REF] Harsha | Nonlinear dynamic response of a balanced rotor supported by rolling element bearings due to radial internal clearance effect[END_REF] investigated the effect of radial internal clearance for rotor bearing systems in which rolling element bearings show periodic, quasi-periodic and chaotic behavior by considering time response, Poincaré maps and power spectra. Mevel and Guyader [START_REF] Mevel | Experiments on routes to chaos in ball bearings[END_REF] observed experimentally the existence of two different routes to chaos in ball bearing dynamic motion. Nataraj and Harsha [START_REF] Nataraj | The effect of bearing cage run-out on the nonlinear dynamics of a rotation shaft[END_REF] investigated the nonlinear dynamic behavior of an unbalanced rotor-bearing system due to cage run-out. They indicated that the responses can be placed in three categories: periodic with no sensitivity to initial conditions or small variations of system parameters, quasi-periodic and chaotic with extreme sensitivities to both the initial conditions and small variations in the system parameters. Jang et al. [START_REF] Jang | Vibration analysis of a rotation system due to the effect of ball bearing waviness[END_REF] studied the vibration due to ball bearing waviness in a rotating system, taking account of the centrifugal force and gyroscopic moment of the ball.

These studies demonstrated that the nonlinear dynamic response of a balanced rotor supported by rolling element bearings can be very complex and both time consuming and costly to perform when parametric design studies are needed. Therefore, due to the fact that the non linear behavior of rotor systems can be complex, much work has been done on the there has been on the treatment of nonlinear differential equations. More particularly, non-linear methods have been developed to analyse vibration problems [START_REF] Nayfeh | Nonlinear oscillations[END_REF][START_REF] Sinou | A non-linear study of a cracked rotor[END_REF]. One of the classical approaches for obtaining the non-linear response of systems is the numerical integration procedure. However, this approach for non-linear models with many degrees of freedom or strong non-linearities can be rather expensive and requires considerable resources both in terms of computation time and data storage. The well-known approximation technique requires an initial assumption about the form of the solution of the non-linear system, i.e. the non-linear solution is approximated by finite Fourier series. In order to study the non-linear dynamics problems of systems with strong non-linearities, the numerical tools most often used are Harmonic Balance Methods and continuation schemes [START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF][START_REF] Cameron | An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems[END_REF]. Moreover, if this system can be considered as a linear structure with few additional non-linear elements, it may be of great interest to keep only the non-linear degrees of freedom, by introducing a condensation process to solve only the nonlinear equations associated with the non-linear components of the system [START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF].

In the present analysis, we propose to demonstrate the efficiency of the Harmonic Balance Method with a condensation process on the non-linear terms in order to study the unbalance responses of a non-linear flexible rotor system with a radial clearance and a Herztian contact between the races and the rolling elements. Secondly, we investigate in particular the evolution of the restoring nonlinear forces on each ball bearing for different rotating speeds. Th effects of the radial clearances and unbalance are investigated and both hardening-type nonlinearity and softening-type nonlinearity are examined.

The paper is divided into three parts. Firstly, the basic concept of the Harmonic Balance Method, the path following continuation based on Lagrangian extrapolation, and the condensation on the non-linear degrees-of-freedom are introduced. Secondly, the non-linear rotor system under study is described and developed. The third part presents non-linear analysis of the rotor-bearing system. The evolution of the non-linear responses within the speed range of interest is investigated for different unbalance masses and radial clearances. The associated orbits of the rotor and the stator are calculated. Attention is finally focused on the contact at the bearing-elements: the non-linear contacts for each ball of the bearing are evaluated and discussed for different rotating speeds of the rotor system and for different radial clearances and unbalance masses.

General theory of the harmonic balance method and condensation process

In this section, the harmonic balance method with a condensation process on the non-linear degrees of freedom will be presented.

Non-linear method: the Harmonic Balance Method

For a flexible non-linear rotor, the equations of motion may be written as

M Ẍ + D Ẋ + KX = F L X, Ẋ, ω, t + F N L X, Ẋ, ω, t = F X, Ẋ, ω, t (1)
where F L and F N L are the linear and non-linear terms of the rotor system. In order to estimate the response of the non-linear system as a truncated Fourier series (if this solution exists), the right-hand side of the system is assumed to be a function that is periodic in time with period T . Thus we assume that the non-linear dynamical response of the rotor may be approximated by finite Fourier series with ω = 2π T the fundamental frequency:

X (t) = B 0 + m k=1 (B k cos (kωt) + A k sin (kωt)) ( 2
)
where m is the order of the Fourier series. B 0 , A k and B k define the unknown coefficients of the finite Fourier series. It can be seen that the Fourier series considered in this study are developed for harmonic and super-harmonic responses of the non-linear rotor system. The number of harmonic coefficients is selected on the basis of the number of significant harmonics expected in the non-linear dynamical response. Generally speaking, harmonic components become less significant when m increases. Moreover, we assume that the vector force F X, Ẋ, ω, t can be solved in finite Fourier series of order m

F X, Ẋ, ω, t = C 0 + m k=1 (C k cos (kωt) + S k sin (kωt)) (3)
Substituting equations 2 and 3 in equation 1 yields a set of (2m + 1) * n equations (where n is the number of the degrees-of-freedom for the complete rotor bearing system).

Using the first n th equations, the constant terms B 0 can be determined. We obtain

KB 0 = C 0 (4)
The 2m * n remaining equations defining the k th Fourier coefficients A k and B k are given by

K -(kω) 2 M -kωD kωD K -(kω) 2 M A k B k = S k C k (5)
However, it may be difficult to determine the Fourier coefficients C 0 , S k and C k (with 1 ≤ k ≤ m) from B 0 , A k and B k (with 1 ≤ k ≤ m) directly due to the complexity of the non-linearities. Cameron and Griffin [START_REF] Cameron | An alternating frequency time domain method for calculating the steady state response of nonlinear dynamic systems[END_REF] suggested that the truncated Fourier expression of F should be calculated by applying an Alternate/Frequency Time domain method (AFT-method)

[B 0 A 1 B 1 • • • A m B m ] F F T -1 → X (t) ⇒ F X, Ẋ, ω, t F F T → [C 0 S 1 C 1 • • • S m C m ] (6)
Following this, the (2m + 1) * n non-linear equations of motion 4 and 5 can be solved by using a non-linear system of equations solver such as the Broyden method [START_REF] Broyden | A class of methods for solving nonlinear simultaneous equations[END_REF].

The path following continuation

In the field of rotating machinery, the behaviour of systems is often calculated for different operational speeds of interest while all the other parameters are kept constant. In this case, it may be useful to apply predictor and corrector mechanisms in order to estimate the non-linear response of the rotorbearing system when rotation speed increases or decreases. We assume that four previous non-linear responses of the rotor-bearing system are obtained. Based on these points, the estimated point on the solution branch can be predicted at a given arc length by applying the Lagrangian polynomial extrapolation method. Any point on the solution branch is represented by (X i , ω i ) where X i and ω i define the Fourier coefficients and the parameter (i.e. the rotating speed of the rotor-bearing system). Then, the arc length between two consecutive points (X i , ω i ) and (X i+1 , ω i+1 ) can be given by

∆S i+1 = (X i+1 -X i ) T (X i+1 -X i ) + (ω i+1 -ω i ) 2 i = 0, . . . , 2 (7)
with the arc length parameters given by ∆S 0 = 0, S 1 = ∆S 1 , S 2 = S 1 + ∆S 2 , S 3 = S 2 + ∆S 3 and S 4 = S 3 + ∆S.

The estimation of the following point at distance ∆S can be predicted by using the Lagrangian extrapolation scheme

X 4 ω 4 = 3 i=1 3 j=0,i =j S 3 -S j S i -S j X i ω i i = 0, . . . , 3 (8)

Condensation procedure

If a non-linear system consists of an n-degree of freedom system with non-linear forces associated with q of these degrees of freedom, this system may be considered as a linear structure with p = nq degrees of freedom and having several additional non-linear elements. Therefore it may be of great interest to keep only the q non-linear degrees of freedom [START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF][START_REF] Sinou | Detection of cracks in rotor based on the 2X and 3X super-harmonic frequency components and the crack-unbalance interactions[END_REF]. Equation 1 can be re-ordered by considering the linear transformation X = PY = P [Y p Y q ] T where Y p and Y q contain the p linear degrees of freedom and the q non-linear degrees of freedom, respectively.

p Y q = Fp Fq (9)
Here, we present the condensation procedure used to obtain the Fourier coefficients associated with the non-linear and linear elements of the complete system.

2.

3.1 Determination of the constant terms B q 0 associated with the non-linear elements By substituting equations 2 and 3 in equations 9, and by only equating coefficients for the constant terms corresponding to the first linear equation of the system 4, we obtain Kpp Kpq Kqp Kqq

B p 0 B q 0 = C p 0 C q 0 ( 10
)
where B p 0 and B q 0 are the vectors with the p linear and q non-linear degrees of freedom of the system, respectively. By eliminating B p 0 from the latter equation, B q 0 is given by

B q 0 = Kqq -Kqp K-1 pp Kpq -1 C q 0 -Kqp K-1 pp C p 0 (11)

Determination of the sine Fourier coefficients

A q k associated with the non-linear elements Then, by equating coefficients for the cosine terms for the k th harmonic of equations 5 we obtain

kω DA k + K -(kω) 2 M B k = C k (12)
Similarly, by equating coefficients for the sine terms for the k th harmonic of equations 5 we obtain

K -(kω) 2 M A k -kω DB k = S k (13)
Then, the determination of the Fourier coefficients A k can be obtained by eliminating the Fourier coefficients B k from equation 12 and 13. By premultiplying equations 12 by K -(kω) 2 M -1 and substituting in equations 13, we obtain

K -(kω) 2 M + kω D K -(kω) 2 M -1 kω D A k = S k + kω D K -(kω) 2 M -1 C k (14)
By introducing

T k = K -(kω) 2 M + kω D K -(kω) 2 M -1 kω D (15)
W k = S k + kω D K -(kω) 2 M -1 C k (16)
equation 14 may be rewritten in partitioned form (as done previously for the constant terms)

T k,pp T k,pq T k,qp T k,qq A p k A q k = W p k W q k (17)
Finally, the Fourier coefficients A q k can be determined by eliminating the Fourier coefficients A p k from the previous equations. After calculation, we obtain

A q k = T k,qq -T k,qp T -1 k,pp T k,pq -1 W q k -T k,qp T -1 k,pp W p k ( 18
)
2.3.3 Determination of the cosine Fourier coefficients B q k associated with the non-linear elements

The same procedure can now be applied to the determination of the Fourier coefficients B q k . By eliminating the Fourier coefficients A k from equation 12 and 13 by premultiplying equations 13 by K -(kω) 2 M -1 and substituting in equations 13, we obtain

K -(kω) 2 M + kω D K -(kω) 2 M -1 kω D B k = C k -kω D K -(kω) 2 M -1 S k (19)
Using the relation 15 and introducing

U k = C k -kω D K -(kω) 2 M -1 S k (20)
equation 19 can be rewritten in partitioned form

T k,pp T k,pq T k,qp T k,qq B p k B q k = U p k U q k (21)
Finally, the Fourier coefficients B q k can be determined by eliminating the Fourier coefficients B p k from equations 21. We obtain

B q k = T k,qq -T k,qp T -1 k,pp T k,pq -1 U q k -T k,qp T -1 k,pp U p k (22) 2.3
B p 0 = K-1 pp C p 0 -Kpq B q 0 ( 23
)
A p k = T k,pp -1 (W p k -T k,pq A q k ) (24)
B p k = T k,pp -1 (U p k -T k,pq B q k ) (25)
Moreover, it can be seen that the same condensation process should be applied to the linear elements if the number of linear components is low compared to the number of non-linear components. After calculation, we obtain the following relations

B p 0 = Kpp -Kpq K-1 qq Kqp -1 C p 0 -Kpq K-1 qq C q 0 (26) A p k = T k,pp -T k,pq T -1 k,qq T k,qp -1 W p k -T k,pq T -1 k,qq W q k (27) B p k = T k,pp -T k,pq T -1 k,qq T k,qp -1 U p k -T k,pq T -1 k,qq U q k (28)
For the latter case, the (2m + 1) * q remaining unknowns Fourier coefficients B q 0 , A q k and B q k associated with the non-linear components can be deduced from equations 10, 17 and 21.

Complete process for the determination of the Fourier coefficients and the non-linear dynamical response of the system

In conclusion, the Fourier coefficients B q 0 , A q k and B q k associated with the non-linear components of the system are first determined by considering the (2m + 1) * q equations 11, 18 and 22. Secondly, determination of the Fourier B p 0 , A p k and B p k associated with the linear components of the system are obtained by considering (2m + 1) * p relations 23, 24 and 25. Then, the calculations of the Fourier coefficients C 0 , S k and C k (with

1 ≤ k ≤ m) from B 0 = [B p 0 B q 0 ] T , A k = [A p k A q k ] T and B k = [B p k B q k ]
T are obtained by using the alternate frequency/time approach (see equations 6).

Therefore the previous equations 11, 18 and 22 can be solved by a solver such as the Broyden method [START_REF] Broyden | A class of methods for solving nonlinear simultaneous equations[END_REF].

To aid the reader's comprehension, it can be seen that for a general non-linear system, certain linear degrees-of-freedom of the vector Y p can be transferred and added to the vector Y q of the non-linear degree-of-freedom without loosing the general process presented previously. This operation can be very interesting if keeping the physical linear degree-of-freedom is necessary for the study (i.e. the non-linear behavior of the chosen physical linear degree-of-freedom can be obtained directly by using relations 18 and 22, without calculating expressions 24 and 25). However, in this case, the size of the vector Y q has been increased, thereby increasing calculation time and storage requirements.

Application to the flexible rotor-bearing system

In this section, the harmonic balance method with the condensation procedure on the non-linear coordinates will be applied to the non-linear rotor-bearing system shown in Figure 1.

Description of the flexible rotor system

Firstly, the complete modeling of the rotor-bearing system is presented. The classical equation of motion of the shaft and disc elements will be briefly described. Then, the description and global expressions of the rolling element bearings will be given in detail by considering the kinematics of the rolling elements, the internal clearance and the Hertz contact nonlinearity.

Shaft elements

The shaft is modeled by 13 Timoshenko beam elements with circular cross sections. Each Timoshenko beam finite element has four degrees of freedoms at each node:

M b T + M b R Ẍb + C b + ωG b Ẋb + K b X b = 0 (29)
where M b T and M b R are the translational and rotary mass matrices of the shaft element. C b , G b and K b are the external damping, gyroscopic, and stiffness matrices respectively. ω is the rotational speed. External damping is taken as classical for the sake of simplicity by considering Rayleigh's expression

C b = α M b T + M b R + βK b
where α and β are constant factors of proportionality and internal rotor damping has been neglected.

Rigid disc

The disk is modeled as a rigid disk and can be written as

M d T + M d R Ẍd + ωG d Ẋd = F d (30)
where M d T , M d R and G d are the translational mass, rotary mass and gyroscopic matrices respectively. F d defines the unbalance or any other excitation on the disc.

Flexible supports and coupling

The flexible supports (see Section A and B in Figure 1) and coupling (see Section C in Figure 1) situated at the two locations of the two rolling-element bearings and the right end of the shaft are modeled as two-node linear elastic spring elements. They are formulated from the following equations of motion:

M s Ẍs + K s X s = F s (31)
where M s , K s and F s are the elementary mass matrix, stiffness matrix, and external load vector.

The rolling bearing

A schematic diagram of the rolling element bearing with its fixed frame of reference and the spin direction of the rotor is shown in Figure 2. The rolling bearing is modeled as a two degrees of freedom bearing with radial clearance and Hertz contact between races and rolling elements. In this study, the outer race of the ball bearing is assumed to be fixed to the flexible support and the inner race is assumed to be fixed to the shaft. The precessional angular position of the center of the ball is given by (Harsha et al. [START_REF] Harsha | Non-linear dynamic behaviors of rolling element bearings due to surface waviness[END_REF] and Tiwari et al. [START_REF] Tiwari | Effect of a ball bearing on the dynamics of a balanced horizontal rotor[END_REF])

ω cage = ω R i R i + R o ( 32
)
where R i and R o are the inner and outer race radii. ω defines the rotational speed of the rotor. The varying compliance frequency is

ω V C = ω cage N b = ωN b R i R i + R o (33)
where N b defines the number of rolling elements for the rolling bearing. By considering that the angular spacing of the rolling elements is constant, the angular space between two balls is equal to ∆θ = 2π N b and each ball is located by its angular position changing with time as the shaft rotates

θ k = 2π N b (k -1) + ω cage t , k = 1, . . . , N b (34)
Then the relative radial distance ∆ k between the inner and the outer races at the k th ball position can be expressed by

∆ k = (x o -x i ) cos (θ k ) + (y i -y o ) sin (θ k ) , k = 1, . . . , N b (35)
where x i , x o , y i and y o define the horizontal and vertical displacements of the inner and outer races at the k th ball.

Considering the position of each rolling element and the local Hertzian contact, the restoring force can be estimated

F k = K H (∆ k -δ) n , ∆ k ≥ δ 0 , ∆ k < δ (36)
where δ is the radial clearance between races and rolling element. K H defines the effective stiffness which is the combined stiffness of a ball in relation to the inner race and outer race contacts

K H = 1 1 K 3/2 i + 1 K 3/2 o (37)
The material deformations are elastic and the dimensions of the contact area are small compared to the curvature radii of the contacting bodies. As indicated in equation 36, restoring force F k appears if the k th ball (at the position θ k ) is loaded (i.e. δ k ≥ δ). If δ k < δ, the k th ball is not loaded and no restoring force is generated. It is noted that loading is assumed to be in normal direction at the contacting surfaces and that the deformations at the contact surface are small compared to the dimensions of the contact area. Finally, the global bearing reaction can be obtained by summing the individual restoring force from each of the rolling-element bearings. Then the total restoring force components F X and F Y in X and Y directions can be expressed as

F X = - N b k=1 F k cos (θ k ) (38)
F Y = - N b k=1 F k sin (θ k ) (39)

The complete flexible rotor system

After assembling, the general dynamics equations of the complete rotor system can be written in the following

M Ẍ + (C + ωG) Ẋ + KX = F L + F N L (40)
where Ẍ, Ẋ and X are the acceleration, velocity and displacement vectors. M is the mass matrix, C is the external damping matrix associated with the non-rotating parts, G is the gyroscopic matrix and K is the stiffness matrix. It should be noted that the internal damping has been neglected in this analysis. However, it may be observed that the non-linear method (i.e. the Harmonic Balance Method) works even if the internal damping of the rotor system is taken into account. F L contains the weight forces and the unbalance forces. F N L corresponds to the complete non-linear forces due to the rolling-element bearings. ω defines the rotating frequency of the system. All the parameters of the rolling-bearing element and the rotor system are given in Tables 1 and2.

Numerical studies

As explained previously, the rolling-bearing system investigated here is composed of a horizontal flexible shaft of 1700mm length with a circular cross section of 40mm. The shaft has one disk with a diameter of 400mm and a thickness of 20mm attached 196mm from the right-end of the shaft. The flexible coupling is added at the right-end of the shaft to connect the motor. The rotor shaft is supported by two bearing supports; the two supports are placed 342mm and 1478mm respectively from the end of the shaft with the flexible coupling. The second bearing support illustrated in Figure 1 is composed of a rolling bearing with radial clearance and Hertz contact between the races and rolling elements, as previously explained in Section 3.1.4. All the parametric values for the rolling bearing element and the complete rotor are given in Tables 1 and2.

Here, the non-linear unbalance responses due to the non-linear and linear components will be investigated, after which the non-linear restoring forces and the associated contact evolution at the rolling bearing element will be calculated.

Non-linear unbalance responses

The non-linear unbalance responses and contacts of the non-linear components (i.e. the rotor and stator coordinates at the non-linear rolling bearing elements) are obtained by solving the non-linear equations 11, 18 and 22, the non-linear unbalance responses of the linear components being obtained by using the relations 23, 24 and 25.

The number of harmonics chosen for the approximated solution is selected on the basis of the number of significant harmonics expected in the non-linear dynamical response. In this study, the number of significant harmonic components varies with the rotating speed of the rotor due to the relative contribution of the unbalance forces, gravity and restoring contact forces on each ball bearing. However, the number of harmonic components retained in the solution for the following numerical analysis is equal to twelve due to the fact that the harmonic components become less significant for upper orders. Firstly, it is shown that the complete non-linear dynamic responses (noted "'CP"' in Figures) at the non-linear and linear rotor-bearing elements are very complex, with significant contributions not only of the first and second orders (i.e. the blue and red lines) but also the sixth orders (i.e. the grey lines) when the rotating speed is in the range [0-500] rpm. However, it can be seen that in some cases the first and second order approximations of the truncated Fourier series should be enough to obtain a good approximation of the non-linear dynamics of the rolling-bearing rotor system (see, for example, the speed range ω = [1000; 4000]rpm where the amplitudes of the rotor are maximal). Moreover, it is clearly shown that the contribution of the n th upper orders can be significant at low amplitudes and when the rotor is passing through the 1 n critical speeds. To facilitate understanding, the forward and backward critical speeds for the first and second modes are situated around 1740rpm and 2050rpm, and 2900rpm and 3520rpm, respectively. For example, the 7×, 6×, 5×, 4×, 3× and 2× harmonic components are predominant at 260rpm, 300rpm, 360rpm, 450rpm, 610rpm and 880rpm (290rpm, 340rpm, 400rpm, 500rpm, 670rpm and 1020rpm, respectively) due to the fact that the rotation speed passes through the subharmonic components of the first backward critical speed (and the first forward critical speed, respectively). The same phenomenon can be observed for the second backward and forward critical speeds: when the rotor passes through the 7×, 6×, 5×, 4×, 3× and 2× harmonic components of the backward or forward critical speeds of the second mode (at 430rpm, 500rpm, 600rpm, 750rpm, 970rpm and 1480rpm for the backward mode, respectively; at 470rpm, 540rpm, 660rpm, 820rpm, 1100rpm and 1720rpm for the forward mode, respectively), the amplitudes of the n th associated upper orders increase.

The non-linear responses of the rotor at the non-linear rotor-bearing element and at the left end of the shaft are illustrated in Figures 5 for a radial clearance of δ = 10µm and an unbalance mass of m u = 4g. It can be seen that the non-linear responses of the rotor system are complex with jump phenomenon for rotor speeds between ω = 2100rpm and ω = 2300rpm. This jump is not only observed for the first order, but also for the second, third, fourth and fifth orders, as illustrated in Figures 5. This phenomenon indicates the "'hardening effect"' at the first critical speed due to the contribution of the non-linear terms of the rolling bearing. Therefore, for an excitation level about twice as high as that of the initial case, the peak of the first forward critical speed appears for a rotating speed higher than about 7% with respect to the initial case. The hardening-type nonlinearity on the maximum amplitude of the harmonic components can also be observed. It is then possible to see that two stable solutions coexist when the rotation speed is around [2140 : 2210]rpm. As previously seen for the unbalance case m u = 2g (Figures 3 and4), the global non-linear responses are very complex with the significant presence of the first, second and sixth orders (see the horizontal amplitudes at the left end of the shaft). Even if the non-linear behavior between the two cases (m u = 2g and m u = 4g) are different (with the appearance of jump phenomenon for example), the global contributions of each order for the two cases are very similar: the maximum amplitude is observed when the rotation speed exceeds the first critical speed (ω = 2000rpm), the most important contributions are given by the first and second orders when the speed reaches about [1000-4000]rpm, the sixth order makes a significant contribution at low speed (around [0-500]rpm), and the n th orders appear when the rotor speed reaches around 1 n of critical speeds. Finally, it can be seen that the nonlinear responses at the left end of the shaft (see the line "'CP"' in Figures 4 and5) are governed by the amplitudes of the sixth order when the rotation speed is between [0-1300]rpm: the contribution of the sixth order is most significant when the rotation speed is between [0-500]rpm and one of the highest amplitudes if the rotor speed reaches around [500-1300]rpm.

The orbits at the disc position, left end and middle of the shaft are shown in Figures 6 for various rotation speeds between [300-1500]rpm when the components of the n th orders have significant contributions. Thus it is clear that the non-linear responses can be very complex with "'multiple inside or outside loops"'. Figures 7 show the complete non-linear amplitudes for the rotor shaft and the two stators of the rolling bearing elements at specific rotation speeds. For each case, the orbits associated with the non-linear degrees of freedom (and thus calculated by using the condensation process and the relation equations 11, 18 and 22) are indicated by red lines. The orbits corresponding to the linear degrees of freedom obtained by using the expressions 23, 24 and 25 are indicated by the black lines. The blue line illustrates the deflection of the rotor shaft. To facilitate comprehension, the top subfigures illustrate the rotor amplitudes of the linear and non-linear degrees of freedom, and the stator amplitudes at the rolling bearing element are presented by the lower subfigures (for the linear and non-linear degrees of freedom). It should be noted that the orbits and the deflection of the shaft are presented without the static deflection of the rotor in order to clearly show the complexity of the non-linear response that generally occurs at low amplitudes. (c) and(d) show the shaft's deflection at the backward and forward first critical speeds (for ω = 1737rpm and ω = 2052rpm). The shaft's deflections for the backward and forward second critical speeds (for ω = 2932rpm and ω = 3526rpm) are given in Figures 7 (e) and (f), respectively. It may be observed that the rotor's orbits are very simple and similar to a simple loop when the rotation speed is passing through the second critical speeds. For the first critical speeds, simple loops are shown at the right-end of the shaft and at the disc position. However, inner loops appear for the orbits at the left position, at the middle of the shaft and at the first bearing support. These results are in perfect agreement with the previous evolutions of the Fourier components: the second order of the Fourier series is sufficient to describe the non-linear behavior and the associated orbits of the nonlinear rotor system at the first critical speeds, and the first order gives an adequate approximation of the rotor orbits at the second critical backward and forward critical speeds. Finally, Figures 7 (a) and (b) illustrate that the orbits can be more complex (as previously seen in Figures 6): the contributions of the n th order are shown for all the shaft orbits. All these orbits shown in Figures 6 and7 clearly indicate the complexity and variability of the non-linear responses of the rotor-bearing system.

Non-linear contact forces and orbits at the rolling-bearing element

In this part of the paper, the contacts between races and balls and the associated restoring forces in fixed coordinates (OXY in Figure 1) are investigated in order to understand the relative contribution of the unbalance and gravitational forces for the rotor, the shaft whirling motion, and the associated nonlinear behavior at the rolling-bearing element. Due to the condensation process previously presented and by using the relations 11, 18 and 22, the non-linear behavior of the rotor system is only estimated at the non-linear degrees of freedom in this part of the study.

In order to carry out this analysis, the evolution of the contact on each ball is first evaluated and discretized over one period of revolution of the rolling-bearing cage

ω cage = ω(R i /(R i + R o ))
where R i and R o are the outer and inner race radii of the rolling bearing element and ω defines the rotation speed of the rotor. Secondly, the global restoring force is calculated in the fixed coordinates not only over one period of revolution of the rolling-bearing cage but also over fifty periods, in order to show the envelops of the non-linear contact forces of the rolling bearing element. Figures 8 and9 show the evolution of the contact restoring forces for each ball of the bearing for a radial clearance of δ = 10µm, and an unbalance mass of m u = 2g and m u = 4g, respectively. The white lines define the limits between the contact and no-contact zones for each ball. Figures 10 and11 give the global restoring forces for both cases in the fixed coordinates. The blue lines represent the evolution of the non-linear contact over one period of revolution of the rollingbearing cage. The green surface corresponds to the envelop of the non-linear contact over fifty periods of revolution. The red line indicates that the rotor and stator are in contact with a non-linear force greater than zero. Figures 9 and11 give contact evolutions for all the speed ranges of interest (i.e. ω = [0 -4000]rpm and for m u = 4g) whereas Figures 8 and10 focus on the first critical speeds for m u = 2g.

For the first unbalance mass m u = 2g, the rotor moves at the bottom of the bearing for different rotation speeds (ω = 1895rpm, ω = 1990rpm, ω = 2215rpm and ω = 2244rpm), as indicated in Figures 10(a,b,e,f) with the red lines that define the non-linear contact of the rolling-bearing element in fixed coordinates. Regarding the associated non-linear contacts on each ball of the bearing (see Figures 8(a,b,e,f)), five or six balls are in contact with the outer ring each time and the number of times of contact and non-contact per ball is similar for each ball. For ω = 1895rpm and ω = 2244rpm, the evolution of contact for each ball is due to the revolution of the rolling-bearing cage: the area within which the rolling elements are still in contact with the raceway, generally referred to as the loaded zone, is situated at the bottom of the ball bearings (in fixed coordinates). The value of the non-linear contact force is still high, as shown in Figure 10(a) and (f) (see the green surfaces that indicate the contact over fifty periods of revolution of the rolling-bearing cage). When the rotor exceeds the first forward critical speed (i.e. ω = 1990 -2130rpm), the contact value at the bottom of the bearing can be equal to zero (see the green surface in Figures 10(b,c,d,e)) and the red lines correspond to the full circle. Moreover, the evolutions of the non-linear contacts on each ball become more complex, as indicated in Figures 8(b,c,d,e). The contact and non-contact times per ball can differ considerably from one ball to another and the associated restoring contact forces can increase or decrease during one period of revolution of the rolling-bearing cage. Finally, Figures 10(c) and 8(c) illustrate the case of a complete whirling contact on the rolling element (i.e. the radial clearance in the fixed coordinates is consumed while the red line defines the full "'circle"' in Figure 10(c)). It is noted that the value of the non-linear restoring forces is still higher at the bottom of the rolling element. Due to cage rotation, the maximum value of the non-linear contact changes from one ball to another, as indicated in Figure 8(c). Finally, it can be seen that the period of contact/no-contact between the rotor and each ball in this case (i.e. ω = 2030rpm and m u = 2g, Figure 8(c)) decreases in comparison to the previous cases (i.e. ω = 1895rpm or ω = 2244rpm for m u = 2g, Figures 8(a,f)).

Then, for the second unbalance mass m u = 4g, the same non-linear evolutions are observed around the first critical speed: each ball participates with almost the same contact time over the period of rotation of the rolling-bearing cage and the non-linear contact forces can increase or decrease during one period. When the rotor whirls on the bottom of the bearing, the evolution of the contact forces corresponds to one period of revolution of the rolling-bearing cage. At the first critical speed (i.e. ω = 2100 -2200rpm), the rotor whirls completely in the rolling-bearing element and the clearance is consumed in the fixed coordinates during its rotation (i.e. the red lines define the "'full circle"', as shown in Figure 11). Moreover, due to the increase of the unbalance force, the contact on all the rotor-bearing elements (i.e. the "'red circle"') is present for a greater speed range and the values of the non-linear contact forces at the top and bottom of the rolling element are more higher (see Figure 11 for ω = 2195rpm). In this case, the unbalance effect is predominant for the non-linear behavior of the contact/no-contact evolutions. As seen in Figure 9 for ω = 2204rpm three or four contacts occur per period of revolution of the rolling-bearing cage for each ball. Each ball participates with almost the same contact time over the period of rotation. As previously shown, the period of contact/no-contact between the rotor and each ball is similar to the period of the rotation of the rotor ω. However, if the rotation speed is around the first critical speed, the contact time and the restoring contact force are very different and complex per period of revolution for each ball (see Figures 9 for ω = 1919rpm ω = 2253rpm) even if the rotor and stator orbits are simple and close to a periodic harmonic form (see results of the previous section). Finally, the non-linear contact forces are given for the speed range of interest ω = 0 -4000]rpm in Figures 11. It can be seen that the non-linear contact over one period of revolution of the cage (i.e. the blue line) or the envelopes of the non-linear contact (i.e. the green surface) are totally different for the different rotation speeds and can be complex in some cases. For example, the non-linear contact can be dissymmetrical on the rotor-bearing element (for example, at the following rotation speeds: ω = [1696; 1787; 2703; 2890]rpm). These dissymmetric non-linear restoring forces at the rolling element are observable at the first and second backward critical speeds where the rotor orbits are elliptical. Then, at low amplitudes, the value of the non-linear force can remain constant due to the fact that the rotor moves on the bottom of the rolling element: see for example Figure 11 for the following rotation speeds ω = [296; 1595; 4000]rpm where the blue lines (i.e. evolution of the non-linear contact forces over one period of revolution of the rolling-bearing cage) are similar to the green surface (i.e. evolution of the non-linear contact forces over fifty periods of revolution of the rolling-bearing cage). As explained previously, the contact and non-contact times per ball and the value of the non-linear forces are similar from one ball to another, in this case due to the rotation of the rolling cage.

Finally, it can be seen that the non-linear contact force at the rolling element can be different for a given rotation speed if in run-up or run-down configurations due to the jump phenomenon. This fact is illustrated in Figure 12 for m u = 4g and δ = 10µm.

Influence of the mass unbalance and the radial clearance

In this section, the influence of the mass unbalance m u and the radial clearance δ on the non-linear responses, the Fourier components and the non-linear contact forces and orbits at the rolling-bearing element are investigated. Figures 13 and14 show the non-linear amplitudes and the Fourier components relating to the unbalance and the radial clearance, respectively. Figures 15 illustrate the associated non-linear contact forces at the rolling-bearing element. Firstly, the unbalance mass is increased from 1g to 4g. In the simulations, a reduction of the unbalance mass appears to be an effective way to reduce the vibration level of the non-linear response and the subharmonic components if the rotation speed is higher than 1000rpm, as illustrated in Figures 13. Initially a softening-type nonlinearity, due to the non-linear contribution of the rolling bearing, is observed for low excitation levels (i.e. m u = [1; 2]g). This "'softening effect"' is observed on the first and second orders around the first forward critical speed. For an excitation level of about 0.5 times the initial case (from m u = 2g to m u = 1g), the peak of the response appears for a frequency lower than about 2% in comparison to the initial case (i.e m u = 2g), even if the amplitude increases. For the two highest excitations (i.e. m u = 3g and m u = 4g), a hardening-type nonlinearity is obtained with jumps. Jumps can be observed not only for the first order and the complete non-linear response (i.e. Composite Power), but also for the second harmonic component, as shown in Figure 13. For these two highest excitations, increases of about 2% and 8% respectively of the frequency of the maximum amplitude are observed. Moreover, an increase of the rotation speed interval where two stable solutions coexist is obtained due to the hardening-type nonlinearity of the rolling bearing. For low rotational speeds (i.e. ω = [0 -1000]rpm), it is noted that the non-linear responses and the most important sub-harmonic component in the amplitudes (i.e. order 6) are not changed. Secondly, Figures 14 give the non-linear amplitudes and the Fourier components for the radial clearance. In this example, the radial clearance is decreased from 30 µm to 0µm. It can be seen that that the radial clearance influences the critical speeds of the rotor. By increasing the radial clearance, a small softening-type nonlinearity can be observed. For the largest radial clearance (30 µm), the peak of the response appears for a frequency lower by about 4% in comparison to the initial case (10 µm). However, the non-linear response and the subharmonic components appear to be very similar for the different radial clearances under study. Finally, Figures 15 illustrate the non-linear contact forces for the different unbalances and radial clearances for the rotation speed around 2000rpm. It should be noted that the unbalance mass influences the non-linear restoring forces at the rollingbearing element: increasing the mass unbalance allows the rotor to whirl in the cage with an increase of the non-linear contact forces. A reduction or increase of radial clearance has very little effect on the non-linear contact forces.

b) ω = 1595rpm (c) ω = 1696rpm (d) ω = 1797rpm (e) ω = 1916rpm (f) ω = 1996rpm (g) ω = 2195rpm (h) ω = 2253rpm (i) ω = 2703rpm (j) ω = 2890rpm (k) ω = 3004rpm (l) ω = 4000rpm (

Conclusion

This Harmonic Balance Method with a condensation procedure was developed in order to analyze the non-linear behavior of a rotor system with a rolling bearing elements. In this study, twelve harmonic components were chosen in order to approximate the non-linear response. The non-linear elements are both due to the radial clearance and the Herztian contact between races and rolling elements. The results demonstrated that the non-linear unbalance response and the associated orbits can be very complex due to the non-linear elements and the relative contribution of unbalance forces and radial clearance: the n th Fourier components can make a significant contribution when the rotor passes through the 1 n sub-critical resonances and critical speeds. In particular, the 2× and 3× and 6× superharmonic frequency components and 1 2 , 1 3 and 1 6 sub-critical resonances are significant. It appears that the vibration amplitudes in the sub-critical and critical resonances depend on the radial clearance and unbalance mass. If the unbalance mass and radial clearance remain constant during operation of rotor, the 1×, 2× and 3× amplitude may change during run-up or rundown due to the occurrence of the jump phenomenon. The changes in the non-linear contacts at each rolling bearing element were examined in details, with particular attention given to when the rotor exceeds critical speeds. It was shown that the contact evolution for each ball-bearing can be very simple even if the non-linear behavior of the rotor and the associated orbits are complex. Finally, it was shown that mass unbalance and radial clearance affect the non-linear contact forces and the whirling motion of the rotor at the rolling-bearing element. In particular, it was shown that increasing small excitation levels leads to an initial softening-type nonlinearity, which turns into hardening-type for higher excitations.

Nomenclature

X displacement vector Ẋ velocity vector Ẍ acceleration vector M mass matrix of the rolling-bearing rotor system K stiffness matrix of the rolling-bearing rotor system G gyroscopic matrix of the rolling-bearing rotor system C damping matrix of the rolling-bearing rotor system F L vector of the weight and unbalance forces F N L vector of the non-linear forces due to the rolling-element bearing m order of the Fourier series A k Fourier coefficients of the sinus function for the k th order B k Fourier coefficients of the cosinus function for the k th order Y p p linear degrees of freedom of the rolling bearing rotor system Y q q non-linear degrees of freedom of the rolling bearing rotor system B q k Fourier coefficients of the sinus function for the k th order of the non-linear dof A q k Fourier coefficients of the cosinus function for the k th order of the non-linear dof B p k Fourier coefficients of the sinus function for the k th order of the linear dof A p k Fourier coefficients of the cosinus function for the k th order of the linear dof

Figure 1 :Figure 2 :

 12 Figure 1: Rotor system

  The complete vertical and horizontal unbalance responses at the non-linear rotor-bearing element and the contribution of the first twelve orders for the non-linear components are shown in Figures 3 (for a radial clearance of δ = 10µm and an unbalance mass of m u = 2g). The non-linear unbalance responses of the linear components of the rotor system at the left end of the shaft are given in Figure4. These non-linear amplitudes are obtained by considering the previous relations 23,24 and 25.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Non-linear unbalance response at the non-linear rotor-bearing element for δ = 10µm and m u = 2g (a) horizontal (b) vertical

Figure 6 :

 6 Figure 6: Evolution of complex orbits for δ = 10µm and m u = 2g (b,c) disc position (d,e,g,h) left end (a,f) middle of the shaft

Figure 7 :

 7 Figure 7: Evolution of the orbits for δ = 10µm and m u = 2g

  0

Figure 8 :

 8 Figure 8: Evolution of the contact on the balls for δ = 10µm and m u = 2g

Figure 9 :Figure 10 :

 910 Figure 9: Evolution of the contact on the balls for δ = 10µm and m u = 4g

Figure 11 :Figure 12 :

 1112 Figure 11: Evolution of the contact at the rolling-bearing for δ = 10µm and m u = 4g

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Influence of the unbalance mass on the rotor unbalance responses (black=composite power, blue=order 1, red=order 2, magenta=order 3, green=order 4, cyan=order 5, grey=order 6, -m u = 1g, --m u = 2g, -. m u = 3g, • • • m u = 4g)

.4 Determination of the Fourier coefficients B p 0 , A p k and B p k associated with the linear ele- ments

	Now, the (2m + 1) * p remaining unknown Fourier coefficients B p 0 , A p k and B p k can be estimated from
	equations 10, 17 and 21. After calculations, we obtain

Table 1 :

 1 Geometric properties of the rolling-bearing element

	10

Table 2 :

 2 Parametric values of the flexible rotor system


J-J Sinou
email: jean-jacques.sinou@ec-lyon.fr

An experimental investigation of condition monitoring for notched rotors through transient signals and wavelet transform

Experimental studies considering the dynamic transient signals of a notched rotor (i.e. a rotor with a transverse open crack) are investigated in order to examine the possibility of detecting the presence of open cracks in rotating machinery for low or high rotor accelerations. Firstly, the orbital patterns, spectrum cascade plots of the Power Spectral Density, evolutions of the 2× amplitudes of rotor with and without an open transverse crack are compared and discussed in the case of low accelerations of the rotor system. Secondly, the non-stationary vibration signals of the rotor with an open transverse crack at high accelerations are investigated. The timefrequency features and a tool based on the continuous wavelet transform (CWT) to detect open cracks in a rotor system are discussed. It is demonstrated that both the continuous wavelet transform (CWT) and changes on the 2× harmonic components are robust indicators for crack detection.

Introduction

Structural damage detection and condition monitoring have been receiving increased attention in recent years due to their practical importance in many industrial applications. Extensive reviews have been compiled by Wauer [START_REF] Wauer | Dynamics of cracked rotors: Literature survey[END_REF] and Dimarogonas [START_REF] Dimarogonas | Vibration of cracked structures: a state of the art review[END_REF]. Generally, two types of approaches are suggested by different researchers. One monitoring changes in the resonant frequencies and in the mode shapes and other based on the change in non-linear dynamic responses during rotation of the notched shaft. The first approach is an off-line process for rotating machinery. The modal tests and detection based on changes in modal parameters have to be carried out during static condition of rotor. Reviews of research works dealing with the problem of crack detection based on changes in modal parameters can be found in [START_REF] Montalvao | A review of vibration-based structural health monitoring with special emphasis on composite materials[END_REF][START_REF] Doebling | Damage Identification and Health Monitoring of Structural and Mechanical Systems from Changes in Their Vibration Characteristics: A Literature Review[END_REF][START_REF] Sohn | A Review of Structural Health Monitoring Literature : 1996-2001[END_REF][START_REF] Salawu | Detection of structural damage through changes in frequency: a review[END_REF]. However, several researchers [START_REF] Imam | Development of an on-line rotor crack detection and monitoring system[END_REF][START_REF] Sinou | Damage assessment based on the frequencies' ratio surfaces intersection method for the identification of the crack depth, location and orientation[END_REF] have demonstrated that changes in modal properties are not always reliable and this type of damage identification technique can become very difficult. The second approach focuses on the non-linear dynamic vibration of rotor during its rotation. Many researchers [START_REF] Davies | The vibrational behaviour of a multi-shaft, multibearing system in the presence of a propagating transverse crack[END_REF][10][START_REF] Gasch | A survey of the dynamic behavior of a simple rotating shaft with a transverse crack[END_REF][START_REF] Pu | Quasi-periodic vibration of cracked rotor on flexible bearings[END_REF][START_REF] Zhu | The dynamics of a cracked rotor with an active magnetic bearing[END_REF][START_REF] Sinou | Influence of cracks in rotating shafts[END_REF] indicated that the diagnosis of crack presence can be done by considering the appearances or/and changes in the steady-state vibration amplitudes of 1×, 2× and 3× harmonic components of steady-state response of the rotor system. They also proposed to exploit orbital patterns of the notched rotor at one-half of the first critical speed [START_REF] Adewusi | Experimental study on the vibration of an overhung rotor with a propagating transverse crack[END_REF][START_REF] Sinou | Effects of a crack on the stability of a non-linear rotor system[END_REF]. On the other hand, Bachschmid et al. [START_REF] Bachschmid | Identification of multiple faults in rotor systems[END_REF][START_REF] Pennacchi | A model-based identification method of transverse cracks in rotating shafts suitable for industrial machines[END_REF] proposed a model-based diagnostic approach with a least-squares identification for the detection of cracks in rotors. Darpe [START_REF] Darpe | Experimental investigations of the response of a cracked rotor to periodic axial excitation[END_REF] proposed to use coupling of bending and longitudinal vibrations with periodic axial excitation to detect the crack presence. Han [START_REF] Han | Vibration analysis of periodically time-varying rotor system with transverse crack[END_REF] illustrated that directional frequency response functions (d-FRFs) is a robust tool for investigating the presence or the propagation of the crack. Then, non destructive detection of cracks in rotors based on transient non-linear signals have been developed [START_REF] Prabhakar | Transient lateral analysis of a slant-cracked rotor passing through its flexural critical speed[END_REF][START_REF] Sekhar | Detection and monitoring of crack in a coast-down rotor supported on fluid film bearings[END_REF][START_REF] Darpe | A novel way to detect transverse surface crack in a rotating shaft[END_REF]. For example, Sekhar [START_REF] Sekhar | Condition monitoring of cracked rotors through transient response[END_REF] investigated the effects of different factors such as crack depth, unbalance eccentricity and acceleration influencing transient vibrations and concluded that the measurement of transient response when the rotor is passing through the critical speed or one-half of the critical speed can be used for the non destructive detection of cracks in rotors. Finally, some researchers [START_REF] Sekhar | Crack detection and monitoring in a rotor supported on fluid film bearings: Start-up vs run-down[END_REF][START_REF] Sekhar | Crack detection through wavelet transform for a run-up rotor[END_REF][START_REF] Adewusi | Wavelet analysis of vibration signal of an overhang rotor with a propagating transverse crack[END_REF] proposed to introduce time-scale signal processing tools based on wavelet transform due to the limitation of the conventional Fourier analysis that is suitable for steady state vibration signals, but provides a poor representation of signals well localized in time. This method for non-stationary signal analysis has been recently used for various fault diagnosis and damage localization [START_REF] Yanping | Continuous wavelet grey moment approach for vibration analysis of rotating machinery[END_REF][START_REF] Gaul | Identification of the impact location on a plate using wavelets[END_REF][START_REF] Junsheng | Application of an impulse response wavelet to fault diagnosis of rolling bearings[END_REF]. However, few experimental researches could be found on the detection of cracks in rotating machinery with wavelet techniques.

So, the objectives of this work is to present an experimental study on the dynamic response of a rotor with a transverse open crack for low and high rotor's accelerations and to show the ability and feasibility of the application of Continuous Wavelet Transform (CWT) in the diagnostic of cracks in the rotating machinery using the non-stationary vibrations signals during machine run-up.

In this study, the type of crack that has been chosen is a transverse open crack: the crack that always remains open is perpendicular to the shaft axis. As explained by Sabnavis et al. [START_REF] Sabnavis | Cracked shaft detection and diagnostics: A literature review[END_REF], most experimental studies are focused on open cracks (called notches or gaping cracks) due to the fact that this crack type are easy to mimic in a laboratory environment.

Firstly, the conventional fast Fourier transform (FFT)-based spectral analysis method is used to analyze steady-state vibration signals of the rotor with and without crack. Secondly, the non-stationary transient signal of rotor with high accelerations is analyzed by using the continuous wavelet transform. The effectiveness of wavelet transforms for crack detection and monitoring in rotors will be illustrated for various experimental tests.

Experimental Setup

The test equipment used for the experiment include a test rig, motor speed control, eddy current displacement/proximity probes, and data acquisition interface unit that is controlled by the I-DEAS Software. The test rig as Figure 1 illustrates consists of an elastic shaft and a rigid disk. The rotor's shaft is simply supported at both ends (designed by Support 1 and Support 2 in Figure 1). The length of the shaft is 480mm and its diameter is 10mm. The disc is mounted at the position 260mm from the left end of the rotor. The diameter and thickness of the disc-rotor are 75mm and 35mm, respectively. The rotor system is driven by an electric motor incorporated with a shaft through a flexible coupling. The first critical speed of the crack-free rotor is 1980 rpm. The rotor's crack is made at the position 240mm by cutting a piece with width 0.2mm and a non-dimensional depth of µ = 1 (corresponding to the loss of half the shafts area). It must be noted that the cut that has been made with a jeweler's saw introduces an open crack in the rotor. Indeed the lateral deflections being small, the slit may not completely close and would The vibration measurements are conducted by using six proximity probes that are mounted in both the horizontal and vertical directions, as indicated on Figure 1. They consist of two orthogonal sensors for each plane and are placed at 70mm (plane A1), 210mm (plane A2) and 440mm (plane A3) from the left end of the shaft. Planes A1 and A3 are situated near the left and right supports respectively, and the plane A2 is near the crack's position. The probes are connected to the National Instruments PCI-4472 channel dynamic signal acquisition devices, which is in turn connected to the PC. The data acquisition interface unit operation is controlled by the I-DEAS Software to collect and store vibration data. Then the data are transferred from I-DEAS Software to Matlab in order to be processed.

Results and discussions

Steady-state dynamic responses in time domain

Orbital patterns and 2× amplitudes at low acceleration

First of all, the dynamic responses of the healthy and notched shafts is analyzed for low acceleration of rotor (20rpm.s -1 ). Figure 2 shows the vertical and horizontal vibrations of the rotor without and with open crack passing through its first critical speed in time domain.

An increase in amplitudes of vibration is observed with the notched rotor compared to the healthy ro-tor in both the vertical and horizontal amplitudes when the rotor is passing through its first critical speed (for t = 100s). Moreover, there is large increase of vibrations near the rotational speeds at one-half of the first critical speed (for t = 50s).

Evolutions of the orbital patterns of the healthy and notched rotors when the rotor is passing through one half of the first critical speed are shown in Figures 3 and4. For the reader comprehension, the rotating speed is increasing from Figures 3(a) to 3(h) (from 4(a) to (h) respectively). It can be clearly observed that the presence of a transverse open crack induces a loop containing an important loop inside. When the rotating speed of the notched rotor reaches through half of the first critical speed, a rotation of the inside loop indicating a significant change in the phase and amplitude of the 2× vibration appears. For the healthy rotor, the orbit firstly changes from double loops to an outside loop (like a "'eight"') and then a small inside loop is observed. This non-linear vibrations of the healthy rotor and the small contribution of super-harmonic components can be attributed to the bow-related side load on the shaft, the rotor anisotropy and the unbalance. Comparing the orbital patterns of the healthy and notched rotors (3 and 4), it can be concluded that the size of the inside loop and the orbital movements of the rotor through the passage of one-half of the first critical speed can be used as the signature of the presence of a open transverse crack in rotating machinery.

To identify the presence of sub critical response peaks of the dynamic responses of the healthy and notched shafts, spectrum cascade plots of the Power Spectral Density (PSD) which describes how the power of time series is distributed with frequency, are calculated (as previously explained the signal is assumed to be stationary for small intervals of time in this first section of the study). Figures 56(d) and 5-6(b) illustrate the power spectral density when the rotor is passing through the first critical speed (t=100s) and one-half of the first critical speed (t = 50s). It is clearly observed that the 2× super-harmonic frequency component is more prominent than the rotational frequency near the rotational speeds at the 1 2 of the first critical speed for the rotor without and with crack. However, the maximum value of power spectral density is obtained for the notched rotor. So an increase of the power spectral density at one-half of the critical speed can be used as an indicator of the presence of a transverse open crack in rotors. Then, for the both cases, only the 1× revolution vibration amplitude is present when the rotor is passing through its critical speed. Finally, before the 1 2 of the critical speed (for t = 40s) and between 1 2 and 1 of the critical speed (for t = 80s), the motions of the notched rotor are complicated due to existence of the 1× and 2× frequencies of the rotational speed in the vibration signals (see Figures 6(a) and (c)), whereas the motions of the crack-free rotor are principally described by the 1× revolution vibration amplitudes (see Figures 56(a) and (c)). To clarify the contributions of the n× super-harmonic frequency components in run-up stages, the normalized conventional power spectral density cascade plots are drawn in Figures 7 for the rotor without and with an open crack. It is noted that the characteristic sub-critical response peaks and an average value of power spectral density can be estimated at various regular intervals of times due to the low acceleration of the rotor during the transient analysis (i.e. for each interval, the transient vibrations are assumed to be stationary). The normalized spectral density are defined by

P SD normalized = log (P SD) -log (P SD min ) log (P SD max ) -log (P SD min ) (1)
where P SD min and P SD max define the minimum and maximum values of the power spectral density on the speed range. The power spectral density is normalized in order to be able to compare the contributions and evolutions of the n× super-harmonic frequency components for the rotor with and without crack. The values of the normalized power spectral density varies from zero to unity. The physical values have been given in the previous power spectral density plots (see Figures 5 and6).

As indicated in Figures 7, the contribution of 2× harmonic amplitudes for the notched rotor is more significant and greater than that of the healthy rotor when the rotor reaches one-half and one-third of the first critical speed for both the vertical and horizontal directions. Moreover, it can be observed that the 3× harmonic components are present in the horizontal direction when the rotor is passing through nearly the first critical speed. However, the differences of the 3× amplitudes for the healthy and notched rotors are not significant enough in order to be used as a diagnostic tool.

In conclusion, these experimental results reflect the fact that the contributions of the 2× amplitudes and the orbital patterns at one-half of the first critical speeds could be key indicators for the detection of transverse open cracks in a rotating shaft.

Amplitudes for low and high accelerations of the rotor

If the rotor's acceleration is very high, the vibration responses are typically transient along with frequency and amplitude modulation because of high acceleration of the rotor system, as illustrated in Figures 8. In this case, the conventional Fourier analysis can provide a poor representation of signals well localized in time. So, the evolutions of the orbit patterns or the conventional spectrum cascade plots can not be used for analyzing vibration characteristics and for detecting the presence of a transverse crack in rotor clearly. However, the maximum amplitudes for one-half of the first critical speed and the critical speed can be investigated for various accelerations of the rotating machinery. The results are presented in Figures 9. The various accelerations are given in Table 1 for each case. First of all, it clearly appears that increasing the rotor's acceleration decreases the amplitudes at the first critical speed both in the vertical and horizontal amplitudes (see Figure 9(a)).

Then, it is observed that the maximum horizontal and vertical amplitudes at one-half the first critical speed are not greatly affected by the rotor's acceleration (see Figure 9(b)). This reflects the fact that the 2× super-harmonic frequency components for the notched rotor are predominant not only at one-half of Finally, Figures 10 show the evolutions of the maximum amplitudes versus unbalance : it clearly appears that the maximum amplitudes at the first critical speed corresponds mainly to the unbalance effect whereas the maximum amplitude at one-half of the first critical speed correspond to the 2× subcritical resonance and mainly to the crack effect.

So the comparisons of maximum amplitudes at one-half of the first critical speed between the rotor with and without an open crack can be applied to prognosis diagnosis of crack in a rotor system for low or high rotor's accelerations.

Wavelet analysis with Continuous Wavelet Transform

Due to the fact that the characteristics of transient responses cannot be precisely described and analyzed by means of the conventional spectrum cascade plots, the experiments will now be conducted at higher rotor accelerations to explore the effectiveness of the crack detection via the Continuous Wavelet Transform. Firstly, a brief basic theory of the wavelet analysis is given. Then experimental tests for various high accelerations of the rotor system are investigated and the detection of cracks via the continuous wavelet transform are undertaken.

Basic theory

The wavelet analysis transforms a signal into wavelets that are well localised both in frequency and time.

The continuous wavelet transform (CWT) of a function f (t) is a wavelet transform defined by

W (a, b) = +∞ -∞ f (t) ψ * a,b (t) dt (2)
where

ψ a,b (t) = 1 √ a ψ t -b a (3)
are the daughter wavelets (i.e. the dilated and shifted versions of the "'mother"' wavelet ψ that is continuous in both time and frequency). a is the scale parameter, b is the time parameter. The asterisk ψ * a,b

indicates the complex conjugate of ψ a,b .

The following admissibility condition has to be satisfied 0 < C ψ < +∞ where C ψ defines the admissibility constant

C ψ = +∞ -∞ | ψ (ω) | 2 |ω| dω ( 4
)
and ψ is the Fourier transform of ψ ( ψ = +∞ -∞ f (t) e iωt dt). In this study, the Morlet mother wavelet has been chosen due to the fact that it is one of the most commonly used CWT wavelet (quite well localized in both time and frequency space). It is defined as following in the time domain: where m is the wavenumber and η is non-dimensional time parameter. For a time signal f (t) represented by N sampled data points (with uniform time step δt), the Continuous Wavelet Transform of Equation 2 is a convolution of the data sequence f (n ′ ) (with n ′ = 1, . . . , N ) with a scaled and normalized wavelet. It can be represented as follows:

ψ 0 (η) = π -1 4 e imη e -
C (a, n) = N -1 n ′ =0 f n ′ δt a Ψ * 0 (n ′ -n) δt a (6)
where n defines the localized time index and δt is the sampling interval. The wavelet function contains unit energy at every scale due to the normalization of the mother wavelet. Then, the wavelet power is defined as

|C (a, n) | 2 .

Crack detection using wavelet transform

First of all, the wavelet analysis is used to compare the transient responses of the rotor with and without an open crack. Figure 11 shows the normalized wavelet power spectrum

|C(a,n)| 2 σ 2
where σ 2 is the variance [START_REF] Torrence | A practical guide to wavelet analysis[END_REF] for the two cases previously detailed in Section 3.1). Comparing Figure 11(a) and (b), the conclusion is that the continuous wavelet transform can identify the notched rotor effectively. The top figure correspond to the horizontal direction and the bottom to the vertical direction. It clearly appears that the presence of a transverse crack induces the appearance of a second maximum value of the normalized wavelet power spectrum (which corresponds to the sub-harmonic resonance of the notched rotor) when the rotor is passing through one-half the first critical speed (for t ∈ [45 -50]s). These results reflect the previous observations done in Figures 7. Now the experiments are conducted at higher rotor accelerations to explore the effectiveness of the crack detection via the Continuous Wavelet Transform. The measured vibration behaviour of the rotor both in vertical and horizontal directions in time are shown in Figures 12(a),(c) and (e) for Cases 6, 8 and 11 (see Table 1). The associated normalized wavelet power spectrum are given in Figures 12(b), (d) and(f). It clearly appears that there is no significant change in performance of the detection of crack via the continuous wavelet transform with high rotor accelerations : near the rotational speeds at the 1 2 of the first critical speed for the rotor, the normalized wavelet power spectrum increases drastically indicating the presence of the transverse open crack in the rotating machinery. The wavelet power spectrum consists of both unbalance and crack effects : the first peak corresponds mainly to the crack effect and the associated 2× sub-critical resonance and the other, mainly to the unbalance and the interaction with the crack. Finally, Figures 13 give the normalized wavelet power spectrum for two unbalance cases. Even if the unbalance is important (see Figure 10 (a) at the critical speed (b) at 1 2 of the critical speed Figure 9: Maximum amplitudes at one-half of the first critical speed and at the critical speed for various rotor's accelerations (Aiv: vertical amplitude for the i th plane, Aih: horizontal amplitude for the i th )

and the wavelet analysis can be used to identify the presence of a crack in rotating machinery.

In conclusion, the change in the wavelet power spectrum and the transient dynamic behavior of the notched rotor can be used as a diagnosis feature for low and high rotor's accelerations.

Conclusion

In this paper, it is demonstrated that the diagnosis of the presence of a transverse open crack based on the transient vibration signals during machine run-up can be easily done by using various appropriated tools. For steady-state vibration signals (low rotor's accelerations), the traditional spectrum cascade plots of the power spectral density and the evolution of orbits patterns with the classical two-loops orbit can be

Figure 1 :

 1 Figure 1: Test rig of the rotating machinery

Figure 2 :

 2 Figure 2: Dynamic responses of the rotor without and with an open crack passing through its first critical speed in time domain (plane A2)

Figure 3 :

 3 Figure 3: Evolutions of the orbital patterns at one-half of the first critical speed (from (a) to (h)) for the healthy rotor at plane A2

Figure 4 :

 4 Figure 4: Evolutions of the orbital patterns at one-half of the first critical speed (from (a) to (h)) for the notched rotor at plane A2 (case 4)

Figure 5 :

 5 Figure 5: Power Spectral Density in ms 2 /Hz for the healthy rotor (a) before one-half of the first critical speed (t = 40s), (b) at one-half of the first critical speed (t = 50s), (c) after one-half of the first critical speed (t = 80s), (d) at the first critical speed (t = 100s)

Figure 6 :

 6 Figure 6: Power Spectral Density in ms 2 /Hz for the notched rotor (a) before one-half of the first critical speed (t = 40s), (b) at one-half of the first critical speed (t = 50s), (c) after one-half of the first critical speed (t = 80s), (d) at the first critical speed (t = 100s)

  0

Figure 7 :

 7 Figure 7: Spectrum cascade plots of the normalized Power Spectral Density in logscale for the rotor with and without crack for both the vertical and horizontal directions at plane A2 (with an acceleration of 20rpm.s -1 )

Figure 8 :

 8 Figure 8: Transient orbital patterns for t ∈ [4 -5]s with an acceleration of 230rpm.s -1 (case 11)

Figure 10 :

 10 Figure10: Maximum amplitudes at one-half of the first critical speed and at the critical speed for various disc's unbalance (Aiv: vertical amplitude for the i th plane, Aih: horizontal amplitude for the i th )

Acknowledgments

The author would like to thank the financial support of the financial support of the french National Research Agency through the program Young researcher ANR-07-JCJC-0059-01-CSD 2, the Centre National de la Recherche Scientifique and the laboratory LTDS UMR-CNRS 5513 (Laboratoire de Tribologie et Dynamique des Systmes, Unit Mixte de Recherche du Centre National de la Recherche Scientifique 5513) to carry out this research. Some of CWT software includes code originally written by C.Torrence and G. Compo (http://paos.colorado.edu/research/wavelets/).