Touil Abderrahim
email: abderrahim.touil@lec-umc.org

Babaa Fatima

Kratz Frédéric

Bennis Ouafae

An Early Diagnosis of Outer Race Bearing Defects In Induction Machines Based On Line Neutral Voltage Analytical Calculation

This paper is devoted to diagnosis outer race bearing fault of induction motors by studying line neutral voltage analytical development. The study is based on the analysis of the frequency content of the neutral voltage after the appearance of the fault. Despite the minimal cost of bearing maintenance, the most machines failures are caused by this latter. That is why, the investigation of the outer ring bearing faults is essential due to the financial losses and break time resulting from unexpected events caused by this fault. The early detection of the bearing faults becomes necessary in order to increase the yield of the machine as well as its lifespan. Reliable monitoring can only be achieved if the chosen diagnosis method is inadequacy with the nature of the fault. The method chosen must be facilitating the diagnosis as much as possible. In fact, this mechanical defect gives very little significant defect index when using the MCSA method. The most contribution of this paper is to propose an accurate indicator of the outer race bearing fault that is different of the signature engender by the eccentricity fault. The line neutral voltage calculation taking into account the outer race bearing fault is developed and simulated. To support and to confirm the simulation results, experimental study is presented in this paper.

Introduction

Today, maintenance activities are massively oriented towards "FMDS" type services (Reliability, Maintainability, Availability, and Security). These activities permit and meet overall performance objectives to guarantee continuity of the service of the electrical systems while minimizing the maintenance costs.

The new industrial objectives are to switch to Industry 4.0 that consists of monitoring and control machinery and equipment in a real-time. For that, it is necessary to propose an optimal monitoring method that will be able to process the data provided by sensors at each stage of the production process. To obtain the fourth new industrial revolution and to have a connected factory it is necessary to studies separately and efficiency the different elements of the sub-systems of an industrial electrical installation. One of the most important constituents in the process is the induction motor.

In fact, asynchronous machine occupies an ideal position in industrial processes because of the advantage that distinguishes it to other electric machines. Its accuracy, efficiency, simplicity focus on a low maintenance cost [1][START_REF] Silva | Multi-domain model of faulty stator core for thermal effects and losses evaluation[END_REF][START_REF] Fong | An intelligent online machine fault diagnosis system[END_REF][START_REF] Mendes | Fault diagnostic algorithm for three-level neutral point clamped AC motor drives, based on the average current Park's vector[END_REF][START_REF] Alalibo | Bi-Signature optical spectroscopy for online fault detection in electrical machines[END_REF][START_REF] Rahnama | Machine-learning approach for fault detection in brushless synchronous generator using vibration signals[END_REF][START_REF] Vinayak | Wavelet-based real-time stator fault detection of inverter-fed induction motor[END_REF][START_REF] Akbari | Calculation of inductances of induction machines under axial non uniformity conditions[END_REF][START_REF] Khezzar | Induction motor diagnosis using line neutral voltage signatures[END_REF][START_REF] Khezzar | Accurate modelling of cage induction machine with analytical evaluation of inductances[END_REF].During his life, the machine is exposed to various malfunctions. For that, many extensively researchesare proposed to predict and deal well the inappropriate stop of industrial process.Induction machine malfunction can be caused by two principal kinds of failures. These failures can be outside or inside of the machine. The most important outside failure is the unbalanced supply which affects the power supply source [START_REF] Babaa | Combined Electrical Faults Detection and Diagnosis Using Current Signature Analysis[END_REF][START_REF] Morsy | Sensorless speed control of a five-phase induction machine under open-phase condition[END_REF][START_REF] Xu | SVD filtering and TLS-ESPRIT algorithm based on stator fault characteristic detection of doubly-fed induction generator[END_REF][START_REF] Babaa | Steady State Analytical Study of Stator Current Harmonic Spectrum Components on Three-Phase Induction Motor under Unbalanced Supply Voltage[END_REF][START_REF] Moosavi | Comparison of rotor electrical fault indices owing to inter-turn short circuit and unbalanced resistance in doublyfed induction generator[END_REF]. In a long way, this fault provokes a many other faults.

For the inside failures, we distinguish the electrical and mechanical faults. Thestator faults in the machine such as an inter-turn short-circuits which threatens the insulating part in the stator [START_REF] Touil | Analytical calculation of inductances under stator inter-turn short circuits fault condition in operating squirrel cage induction motors[END_REF][START_REF] Han | Modelling, diagnosis, and tolerant control of phase-to-phase fault in switched reluctance machine[END_REF][START_REF] Zeng | Inter-turn fault diagnosis of permanent magnet synchronous machine based on tooth magnetic flux analysis[END_REF][START_REF] Foitzik | Fault tolerant control of a three-phase PMSM by limiting the heat of an inter-turn fault[END_REF][START_REF] Nazemi | Stator and rotor turn-to-turn fault detection in wound rotor induction machines based on the air-gap magnetic field distortion[END_REF][START_REF] Sun | Compound frequency modulation strategy of DFIG-based wind turbines dealing with stator winding inter-turn short-circuit fault[END_REF][START_REF] Chen | Electromagneticthermal coupled modelling and analysis of inter-turn shortcircuit faults of a permanent magnet alternator[END_REF][START_REF] Alvarez-Gonzalez | Permanent magnet synchronous machine stator windings fault detection by Hilbert-Huang transform[END_REF][START_REF] Liu | Stator inter-turn fault detection in closed-loop controlled drive based on switching sideband harmonics in CMV[END_REF][START_REF] Abdallah | Stator winding inter-turn short-circuit detection in induction motors by parameter identification[END_REF][START_REF] Touil | Fault Tolerance Control and Diagnosis of Induction Machine under Stator Fault Conditions[END_REF], or a crack or rupture of the rotor bar [START_REF] Bossio | Rotor fault diagnosis in permanent magnet synchronous machine using the midpoint voltage of windings[END_REF][START_REF] Marzebali | Manipulation of stator current signature for rotor asymmetries fault diagnosis of wound rotor induction machine[END_REF][START_REF] Panagiotou | FEM approach for diagnosis of induction machines nonadjacent broken rotor bars by short-time Fourier transform spectrogram[END_REF][START_REF] Fu | Multiple coupled circuit modelling approach for squirrel cage induction machine under single-broken-bar fault with stator winding functions decomposed in d-q rotor reference frame[END_REF][START_REF] Ramu | Broken rotor bar fault detection using Hilbert transform and neural networks applied to direct torque control of induction motor drive[END_REF][START_REF] Zaggout | Detection of rotor electrical asymmetry in wind turbine doubly-fed induction generators[END_REF][START_REF] Hou | Quantitative broken rotor bar fault detection for closed-loop controlled induction motors[END_REF][START_REF] Kaikaa | Effects of the simultaneous presence of static eccentricity and broken rotor bars on the stator current of induction machine[END_REF][START_REF] Nandi | Detection of rotor slot and other eccentricity related harmonics in a threephase induction motor with different rotor cages[END_REF], or the asymmetry error due to the manufacturing.

As it concerns the mechanical faults, the rolling elements fault like the bearing fault, the ball defect and the cage of ball defect [START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF][START_REF] Jafarian | Spectral analysis for diagnosis of bearing defects in induction machine drives[END_REF][START_REF] Vilchis-Rodriguez | Wound rotor induction generator bearing fault modelling and detection using stator current analysis[END_REF][START_REF] Mcfadden | Model for the vibration produced by a single point defect in a rolling element bearing[END_REF][START_REF] Blodt | Models for bearing damage detection in induction motors using stator current monitoring[END_REF][START_REF] Schoen | Motor bearing damage detection using stator current monitoring[END_REF], are the defects that appears the most. This fault is definedas an eccentricityproduced by the axial overload on the axes of rotation [START_REF] Dorrell | Analysis of airgap flux, current, and vibration signals as a function of the combination of static and dynamic airgap eccentricity in 3-phase induction motors[END_REF][START_REF] Hadjami | Analytical model of cage induction machine dedicated tothe study of axial non-uniformities[END_REF][START_REF] Faiz | Eccentricity fault diagnosis indices for permanent magnet machines: state-of-the-art[END_REF][START_REF] Mazaheri-Tehrani | Airgap and stray magnetic flux monitoring techniques for fault diagnosis of electrical machines: An overview[END_REF][START_REF] Faiz | Detection of mixed eccentricity fault in doubly-fed induction generator based on reactive power spectrum[END_REF][START_REF] Afshar | Eccentricity fault detection in brushless doubly fed induction machines[END_REF][START_REF] Ehya | Static and dynamic eccentricity fault diagnosis of large salient pole synchronous generators by means of external magnetic field[END_REF][START_REF] Nasiri-Gheidari | New design solution for static eccentricity in single stator-single rotor axial flux induction motors[END_REF][START_REF] Faiz | Dynamic air gap asymmetry fault detection in single-sided linear induction motors[END_REF][START_REF] Oumaamar | Static air-gap eccentricity fault diagnosis using rotor slot harmonics in line neutral voltage of threephase squirrel cage induction motor[END_REF], that are divided in the outer and the inner bearing ring defects.These latterare known as a temporary static and dynamic eccentricity defects respectively [START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF][START_REF] Mcfadden | Model for the vibration produced by a single point defect in a rolling element bearing[END_REF]. Bearings are a very most important part of the machine because they are the force transmission elements between the moving and static components. As such, they contribute considerably to the proper functioning of equipment. Bearings deteriorate as a result of manufacturing errors, incorrect assembly, static and dynamic loads, deficient lubrication or because of the harshness of the environment. However, even if a bearing is perfectly manufactured, assembled, etc., it will eventually present a defect caused simply by normal wear of the elements that compose it.Studying bearings permits to understand their shaping in order to predict their operation and make a distinction between theimpacts on the behaviour of the machine and also understand the signatures allowing going back to the main cause of the failure. The one of the first works to have studied this defect is McFadden and al. [START_REF] Mcfadden | Model for the vibration produced by a single point defect in a rolling element bearing[END_REF]. This paper is considered as a key reference in bearing fault detection. In fact, bearing faults provoke the misalignment and eccentricity in the rotating part that can be divided into several parts, static eccentricity, dynamic eccentricity, but in the reality, we have the co-existence of the both faults in the same time which called by the mixed eccentricity.

There are various diagnostic methods to monitor the conditions of the working system and to prevent the industrial system from being interrupted [START_REF] Touil | Analytical calculation of inductances under stator inter-turn short circuits fault condition in operating squirrel cage induction motors[END_REF][START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF].

Depending on the measurable quantities, the diagnostic methods know is: Stator current analysis, vibration analysis, Flux analysis, Electromagnetic torque analysis and neutral to neutral analysis [START_REF] Khezzar | Induction motor diagnosis using line neutral voltage signatures[END_REF].

Line neutral voltage is defined as the voltage present between the neutral of the power source and the neutral of theasynchronous motor.

In this paper we study the outer ring bearing fault by the line neutral analysis method. Through an analytical development of line neutral voltage, the paper proposes also to separate and dissociate the indices of the apparition of this fault with the others related to the eccentricity fault. We notice that, in previous works that, line neutral voltage has been previously used and that it gave better results in comparison with the Motor Current Signature Analysis method (MCSA) [START_REF] Oumaamar | Static air-gap eccentricity fault diagnosis using rotor slot harmonics in line neutral voltage of threephase squirrel cage induction motor[END_REF].The development of our model is based on the sensitivity of the magnetic field in the air gap which separates the stator from the rotor. In fact, the variation of this air gap as a function of time will be changed in presence of fault and stay unchanged in healthy conditions. In some fault condition, the minimum thickness of the air gap will vary over time. This variation influences considerably the values of the different inductances of the asynchronous machine [START_REF] Touil | Analytical calculation of inductances under stator inter-turn short circuits fault condition in operating squirrel cage induction motors[END_REF][START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF], in particular the mutual inductances which depends on the permeability of the air gap. This analytical study takes into account the variation that appear in the mutual inductances and propose a good and accuracy formula,that is using to the inductances calculation.To support our study and the proposed analytical development, experimental results are also presented in the paper.

Analytical multi-loop model for squirrel cage induction motor.

Electrical Equations

In our squirrel-cage induction motor modelization, we assume that the distribution of the windings, in the stator phases and the rotor cage, is identical.

We also suppose that bars are connected electrically with the end-rings to produces an identical spaced multiloop. Under these assumptions, the system equations in vector matrix form, including line neutral voltage is:

Where:

, and are the stator current, voltage and flux vectors respectively.

, and are the rotor current, voltage and flux vectors respectively, and for , and are the correspondent end ring voltage, current and flux.

, are the stator and rotor resistance matrices respectively.

, are the stator and rotor inductances respectively.

For the neural voltage there is whith is placed between the supply neutral set point and the stator neutral set point.

As there is no way to quantify and to know the behaviour of line neutral voltage, we put the previous equation in the line-to-line voltage form:

Where:

, and From with .

Mechanical Equations

According to the principle of the energy conversion and considering change of co-energy of the system, the torque variation taking into account a slight change in the rotor's position, when the currents are held constant is:

The equations describing theelectromagnetic torque, the rotor's position and the rotor's angular speed as follow: The angular speed of the rotor can be calculated from the variation of the electromagnetic torque with the presence of a load as in equation [START_REF] Babaa | Combined Electrical Faults Detection and Diagnosis Using Current Signature Analysis[END_REF].

The minimum length position of the radial air

gap assumption Starting from the common assumption that the air gap is constant in the construction of the asynchronous machine, the movement of the rotor inside the machine with the adequacy of the rotor axis is better. In healthy monitoring conditions,in the proposed model, we have chosen is the axis of rotation. Inthe case, of any abnormal factors such as poor manufacturing, nonuniformly circular stator core or the bearings fault, we remarkthe existence of a permanent static eccentricity. Most referenced researches about static eccentricity rises the idea of a misalignment between the rotor centre and the axes of rotation.Regardless of mixed eccentricity, that is the most famous one among eccentricity faults, the pure static eccentricity is defined as a fixed position of the minimal radial air-gap length in space(see Fig. 1)expressed with the following formula as:

When the outer race bearing fault occurs in presence of a small gap or slit as in Fig. 1,the rolling elements pass by the irregularityandproduces a temporary mismatch between the rotor centre and the axes of rotation. This phenomenon leads to a vibration considered as temporary static eccentricity as in [START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF]. This vibration signal appears with a periodic frequency called The Bale Pass Frequency (BPFO) (see Fig. 2) which is described as a train of impulses. This signal is expressed by: Fig. 3 shows the shape of the length of the air gap as a function of the angle of rotation in different states, healthy condition which is constant, in the static eccentricity condition referenced to equation ( 13)and in the outer ring bearing fault with equation ( 15).

Inductances calculation

Dimensioning inductors theoretically consists in the study of the winding function theory which is the development of mutual inductance between two windings in the electric machines. -a-

-b-

In healthy state the air gap is uniform, sothe winding stator functions is calculated as follows:

To get a reality of the machine, it is necessary to considering the non-uniform of the air gapwhen we have a faulty condition. Then it is necessary to modify the winding function.By using the formula of the winding function theory [START_REF] Khezzar | Accurate modelling of cage induction machine with analytical evaluation of inductances[END_REF][START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF].The different inductance calculations (stator, rotor, mutual stator-rotor) depend h val f "p" mb f p l pa W our case the condition p Inprevious research we calculate the difference inductors under healthy state and both of the faulty parts (static eccentricity and outer race bearing fault) [START_REF] Abderrahim | Analytical Model for Separated Frequency Signature of Outer Race Bearing Fault from Static Eccentricity[END_REF].

The choice to work with the presence of static eccentricity is due to the proposed approach that the condition for the occurrence of bearing outer ring faults is a temporary fault which affects the performance of the machine instead of a permanent fault when a static eccentricity occurs. The changes that occur in inductance values led to the conclusion that there is a small change in stator coils values and of air gap length change affects ' l .On the other hand, the air gap length variation affects the distribution of the mutual induction by causing an imbalance in the distribution of the magnetic field around the rotor.

So, the mutual inductance between stator coils and rotor loops with the consideration of space harmonics is given by: Healthy condition:

Static air gap eccentricity condition:

Where:

Outer race bearing fault condition:

Fig. 4 shows different presentations of mutual inductance in healthy and faulty conditions. The influence of the outer bearing fault is considered as a mp a y a y ma k by a fa l a l γ and a static eccentricity level .

Line neutral voltage

Healthy Condition

As known, the multi-loop model of induction machine supplied by an asymmetrical three-phase network voltage produces a rotating field and causes the creation of current in the rotor loops resulting the rotation of the shortcircuit cage.

The multi loops currents expressionis: Where:

The expression of line neutral voltage is a two-part formula, withshape of a sinusoid that is the sum of an N b series of equations shifted by This expression equals zero except when the space harmonics on the part of the rotor cage have a common point with . This assumption leads to a specified setwhich represents the order of a high frequencies component in the induction machine line neutral voltage.

Under this condition sets, the harmonics components that appear in the line neutral voltage are a consequence of rotor loops distribution called the RSHs. These harmonics and each related orders are shown in Table 1. Fig. 5 represents the spectral content in for line neutral voltage for a healthy machine.The simulation was made for an induction machine whichwere made for a 1.1-kW motor with 24 stator and22rotor slots, 2 pair poles.

The RSHs which can be seen onthe spectrum curve at related to the order h b with and , the second one at related to the orderh b with . It is evident that for the case of or theRSHs harmonics doesn't appear because of the correspondent harmonic order does not belong to 'G'.

Under Static Eccentricity Action

Permanent misalignment of the axis of rotation caused by an unequal distribution of the air gap between the stator and the rotor leads to the existence of field variation that crosses the air gap cause by mismatch of the minimal space. Also, a permanent misalignment of rotation axis caused by a small eccentricity affects the inductances values which are dependent on the winding and the permeance functions (the inverse air gap function).In this case the change in the flux which crosses the rotor introduces a change in currents of rotor loops. The resulting currents are defined by the interaction of the stator current with the new part of the mutual inductance.

The rotor loops current become:

The expression for line neutral voltage:

Where: is the additional part of the mutual induction that defines the eccentricity appearance.

The first term represents the interaction of healthy mutual inductance with the healthy rotor loops current.The first additional part of line neutral voltage term is related to the interaction of healthy mutual inductance with the additional rotor loops current that is caused by the eccentricity appearance describe in [START_REF] Hadjami | Analytical model of cage induction machine dedicated tothe study of axial non-uniformities[END_REF].

Where:

Thisexpression is a sumof multiple waves shifted either by or by . Equation(37)equals zero except when the space harmonicsof the rotor cagepart have a common point with . Undid, in this casethe assumption leads to a specified order of harmonics due to static air gap eccentricity related to space harmonics with the following set:

The second additional part of line neutral voltage term is related to the interaction of the faulty part of the mutual inductance with the healthy part of rotor loops current (the term is described in [START_REF] Mazaheri-Tehrani | Airgap and stray magnetic flux monitoring techniques for fault diagnosis of electrical machines: An overview[END_REF]).

Where:

hp hp
This expression is a sumof multiple waves shifted either by hp or by hp . Equation(39)equals zero except when the space harmonicsof the rotor cagepart have a common point with . In this casethe assumption leads to a specified order of harmonics due to static air gap eccentricity related to space harmonics with the following set:

The third additional part of line neutral voltage term is related to the interaction of faulty part of the mutual inductance with the additional rotor loops current caused by the eccentricity appearance, describe in [START_REF] Afshar | Eccentricity fault detection in brushless doubly fed induction machines[END_REF].

Where:

hp hp hp hp
This expression is a sumof multiple waves shifted either by hp or by hp 1 .

Equation ( 41)equals zero except when the space harmonicsof the rotor cagepart have a common point with . In this casethe assumption leads to a specified order of harmonics due to static air gap eccentricity related to space harmonics presentedin the following set:

hp hp hp hp
Under this condition sets, the harmonics components that appears in the line neutral voltage under static eccentricity fault are shown in Table 2.

The spectral content for line neutral voltage of faulty machine (see Fig. 6), and the RSHs which can be seen additional to the ones in healthy condition is according to assumption the ones in the third additional term. So, the spectrum curve at related to the order h b with and , the second onesat related to the order h b and related to the orderh b . For the case of =2 or =4 the RSH harmonics appear in the faulty case with the same corresponding harmonic order 'G' as in healthy condition.

Presence of outer race bearing fault

Since theouter race bearing fault is the fault that appears by any cracks or deformation or breakage, or irregularity in the outer race of the ball bearing, each contact between the ball and the outer ring produces a temporary vibration can be considered as temporary static eccentricity. The change of the flux which crosses the rotor provokes a change in currents of rotor loops.Then, anew resulting current defined by: The line neutral voltage expression is:

Where: In line neutral voltage expression two terms are exist. The first is related of healthy state, the second related to the static eccentricity. However, another term appears. This latter is related to the bearing fault of the outer ring from the interaction of healthy mutual inductance with the additional rotor loop current. It is also from the interactionof the additional part of the mutual inductance caused by the occurrence of the fault with the healthy currenttaking into account the interaction with the static eccentricity expressions.

The first additional part of line neutral voltage term is describedin [START_REF] Jafarian | Spectral analysis for diagnosis of bearing defects in induction machine drives[END_REF].

Where:

The expression is a sumof multiple waves shifted either by or by . This equation equals zero except when the space harmonicsof the rotor cagepart have a common point with . The assumption made above leads to define the order of harmonics generated by the bearing faults:

The second additional part of line neutral voltage term is related to the interaction of the faulty part of the mutual inductance caused by bearing fault appearance with the healthy part of rotor loops current which is describe in [START_REF] Mcfadden | Model for the vibration produced by a single point defect in a rolling element bearing[END_REF].

Where:

hp hp
The expression is a sumof multiple waves shifted either by hp p or by hp p . This equation ( 49) equals zero except when the space harmonicsof the rotor cagepart have a common point with . So, theorder of harmonics due to bearing fault appearance related to space harmonics is:

The third additional part of line neutral voltage term is related to the interaction of faulty part of the mutual inductance caused by the eccentricity appearance with the additional rotor loops current caused by bearing fault appearance which is describe in [START_REF] Schoen | Motor bearing damage detection using stator current monitoring[END_REF].

The expression is a sumof multiple waves shifted either by hp or by hp ±1 .

This equation ( 51) equals zero except when the space harmonicsof the rotor cagepart have a common point with . In this case, the orders of harmonics are:

As it concerns the forth additional part of line neutral voltage term, is define as the interaction of faulty part of the mutual inductance caused by bearing fault appearance with the additional rotor loops current caused by the eccentricity appearance which is describe in (53).

The expression is a sumof multiple waves shifted either by hp or by hp ±1 .

This equation equals zero except when the space harmonicsof the rotor cagepart have a common point with . In this case theorder of harmonics is define as:

The last additional part of line neutral voltage term is related to the interaction of faulty part of the mutual inductance with the additional rotor loops current caused by bearing fault appearance which is described in (55).

Each part of this expression is a sum of multiple waves shifted either by hp or by hp . The equation (55) equals zero except when the space harmonicsof the rotor cagepart or have a common point with . We put theorder of harmonics in the following equation.

The harmonics components that appear in the line neutral voltage under outer race bearing fault are shown in Table 3. The RSH orders for the frequencies engender by outer race bearing fault indicate that only the third and the forth terms of line neutral voltage contribute to the appearance of some harmonics related to RSH components with the formula f The fifth term contribute with the appearance of the RSH components define by the formula f .

Simulation results

For simulation, we use a machine with 2 pairspoles , 3 phases, 50Hz with 22 rotor bars.The spectrum results were computed with Fast Fourier Transform, all results were obtained with MATLAB software.The line neutral voltage spectrum is simulated in static eccentricity fault condition. We confirm the appearance of several RSHs related to the fault. Another RSHs harmonicsare also evident and are related tothe healthy condition.By using [START_REF] Panagiotou | FEM approach for diagnosis of induction machines nonadjacent broken rotor bars by short-time Fourier transform spectrogram[END_REF]healthy condition contributes to the emergence of RSHs with the orderh b p andh b p respectively when and with , While and do not appear, they indicate that the harmonic order does not belong to "G".

The RSHs exist for and due to the presence of SE.The harmonics with the order h b p andh b p belong to " ".

Table 2 summarize the frequencies shown on line neutral voltage spectrum in the case of SE existence.Under the same testing conditions with the slip the RSHs that can be seen additionally to the othersare related to the orderh b with , the second oneat related to the orderh b and related to the orderh b . In bearing fault condition, the sum of the frequencies which appear can be divided into three parts. The first one when harmonic order belongs to "G", that is the harmonics of the healthy part. These frequencies appear also in spectrum curve of outer race bearing faultTable 1.

The second part is related to the RSHs with and these latter are appeareddue to the presence of a temporary SE when the ball passes by the irregularity and the harmonic order belong to " "Table 2.

The third part, harmonics represent the RSHs related to bearing fault exclusively.The harmonic order belongs to " " or " "with an offset of , and with when harmonic order belongs tothe association of" " and " " Table 3.

h f h f h f h f h f h f h f h f h f f h f f h f f h f f 7.

Experimental results and discussion

A test bench has been set up at Constantine Electrical Engineering Laboratory, "LEC"University Constantine 1, 25000, Constantine,Algeria.In order to support the suggested model for outer ring bearing fault diagnosis using line neutral voltage, we synthesize data from experimental tests with a 5-kW squirrel-cage induction motor. The motor is loaded with magnetic powder brake. The tested motor is a 5 KW three-phase induction motor with28 rotor bars, 2 pole pairs of FIMET manufactures, with accessible neutral point(see Fig. 8). The line neutral voltage acquisition has been accomplished using dSPACE with (50 kHz) sample rate and data signal length of (5s).

The outer ring defect was created by grooving the outer ring of the bearing(see Fig. 9). The results obtained for the outer ring bearing fault areobtain after several tests for several load speeds of 1421 rpm and 1422rpm for a defective bearing and 1420 rpm /min for a healthy bearing -Series 6208 bearing (Bd= 10.49 mm, Pd= 61 mm, N=9).

This type of bearing series feature is low friction and designed to generate low noise and vibration levels at high rotational speeds. They support radial and axial loads in both directions, are easy to mount, and require little maintenance.

We made a diagnosis by analyzing the line neutralvoltage. The line neutral voltage spectrum was studying and analyzefor thehealthy and the defective cases. One can already notice spectral components at f s ± f r in the healthy state because of the original eccentricity. In the faulty case, the component f s + f r shows a tiny growth in value, especially for the first harmonic as shown in (see Fig. 10).

Like described in many research before, these components are related to the mixed eccentricity that exist on the sideband of the low and high frequencies with the relation f s + f r. They appear in a healthy state as a consequence of a small vibration caused by the asymmetry in rotor construction and load torque oscillation.In outer race bearing fault,we clearly notice the increase of these components.

The machine being directly connected and tested through the power source,an unbalanced exist at the beginning. Following this, some specific harmonics exist in both healthy state and outer race bearing fault. These harmonics exist in all induction motors situation due to an existence of a small unbalance provide by the means supply (see Fig. 11).

The spectral analysis of the line neutral voltage always shows distinct frequencies linked to the distribution structure of the rotor bars(Harmonics related to RSHs). Indeed, these frequencies still exist but with low perceptible eccentricity(see Fig. 12), which explains their existence in the outer race bearing fault condition, as shown in [START_REF] Ehya | Static and dynamic eccentricity fault diagnosis of large salient pole synchronous generators by means of external magnetic field[END_REF].

We remark, the same harmonics appear in a healthy state because of the presence of static eccentricity due to a small asymmetry in the induction machine construction.

To find the different bearing fault signatures of the outer ring, it is necessary to identify the harmonics using the formulas in Table 3. Fig. 13 shows some of specific harmonics amplitudes related to the outer race bearing fault have increased due to the small eccentricity resulting by the ball movement and the irregularity in the outer ring in fully loaded machine condition «see Appendix 1 for table of results". Fig. 14 represents some harmonics predicted in literature related to f 0 the characteristic frequency of bearing outer ring failureTable 4. The tests are run several times and the induction motor is tested under the same load conditions, with different fault level. The specific harmonics of the outer ring fault increase with increasing fault level. This increase is justified by the increase in degree of eccentricity. Fig. 15 shows some examples explaining the differences.

Conclusion

To make a good diagnosis, there are several methodswhich have been explored and tested in the past, with varying degrees of effectiveness. All methods aim to provide a reliable and accurate fault indicator that can detect, at early stage, any fault that occurs in the machine.

Asophisticated monitoring method is vitally important to ensurean optimal maintenance operation of the machine. Thus, it is necessary to provide a reliable and easily interpretable diagnosisto detect the appearance and evolution of most mechanical and electrical faults. This early detection allowsavoid unplanned production stoppages.Studying of outer race bearing fault occupies a most important consideration in the new researches. It is necessary to have a strongly information of the defecting part of the machine to propose an adequate diagnosis method. The MCSA method is very poor in information for the outer bearing faults. The proposed signature is not clear and it is difficult to interpret it.

In this paper an analytical development of line neutral voltage is proposed. The developed model takes into account the healthy and outer bearing fault conditions. The simulation results of our equations show that the neutral voltage contains very rich information about this fault. Following the analytical results obtained, a reliable and specific bearingdefect index has been proposed.The proposed index permits also to separate the outer race bearing fault to another eccentricity faults. To confirm our analytical development, a series of experimental results are presented and discuss.

Fig. 3 .

 3 Fig. 3.Air gap length variation in different conditions.

Fig. 1 .

 1 Fig. 1. Air gap length in different conditions (a) Healthy (left)-static eccentricity SE (right), (b) Outer race bearing fault

Fig. 2 .

 2 Fig. 2.Vibration signal (train of impulses).

Fig. 4 .

 4 Fig. 4. Mutual inductance Msr, for Healthy situation, for static eccentricity and outer race bearing fault.

Fig. 5 .

 5 Fig. 5.Simulated, FFT spectrum of line neutral voltagein healthy induction machine state.

Fig. 6 .

 6 Fig. 6.Simulated, FFT spectrum of line neutral voltageunder static eccentricity induction machine state.

Fig. 7 .

 7 Fig. 7.Normalised FFT spectrum of the line neutral voltage for deferent states on induction machine performance. (Top) Healthy, static eccentricity SE (middle), (bottom) Outer race bearing fault ORBF.Table 3The frequencies shown on line neutral voltage

Fig.

  Fig. 8.Test bench dedicated to the study of outer ring bearing fault diagnosis using line neutral voltage

Fig. 12 .

 12 Fig. 12. Experimental FFT spectrum of line neutral voltage (Harmonics related to RSHs).

Fig. 11 .

 11 Fig. 11.Experimental FFT spectrum of line neutral voltage (Harmonics related to small unbalance provide by the means supply).

Fig. 10 .

 10 Fig. 10.Experimental FFT spectrum of line neutral voltage (Harmonics related to the component f s + f r).

Fig. 9 .

 9 Fig. 9. Some bearings for the experiments

[ 1 ]

 1 Khan, F., Sulaiman, E., & Ahmad, M. Z.: 'Review of switched flux wound-field machines technology', IETE Technical Review, 2017, 34, (4), pp. 343-352

Fig. 15 .

 15 Fig. 15. Deferent angle of defected outer ring (A1>A2).

Fig. 14 .Fig. 13 .

 1413 Fig. 14.Experimental FFT spectrum of line neutral voltage (Harmonics related to f 0 )

Table 1

 1

		Rotor slot harmonics for the healthy induction
	machine in the line neutral voltage
	RSHs orders		RSHs frequencies
	h	b p	h	b p	f
	h	b p	h	b p	f

Table 2

 2

		The frequencies shown on line neutral voltage
	spectrum in the case of SE existence
	Neural	RSHs orders	RSHs frequencies
	voltage
		h	h	f
		h	h	f
		h	h	f
		h	h	f
		h
			h	f
			h	f

Table 3The

 3The

		frequencies shown on line neutral voltage
	spectrum in the case of ORBF existence
	Neural	RSHs orders	RSHs related frequencies
	voltage

  8.Test bench dedicated to the study of outer ring bearing fault diagnosis using line neutral voltage

					h					hp
					h		hp
				h				hp
			h			hp
				h			hp
			h					hp
				h			hp
			h				hp
					h			hp
			h				hp
				h			hp
			h				hp
				h					hp
			h				hp
				h	hp		p	hp		γ
	hp	p	h	ka		hp
			hp		h	p			γ	hp hp	p	ka
					h		hp
				h	hp			p hp			γ
	hp	p	h	ka		hp
			hp h	h		p	hp	γ hp	hp	p	ka
				h			hp
			h		hp		hp	p			γ
	hp	p		ka	h			hp
			hp h			p	hp		γ	hp	p	ka
				h			hp
			h	h	hp		hp	p	hp		γ
	hp	p	hp h	ka h		p	hp	hp	γ	hp	p	ka
			h			hp

Table 4

 4 Harmonics related to f 0 the characteristic frequency of bearing outer ring failure

	Outer race bearing failure characteristic
	harmonics.
	87,36
	174,71
	262,07
	349,43
	436,78
	786,21
	1135,64
	1222,99

Table 6Table

 6Table

		of results
		f
		k
		k
		k
		k
	Table 5Table of results	k
	f
		k
	k
		k
	k
	k
		k
	k
		k
	k
	k


Olivier Rousseau

Yves Bourgault

HEART SEGMENTATION WITH AN ITERATIVE CHAN-VESE ALGORITHM

This paper presents 2D and 3D applications of the Chan-Vese model to heart and trachea segmentation. We improved the multi-phase Chan-Vese model by introducing an iterative method, by choosing an appropriate L 1 fidelity term as well as an efficient and prior free initial condition. For 3D applications, the algorithm is parallelized in order to speed up the computations. We provide extensive information on computational details, on the convergence times and on the quality of segmentations obtained. The results of the segmentations are then meshed to be used for finite element simulations.

Introduction

The recent increase of computers capacities now allows for realistic simulations of human organ physiology. Mathematical modeling of organ functions opens a wide range of research, allowing diagnostics and understanding of malfunctions and diseases. Such simulations usually require a mesh of the given organ. However, most computations are made on meshes of idealized geometries and there is a real lack of accurate 3D models. Realistic geometries should be extracted from medical images, this is known as the segmentation process.

The goal of this work is to segment the heart muscle from high resolution CT scans of the thorax and to produce meshes that are adequate for numerical simulations in electro-physiology.

Most existing methods for heart segmentation involve a prior knowledge of the heart's shape (see for example [START_REF] Zhukov | Dynamic deformable models for 3D MRI heart segmentation[END_REF][START_REF] Ecabert | Towards automatic full heart segmentation in computed-tomography images[END_REF]). We intend to use a modified version of the Chan-Vese model [START_REF] Tf Chan | Active contours without edges[END_REF], also known as Active Contours without Edges to segment the heart, since it has no geometrical or topological a priori. This method has been successfully used for brain segmentation [START_REF] Vese | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF][START_REF] Chung | Energy minimization based segmentation and denoising using a multilayer level set approach[END_REF]. Brain images usually have nice contrasts between gray and white matter. To our knowledge, this method has not been attempted yet for segmenting 3D scans of the heart and lungs. This is a real challenge since the images are more diversified. They contain many objects of similar grey levels that the method should separate.

The main contributions of this paper are:

(1) Application of the Chan-Vese algorithm to trachea and heart segmentation. (4) Analysis of parallelization of the Chan-Vese algorithm to fit needs of large applications. (5) Creation of meshes from the segmentation results.

In the first section, the Mumford-Shah functional and the Chan-Vese model are presented as well as the numerical methods for solving the latter. In section 2, the improvements made to the original Chan-Vese model are explained and justified. Sections 3 and 4 respectively present 2D and 3D applications of this modified Chan-Vese algorithm. where Ω is some region of R n , typically a square or a cube. In the segmentation process, the goal is to split the image g into its constituting objects. Mumford and Shah proposed to minimize the following functional

(1) E MS (u, K) = Ω |∇u| 2 dx + λ Ω\K |g -u| 2 dx + µH N -1 (K)
over pairs (u, K). K is a compact subset of Ω representing edges of objects in g, and u ∈ H 1 (Ω \ K) is the intensity of the image. This intensity varies smoothly inside the connected components of Ω \ K [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. The middle term in the equation is called the fidelity term, while the combination of the two others form the regularity terms that do not depend on the underlying image g. The weights µ and λ should be adjusted in accordance with the noise level of the image to be segmented. It is a well known result that minimizers of the Mumford-Shah energy exist (see for example [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]). It is a hard task to directly find the actual minima of the functional. To achieve this, one can approximate the Mumford-Shah functional by a sequence of elliptic functionals easier to solve [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF][START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]. An other approach is to minimize E M S over a restricted domain. One example is the Chan-Vese model. It seeks for a minimum over functions that take only two values c 1 and c 2 [START_REF] Tf Chan | Active contours without edges[END_REF].

Such functions can be written as

u = c 1 1 F + c 2 (1 -1 F ) for some F ⊆ Ω. The Mumford-Shah energy then rewrites as (2) E CV (F, c 1 , c 2 ) = F |g -c 1 | 2 dx + Ω\F |g -c 2 | 2 dx + H N -1 (∂F ∩ Ω).
The set F can be described via a level set function φ : Ω→ R, that is F = {φ ≤ 0}.

Then the Level Set formulation [START_REF] Osher | Fronts propagating with curvature-dependent speed-Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Tf Chan | Active contours without edges[END_REF] of the problem 2 is to find φ that minimizes (3)

E CV (φ, c 1 , c 2 ) = λ Ω |g-c 1 | 2 (1-H(φ)) dx+λ Ω |g-c 2 | 2 H(φ) dx+µ Ω |DH(φ)| dx,
where H(•) stands for the Heaviside function. |DH(φ)| = δ(φ)|∇φ| is the derivative of H(φ) in the sense of distributions. This is a Dirac measure with support on the discontinuity set {x : φ(x) = 0} of u. The two-phase image then rewrites as

u = c 1 (1 -H(φ)) + c 2 H(φ).
We can deduce that if (φ, c 1 , c 2 ) is a minimum of E CV , then it is one of the many solutions of the Euler-Lagrange equation. To avoid difficulties related to the non-uniqueness of φ, a time dependence is added and the following initial problem is solved:

(4)        ∂φ ∂t = δ(φ) µdiv ∇φ |∇φ| + λ (g -c 1 ) 2 -(g -c 2 ) 2 φ(x, 0) = φ 0 (x), ∂φ ∂n = 0 on ∂Ω,
for some initial condition φ 0 which describes an initial curve C 0 = {φ 0 = 0}. The values c 1 and c 2 are also given by ( 5)

c 1 = Ω g(1 -H(φ)) dx Ω (1 -H(φ)) dx , c 2 = Ω gH(φ) dx Ω H(φ) dx ,
That is, c 1 and c 2 are mean values of g inside the regions {φ < 0} and {φ ≥ 0} respectively.

While solved with a time-stepping scheme, the problem 4 can be interpreted as a gradient method applied to the minimization of the functional in 3. Any steady state of 4 will satisfy the Euler-Lagrange equation.

In order to get a splitting of the image into more than 2 phases (multi-phase segmentation), Chan and Vese proposed to have several level set functions φ i [START_REF] Vese | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF].

A set of n level set functions can describe up to 2 n different regions of Ω. In this case n coupled PDE's must be solved simultaneously. An other approach is to use several level curves of the level set function [START_REF] Chung | Energy minimization based segmentation and denoising using a multilayer level set approach[END_REF]. In this case the number of regions has to be chosen in advance. This method does not allow to have junctions of three different regions as level curves are parallel to each other. 1.2. Numerics. The equation 4 can be solved via standard finite difference scheme ( [START_REF] Tf Chan | Active contours without edges[END_REF], [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]) using a C ∞ regularization δ of δ. For example, in 2D the first term of div ∇φ |∇φ| can be discretized as

(6)   ∂ ∂x φ x φ 2 x + φ 2 y   i,j ≈ D - x     D + x φ i,j (D + x φ i,j ) 2 + D 0 y φi,j +D 0 y φi+1,j 2 2 + ε 2     . D + x φ, D 0 x , D -
x respectively stand for forward, centered and backward finite difference approximations. Hence an explicit discretization of equation 4 can be written as

φ n+1 i,j = φ n i,j + ∆tδ (φ i,j )[d 1 (φ n i+1,j -φ n i,j ) + d 2 (φ n i,j -φ n i-1,j ) + d 3 (φ n i,j+1 -φ n i,j ) + d 4 (φ n i,j -φ n i,j-1 ) (7) + λ (g -c 1 ) 2 -(g -c 2 ) 2 ].
This is the form of a non linear diffusion with coefficients

d 1 = 1 h 2 x µ (D + x φ i,j ) 2 + D 0 y φi,j +D 0 y φi+1,j 2 2 + ε 2 d 2 = 1 h 2 x µ (D - x φ i,j ) 2 + D 0 y φi,j +D 0 y φi-1,j 2 2 + ε 2 (8) d 3 = 1 h 2 y µ D 0 x φi,j +D 0 x φi,j+1 2 2 + (D + y φ i,j ) 2 + ε 2 d 4 = 1 h 2 y µ D 0 x φi,j +D 0 x φi,j-1 2 2 + (D - y φ i,j ) 2 + ε 2
However, for an explicit time discretization, the time step constraint for stability may be quite restrictive, leading to a large number of time steps to reach a steady solution.

It has been noticed that replacing δ(φ) by |∇φ| increases the stability of the PDE [START_REF] Marquina | Explicit algorithms for a new time dependent model based on level set motion for nonlinear deblurring and noise removal[END_REF], yielding [START_REF] Tf Chan | Active contours without edges[END_REF]

φ t = |∇φ| µdiv ∇φ |∇φ| + λ (g -c 1 ) 2 -(g -c 2 ) 2 ,
which can be interpreted as a motion of the interface in the normal direction with speed V n = -µdiv ∇φ |∇φ| -λ (g -c 1 ) 2 + (g -c 2 ) 2 . Note that equation 9 does not admit any steady state, although the interface {φ = 0} converges to the zero level set of the steady state of equation 4. In fact when equation 9 evolves over time, the level set function φ will converge to -∞ inside the curve and to +∞ outside. To avoid dealing with too large values, φ can be truncated, for example restricted to [-100, 100]. Althought, on must be careful as we remarked that if φ is truncated to a small value (eg. to [-1, 1]), it may lead to numerical artifacts. It is also possible to use a semi-implicit time discretization, which also increases the stability [START_REF] Smereka | Semi-Implicit Level Set Methods for Curvature and Surface Diffusion Motion[END_REF]. In this scheme, the value φ i,j of the central pixel is taken at time n + 1 for the right hand side of equation 7. This yields

φ n+1 i,j = φ i,j + ∆t|∇φ| i,j [d 1 φ n i+1,j + d 2 φ n i-1,j + d 3 φ n i,j+1 + d 4 φ n i,j-1 + (g -c 1 ) 2 -(g -c 2 ) 2 ] 1 + ∆t|∇φ|(d 1 + d 2 + d 3 + d 4 )
.

(

) 10
We will analyze the convergence of equation 9 to a steady state using the semiimplicit time discretization described in equation 10 under several factors. We will study the influence of the initial condition and of the fidelity term chosen. We will also propose an iterative method to obtain multi-phase segmentation.

In [START_REF] Song | A fast algorithm for level set based optimization[END_REF], Song and Chan decoupled equation 9 in order to simplify the resolution. They first smooth the image g to obtain a new image g * , and then solve the simplified equation ( 11)

φ t = (g * -c 1 ) 2 -(g * -c 2 ) 2 .
They proposed a fast algorithm for solving equation 11. In fact, their algorithm corresponds closely to clustering the set {g * (x) : x ∈ Ω} of values of g * into two clusters using a k-means procedure [START_REF] Hartigan | Clustering Algorithms[END_REF]. Osher and He introduced a similar algorithm to solve the multi-phase problem [START_REF] He | Solving the Chan-Vese Model by a Multiphase Level Set Algorithm Based on the Topological Derivative[END_REF]. Again, this is closely related to applying a k-means procedure to the set {g * (x) : x ∈ Ω} with the right number of clusters k. Independently, Gibou and Fedkiw took advantage of this analogy to propose an hybrid algorithm that alternates between smoothing and k-means [START_REF] Gibou | A fast hybrid k-means level set algorithm for segmentation[END_REF]. However, it is not known if any of these fast algorithms lead to minimums of the Chan-Vese energy for some value of the parameters µ and λ.

Improvements

We focus on the resolution of the Chan-Vese model via equation 9 in order to extract the heart shape from a 3D CT scan of a human thorax. We propose improvements on several aspects of the resolution process in order to fit the needs of large 3D applications. These aspects are:

(1) Iterative segmentation instead of multi-phase segmentation.

(2) Choice of a fidelity term in equation 9.

(3) Choice of an initial condition in equation 9.

2.1. Iterative Segmentation. Medical images usually contain many regions of different pixel intensity value. In this case a two-phase segmentation will generally not be able to extract the region of interest.

A solution to the multi-region problem is the algorithm proposed by Chan and Vese for multi-phase segmentation that simultaneously evolves several level set functions, that is several curves [START_REF] Vese | A Multiphase Level Set Framework for Image Segmentation Using the Mumford and Shah Model[END_REF]. There are drawbacks to this method. One is that instability may appear from solving these PDE's simultaneously. Also it may happen that the two curves need to coincide at some places. In this case, if the curvature term is dominant, it may lead to miss-classification in this region, since the curves may not superimpose correctly. Figure 1 shows an example of this phenomena.

The other approach is to look at several level curves of the level set function φ [START_REF] Chung | Energy minimization based segmentation and denoising using a multilayer level set approach[END_REF].

If n different values are chosen, it splits the domain into n+1 different regions. One possible issue is that it can not recover triple junctions, as level curve are naturally parallel to each other.

Instead of these strategies, we will segment by solving the simple two-phase Chan-Vese problem iteratively. The process is as follows:

(1) Split the domain Ω into Ω + and Ω -using the 2-phase Chan-Vese model.

(2) Stop if the object is extracted, that is if the object is either Ω + or Ω -.

Otherwise, decide which of Ω + or Ω -contains the object of interest, and pick this region as a new domain for step 1.

Figure 2 shows the results of an iterative segmentation process on a slice of a CT scan of the heart.

There are many advantages to proceed in this manner. First, if we are interested in a single object in the image, the computations need not to be done on the whole domain after the first step. This clearly saves on CPU time. We can thus focus on the region of interest, rejecting at each step the part of the image that does not contain the given object. However a full 2 n -phase segmentation of the image can also be done with this iterative method. The first step splits the domain into two different sub-domains Ω + and Ω -. In the second step, each of these sub-domains will be split in 2 parts, yielding four sub-domains Ω ++ , Ω +-, Ω -+ and Ω --. This can be done until the requested number of sub-domains is reached. It takes n such steps to get 2 n regions. In the multi-phase case, n different level sets functions need to be evolved. One step of the iterative method has the same computational cost as evolving a single level set function in the multi-phase method. However in practice the multi-phase algorithm requires smaller time step to ensure stability, which makes it slower. Another advantage of the iterative method is that the smoothing parameters can be varied from steps to steps. This saves time since the equation can be solved with larger time steps when the smoothing term is not dominant. Another benefit of our iterative method is that triple junctions are accurately obtained since the boundary conditions force the curve to be normal to the domains boundary. Figure 3 shows how our method performs on the synthetic image of Figure 1. It is clear that the the miss-classification problem is avoided with the iterative method.

2.2. L 1 fidelity term. Replacing the classical L 2 fidelity term by an L 1 fidelity term is an idea that has first been introduced in signal processing [START_REF] Alliney | Digital filters as absolute norm regularizers[END_REF][START_REF] Alliney | Recursive median filters of increasing order: a variational approach[END_REF][START_REF] Alliney | A property of the minimum vectors of a regularizing functional defined by means of the absolute norm[END_REF] and later in image processing [START_REF] Nikolova | Minimizers of Cost-Functions Involving Nonsmooth Data-Fidelity Terms. Application to the Processing of Outliers[END_REF][START_REF] Nikolova | A Variational Approach to Remove Outliers and Impulse Noise[END_REF]. Using L 1 fidelity term makes the problem more robust to noise and outliers in the signal or image.

It has been remarked that for TV denoising, a L 1 fidelity is more natural [START_REF] Chan | Aspects of total variation regularized L1 function approximation[END_REF], since in this case the problem is scale-invariant. Recently, it has been shown that the Chan-Vese model benefits from the same properties when L 2 fidelity is replaced by L 1 fidelity. Indeed if (φ, c 1 , c 2 ) minimizes ( 12) [START_REF] Darbon | A Note on the Discrete Binary Mumford-Shah Model[END_REF]. This seems to be a desirable property since the scaling of the image should not affect the geometry of the segmented objects.

E L 1 CV,g (φ, c 1 , c 2 ) = Ω |g -c 1 |H(φ) dx + Ω |g -c 2 |(1 -H(φ)) dx + Ω |DH(φ)| dx, then (φ, λc 1 , λc 2 ) minimizes E L 1 CV,λg
Computing the Euler-Lagrange equation of 12 yields to the gradient descent equation ( 13)

φ t = |∇φ| div ∇φ |∇φ| + |g -c 1 | -|g -c 2 | ,
where δ(φ) has again been replaced by |∇φ|. Now c 1 and c 2 are the median values of g in {φ ≤ 0} and {φ > 0} respectively instead of the mean values. These values are to be updated at each time step. Computing the median is more demanding than the mean, but this is negligible compared to the time required to compute the update of the level set function at each iteration. How does one decide which fidelity term to chose? The L 2 -fidelity is more sensible to noise as it can be seen in Figure 4. It comes from the fact that the mean value is affected by large values whereas the median is affected by frequent values. Hence, if the object to be segmented has a very distinct color, the L 2 -fidelity term should be preferred. If the object has color close to the rest of the image or if the image is very noisy, the L 1 -fidelity will be more efficient. The CT scan of the thorax falls in the second category. We illustrate this by doing an iterative segmentation of a 2D slice using the L 1 -fidelity, see Figure 5. Note that the number of steps needed for the heart segmentation drops from 5 to 3.

Initial condition.

The Chan-Vese energy may admit some local minimums, depending on the image g and on the scale of the parameter µ. The choice of a good initial condition is then crucial in order for the curve not to get stuck into a local minimum which is not the global minimum. The convergence to the steady state may also be slowed down by a bad choice of initial condition. We tried four different strategies and compared the results obtained. They are summarized in figure 6.

The first strategy consists of taking a circle as an initial curve and its signed distance function as level set function. Many level set functions may describe the curve, but the signed distance function has the advantage of being the most regular. The second strategy is to take the union of many circles spread out over the image as an initial curve. Again the signed distance function to this curve is taken as level set function. The idea is that if the curve is spread out, it should be close to features of interest. The first two strategies are used in the seminal paper of Chan and Vese [START_REF] Tf Chan | Active contours without edges[END_REF].

In the same direction, we propose to take as initial condition a level set function that takes random values in [-1, 1] over the domain Ω. Then there are point inside and outside the curve almost everywhere. The level set function is not regular, but this doesn't seem to be a problem since the level curves gets smooth very fast as time evolves.

The last strategy comes from the fact that when µ = 0, the Chan-Vese problem is very easy to solve via a k-means procedure. This splits Ω into two regions Ω 1 and Ω 2 . One can take as an initial condition φ = χ Ω1 -χ Ω2 . The signed distance function to ∂Ω 1 ∩ ∂Ω 2 would be another option, but numerical experiments suggest that it is a better idea that all points are initially relatively close to 0.

We tried the different strategies on the well known cameraman picture. On Figures 7 and8 we compared the decay of the energy and the L 2 -convergence to the absolute minimum. The energy is approximated using a discretization of the equation 3, the term DH(φ) is simply discretized using a forward Euler scheme. If a regularized version H (φ) of H(φ) is used, the energy becomes less accurate as the function φ evolves since φ develops sharp gradients. Figure 9 shows the semi-log plot of the L 2 distance to the absolute minimizer.

On these figures, it is possible to see two different behaviors. At the beginning of the iteration process the fidelity term is dominant in equation 9. In this case, the random initial curve and the (µ = 0) solution converge very quickly towards the steady state, as seen in Figure 8. As the curves get closer to the steady state, the curvature term gain importance and then every curve approaches the steady state at a similar rate as shown on the semi-log plot in Figure 9. The same kind of phenomena happens for the energy. It is not easy to determine which initial condition is better than the other. It seems that the random curve and the (µ = 0) solution converge quickly to a state that is close to the steady state, which is certainly a nice property. However in practice we favor the use of the random curve since it has no a priori on the position of the interface and the points inside and outside the curve are evenly spread over the image. In this case, there is a lower chance to get stuck in a local minimum.

2D Application

As an illustration of the method, we have presented in Figure 5 the results of the iterative segmentation with L 1 norm on a 2D CT slice. In order to use the results of the segmentation process as a base for numerical simulations of the organ functions, it is usually required to generate a mesh of the given organ. To do so, we used the simple and very efficient 2D/3D mesh generator DistMesh [START_REF] Persson | A Simple Mesh Generator in MATLAB[END_REF]17]. Distmesh is a Matlab tool designed for meshing domains implicitly defined through a level set function. Most mesh generators start by meshing the boundary of the domain (a curve for 2D domains, a surface for 3D domains) and then march to mesh the interior. DistMesh has a different approach: it triangulates the whole domain and then moves vertices to well fit the boundary, with a control on the size and quality of the elements. It can also be modified to accommodate subdomains. For each sub-domain, we have computed the signed distance function to the boundary. The signed distance function is computed using the reinitialization equation ( 14) where S(φ) is a regularization of the sign function [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]. This signed distance function is used to project on the boundary nearby vertices. It is then sufficient to compute the signed distance function only in a neighborhood of the zero level set. With DistMesh, an element diameter has to be specified. This one may differ spatially. We remarked that nice boundaries can be obtain if the element size is decreased continuously when approaching sub-domains boundaries. In practice, the element size is usually decreased to half the size it has in the sub-domain. Figure 10 shows an example of a 2D mesh of the thorax, built from the previous segmentation in Figure 5.

φ t = S(φ)(1 -|∇φ|),

3D Application

We will apply our framework to 3D segmentation of a full CT scan. This CT scan is by courtesy of the Heart Institute of the University of Ottawa. It is of size 512×512×199 and the voxels have a 0.48mm resolution in the transverse plane (the x -y plane) and a 1.25mm. resolution in the direction of the transverse axis (z direction). This image has over 52 million voxels, it is thus a challenge to solve the segmentation problem in an efficient way.

However, since the equation is solved with an explicit time scheme, it is possible to parallelized the solver. The idea is to split the image into pieces that are to be distributed on computing nodes. A simple choice is to split the z direction into as many blocks as we want. At each iteration, the computations can be made independently on each block. Then only the information about shared boundaries of the block need to be exchanged. In the case where the image is split in only one direction, there is only the information about the top slice and the bottom slice to be exchanged, which can be done very quickly. If we apply the Chan-Vese model with L 2 fidelity, we need to compute the average of a function in a given region. This is easy to do in parallel. However, if we are to use the L 1 fidelity, we need to compute the median in parallel, which is trickier. To find the median in the serial case, the C++ template nth element is used [18], which is a linear algorithm for partial sort (it is more efficient than complete sort which is of order n log(n), where n is the number of pixels). Hence for preliminary tests, we used the L 2 fidelity. For all computations, a random initial curve has been chosen. The computations then require less time steps.

The parallel code is implemented in C++ using the openMPI library [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message Passing Interface[END_REF]. The amount of RAM memory required to solve the Chan-Vese problem is then divided into the different nodes. The speed up for this parallelization is nearly perfect: using n processors divides the time per iteration by n, see figure 11.

Most computations are made on a 16 processor SUN cluster with distributed memory. Running the code on 6 processors divides the CPU time by almost 6. As the number of processors increases, the time lost in data transmission becomes more important. Note that many computations have been made successfully on a simple dual-core laptop, cutting the CPU time in 2, although the same amount of RAM is required.

Figure 12 shows the 3D CT image as well as the result of the first segmentation step. Following this first segmentation, one side of the surface must be chosen as the new segmentation domain. The interior is the region that contains the trachea, the exterior contains the heart. To obtain the trachea we segmented the interior region with a high curvature term, since there is much noise and irregularities in the lunges. Once the segmentation process is complete, we just need to chose the connected component that corresponds to the trachea. Figure 13 shows the resulting segmentation. The trachea is then meshed with DistMesh. Figure 14 shows 2 different meshes of the trachea: a coarse one of about 40 000 tetrahedra, and a finer one of about 500 000 tetrahedra.

Now, if we want to segment the heart, we have to take the result of the first segmentation and choose the exterior of the surface. To get to the heart, four more segmentation steps will be required (five steps in total). If this would to be done with several level set functions, it would then take 5 level set functions to get this level of details, which is hard to solve on images of this size. From the second stage, we have decided to do the segmentation not on the exact image but on a blurred version of it, due to the high level of noise. To blur the image, we apply the linear heat equation on the image for a few time steps with an explicit scheme. This is also easily done in parallel. Figure 15 shows the exterior surface of the heart in the final segmentation and Figure 16 shows the segmented heart cavities. The cavities are well segmented, especially the one of left ventricle, which is the most important part in many applications. The position of the mitral valve is also precisely captured as shown in Figure 17. The general shape of the heart is well captured, ventricles and atria are extracted as well as the aorta.

There are some imperfections in the segmentation of the heart surface near the epicardium, as the surface leaks in some areas. This is due to the fact that the heart touches the liver at that place in the image and that the two organs are of same grey level. Further work has to be done to fully extract the heart.

Conclusion

In his paper, the efficiency of the iterative version of the Chan-Vese model have been shown on test cases and successfully applied to large 2D and 3D applications to the human trachea and heart. We also conclude that the L 1 fidelity is more efficient for this kind of segmentation. It is less sensitive to noise and to small regions of different color intensity. Analysis of convergence under various initial conditions showed that our new approach which consists of taking a random initial curve is very efficient. The curve quickly converges to a state close to the steady state and has no shape prior. It has been successfully apply in 2D and 3D for CT scan segmentation. We also explained how the algorithm is implemented in parallel and the computing gain obtained this way. Finally, high quality fine 2D and 3D meshes have been produced. These meshes are actually used for Finite Element and Finite Volume simulations in electro-physiology.

( 2 )

 2 Introduction of an Iterative version of the Chan-Vese algorithm to replace multi-phase segmentation. (3) Analysis of the convergence and efficiency of the Chan-Vese algorithm under different initial conditions and different fidelity terms.

1 . Background 1 . 1 .

 111 The Mumford-Shah energy and the Chan-Vese model. An image can be interpreted as a function g : Ω-→ R,

Figure 1 .

 1 Figure 1. An example of an image for which segmentation with two level set function may lead to incorrect segmentation. (a) is the synthetic image to be segmented. (b) shows a correct segmentation of the image, when the curvature term is not dominant (µ = λ = 1). (c) is the result if the curvature term is more important (µ = 10000, λ = 1): there is a miss-classification of some pixels. (d) shows a close-up on the curves where there is miss-classification: it comes from the fact that level curves are only nearly superposed in these regions.

Figure 2 .Figure 3 .

 23 Figure 2. 2D Iterative segmentation of the heart from a CT scan image. (a) The image to be segmented (b) the result of the first application of the 2-phase Chan-Vese model, the blue part is the region of interest (c) Second step: in green is the region that has been ignored in the segmentation. The red region will be chosen. (d) Red is chosen (e) blue is chosen (f) the blue region is the heart muscle, we stop.

Figure 4 .

 4 Figure 4. Segmentation of the cameraman image using L2fidelity (left) and L1-fidelity (right), with the same weights (µ = 1 and λ = 0.1)

Figure 5 .

 5 Figure 5. Iterative segmentation of the heart from a 2D slice of a CT scan using the L 1 -fidelity. 3 steps are required instead of 5 for the L 2 -fidelity.

Figure 6 .

 6 Figure 6. The different initial curves: (a) one circle, (b) 100 equidistant circles, (c) random values in [-1, 1] and (d) the solution of the problem when µ = 0 and λ = 1.

Figure 7 .

 7 Figure 7. The energy decay under various initial conditions.

Figure 8 .

 8 Figure 8. The L2-convergence under various initial conditions.

Figure 9 .

 9 Figure 9. The semi-log plot of the L2-convergence under various initial conditions.

Figure 11 .

 11 Figure 11. The speed up for the parallel algorithm: the time for doing an iteration with one processor divided by the time for doing an iteration with x processors. It is just below the perfect speedup, that is the line y = x.

Figure 12 .Figure 13 .

 1213 Figure 12. The CT scan to be segmented (a) and the result of the first segmentation step (b)

Figure 14 .

 14 Figure 14. Two meshes of the trachea generated with DistMesh: a coarse mesh of about 40 000 tetrahedra (a) and a finer of about 500 000 tetrahedra.

Figure 15 .

 15 Figure 15. 2 views of the exterior surface of the heart in the final segmentation. The general shape is well recovered.

Figure 16 .Figure 17 .

 1617 Figure 16. 2 views of the interior surfaces of the heart in the final segmentation. The pillars in the left ventricle are well segmented.

Acknowledgment

The authors would like to thank the Ottawa Heart Institute for providing them with quality CT scans of the thorax. The authors acknowledge the support of NSERC under a Discovery Grant and a CGSD (Canada Graduate Scholarship for Doctoral studies) post graduate scholarship for the first author.


J.-L Seichepine

S Albertazzi

P Nardin

Belfort/F C Coddet

D Sporer

Stockport / Gb

Topological characterization of wear tracks in thermally sprayed abradable coatings

Abradable coatings are located on the stationary parts of gas turbines, in front of blades, which cut a track in them. This has to be achieved with minimum wear of the blades, in order to control the over-tip leakage. These coatings are generally deposited by thermal spraying of composite powders comprising a metal base, a polymer filler generating porosities and a dislocator such as hBN. The very demanding properties are nowadays adjusted using rig tests, where samples are rubbed by the contact of a dummy, simulating actual working conditions in an aircraft engine. Several types of behaviour are usually described, but few numerical data are produced from these tests. Only the blade wear (or metal transfer) is generally measured. As the understanding of contact phenomena is fundamental for the development of predictive models allowing the design of more performing materials, a comprehensive characterization process of the rub path was developed. The study was based on a topological survey made by laser profilometry, giving three-dimensional maps. These maps were then processed by image analysis and several parameters were computed, like surface roughness and parameters giving information on the shape and orientation of the holes or grooves in the rubbed surface of the samples.

1

Introduction

Air seals in aircraft jet engines represent one of the most challenging and key technologies world-wide. Well functioning seal systems allow to keep rated pressure ratios, the best efficiency and an optimized specific fuel consumption. For instance, abradable coatings are used to control the over-tip leakage of air in the compressor between rotor blades and casings, labyrinth seals as well as drums and stators. Increasing performances of engines induce more severe working conditions for the coating , which have to be continuously improved. That is why it is important to gain a better understanding of the complex phenomena occurring when a blade contacts an abradable coating. Efforts are presently done in two complementary directions: modelling the thermo-mechanical behaviour of the seal and characterizing it experimentally. This paper presents a route allowing to go further in the analysis of the results given by the industrially used abradability test. It provides additional information useful for understanding and then simulating the abradability.

2

Abradable coatings

Abradable seals allow the blade to cut a track in the coating thickness but this has to be achieved with minimum wear of the blade. Furthermore, the abradable material must not erode due to the hot gas flow otherwise sealing effectiveness is lost. Many such gap-sealing products are available for use in the compressor section of gas turbines [START_REF] Schmid | An Overview of Compressor Abradables[END_REF].

Our study was focused on materials already used or being presently developed to gain the combination of properties required in the high pressure compressor. These abradable coatings are deposited by combustion or plasma spraying of composite powders, comprising a metal base, a polymer filler generating porosities and a dislocator. This combination results in deposits which are very porous, easily removed by the contact of the blade, but able to resist to the severe environmental conditions.

The reference materials for this study, already industrially used were AISi-BN (powder supplied by Sulzer Mateo under the reference SM 320) and NiCrAl-Bentonite (powder supplied by Sulzer Metco under the reference SM 314). The new material being presently developed is NiCrAl-BN, with different rates of BN and sprayed with polyester, in order to increase the porosity. The contact phenomena occurring in an engine are reproduced between a dummy blade and an abradable specimen. The rig is equipped with two independent drive systems, one for the rotor and the other for moving the seal specimen radially into the rotor. The speed of the specimen radial motion, or incursion rate, and the rotor speed are adjusted to obtain a realistic simulation of different typical working conditions. Furthermore, a heating system allows to make tests at increased temperatures [START_REF]Sulzer lnnotec Abradable Facility: Standard Abradability Test Description[END_REF] [3] [START_REF] Schmid | New High Temperature Abradables for Gas Turbines[END_REF]. The main results of these tests are measurements of the blade wear (or metal transfer from the coating to it), records of temperatures, and observations of the rubbed zones on both the blade and the coating. So several typical behaviours have been identified, like pure cutting, smearing, grooving or melting. This information is classically summarized in wear maps [1] [4]. Meanwhile, only few numerical data are usually produced from the analysis of the rubbed surfaces.

Profilometry

A comprehensive characterization process of the rub path cut by the blade in a tested coating was developed.

The study was based on a topological survey made by laser profilometry. By this technique, the distance from a reference plan to the scanned surface was measured with a resolution of 50 points/mm in both directions, giving a three-dimensional map [START_REF] Keyence | Specification of Laser Displacement Sensor LK-031[END_REF].

The general shape of the rubbed zone is cylindrical, because of the circular motion of the blade. This mean shape is removed numerically, giving a table of the distances of the surface points to the mean surface, which can be plotted as a map of grey levels, Fig. 2.

From this data, different parameters can be computed, like the surface roughness Sa.

S3 = --x L Id(x, y) 1 n-1 [m-1 ] nxm x .. o y=O
Where:

-d(x,y) is the distance between the (x,y) point and the mean plan, -n is the number of points sampled in x direction, -m is the number of points sampled in y direction.

5

Image analysis

The analysis of some wear tracks showed different topographies, with craters isotropically distributed or with grooves in the direction of the blade motion, or even with grooves perpendicular to the blade motion.

To get more information, the maps obtained by laser profilometry were processed by image analysis.

In a first step, the grey-level maps were binarized, giving black-and-white pictures where the areas above the mean surface are white and the areas below the mean surface are black. Then the shapes were simplified with the "despeckel" filter, that replaces each pixel by the mean value of its 5 neighbours [6). This filter was run several times, until the stabilization of the contours, Fig. 2.

To go further in the characterization of the peaks shape, the local dimensions of the white areas of these pictures were examined extensively. The whole pictures were scanned in both directions x and y, detecting the white sets of pixels encountered and measuring their lengths. The results are tables of lengths hx and hy, Fig. 3. ( 1 ) I ( 2 ): Segments of white area in x I y direction, having a length of hx I hy.

6

An origin al parameter

The tables of values hx and hy were used to compute a new parameter called "anisotropy coefficient", giving numerical information about the shape and orientation of the grooves or holes in the coating. The idea was to have it equal to -100 for a square picture perfectly banded in the horizontal x direction, equal to 0 for a perfectly isotropic picture and equal to +100 for a picture perfectly banded in the vertical y direction, This coefficient was calculated in several steps. First the mean values fl;;;; and hyO of the lengths of the whole sets of segments are determined. Then, 2 ranges of segments are arbitrarily defined. Those whose length is lower than the mean values are arbitrarily considered as belonging to small insignificant areas. Those whose length is higher than the mean values are considered as describing the main significant white areas or main peaks and the mean values of their length are calculated.

-11.1 is the mean value of the lengths hx higher than their mean value fl;;;; .

-11;1 is the mean value of length hy higher than their mean value n;Q.

The anisotropy coefficient c is finally defined by the following equation:

7
Analysis of a set of wear tracks 26 wear tracks belonging to 9 abradable specimens were characterized by laser profilometry. For each of them the roughness Sa and the anisotropy coefficient were computed. Table 1 summarizes the test conditions (blade tip speed, incursion ration), the measured blade wear (as a percentage of the total incursion) and the results obtained from laser profilometry and image analysis (roughness Sa and anisotropy coefficient c).

.. An original coefficient was defined to give numerical information about the anisotropy of the wear track surfaces. It appeared that this process gives a useful tool to go deeper into the understanding of the complex abradability phenomena. 9

Ill >. CJ "tJ c ftl Ill c. .. Q) .2 Q) � Q) 0 c CJ) Q) Q) c .. Q) c D. c. !!! .. .. ftl g� � E II) i:l &! Q) " §, II) "tJ 0 ftl c. ftl ::J •c: � (,)

3

 3 Abradability testThe performance of the coatings are usually evaluated by abradability rig testing. Tests are conducted in a dynamic abradability rig, Fig.1.

Fig. 1 .

 1 Fig. 1. Schematic of Standard Test Rig.

Fig. 2 .

 2 Fig. 2. Typical wear track scan. y is the direction of blade motion.-Left: grey levels.-Right: Binarized and filtered.

Fig. 3 .

 3 Fig. 3. Image analysis on a binarized wear track scan.

Fig. 4 .c = - 100 Fig. 4 .

 41004 Fig. 4.

Fig. 5 . 3 Fig. 6 . 4 Fig. 7 .

 53647 Fig. 5. Roughness versus blade wear for the 26 wear tracks analysed. On this table as usually a negative blade wear means that metal has been transferred from the coating to the blade tip. The given anisotropy coefficients are mean values of 5 values computed from square areas along each wear track. The computed surface roughnesses Sa have been compared to the blade wears. Some consistency appeared: the surface roughnesses are rather high in case of severe blade wear and rather low in case of metal transfer, Fig. 5.

Table 1 .

 1 Test conditions and results for 26 wear tracks.

		II)	i=	..5	iii	0 et:	c( CJ
	NICrAI -Bentonite	m/s um/s % 305 2,5 3 305 138 33 153 760 71	um 56 48 62	-2 24 40
			305 25	46	93	78
	AISI -BN 305 138	-14	25	65
			153 760	-20	16	64
			305 2,5	0	47	0
	z al 0 ..... ' � (,) z 0 en	1 2 3	305 138 153 760 305 2,5 305 138 153 760 305 2,5 305 138	0 18 0 -2 36 0 -2	44 60 39 49 46 34 34	-8 9 -2 4 1 2 -4
			153 760	3	41	73
			305 25	3	51	-19
		1	250 500	3	44	-15
	z al . � (,) z 0 co	2 3	153 760 305 2,5 250 500 153 760 397 10 397 100 397 100	28 3 50 41 5 -5 0	68 41 35 40 39 37 37	26 0 -1 6 1 -2 5
		4	397 397 100 10	-5 -5	21 18	0 14

0 ....

Literature

Acknowledgements

This study was funded by the European Commission under the FP5 Growth Program. It was completed within Seal-Coat on-going project, which involves the following partners: ;.. ESIL, Dublin/IE ;.. Rolls-Royce pie, Derby/GB ;.. MTU Aero Engines GmbH, Mi.inchen/DE ;.. Euromat GmbH, Heinsberg/DE ;.. RTWH -Thermal Spraying Group, Aachen/DE ;.. Neomet Limited, Manchester/GB ;.. Institute of Plasma Physics, Praha/CZ ;.. RALSA, Langreo/ES ;.. UTBM -LERMPS, Belfort/FR