Elodie Torres
email: elodie.torres@irsn.fr

Jean Desquines

Marie-Christine Baietto

Michel Coret

Franz Wehling

Martine Blat-Trieix

Antoine Ambard

Adsorption and diffusion of hydrogen in Zircaloy-4

Keywords: Zircaloy-4, Hydrogen, diffusion, precipitation, modeling

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Under normal operation in nuclear reactors, zirconium based cladding alloys are affected by water corrosion and hydrogen pick-up according to Equation 1.

Zr + 2H 2 O ZrO 2 + 2H 2 (Eq. 1)
Depending on the corrosion temperature, the dissolved hydrogen can precipitate when exceeding the material solubility limit. In pile, the migration of hydrogen induced by thermal gradient in the cladding radial direction promotes accumulation of hydrides at the cladding outer diameter. After several years irradiation, a hydride rich rim is formed as illustrated in Figure 1. Hydrogen has to be considered in many safety issues. Indeed, it changes many properties of zirconium alloys such as their mechanical or metallurgical behaviors. Hydrogen charged specimens are formed in laboratories to investigate part of the behavior of the irradiated material. Two main hydrogen charging techniques are used: cathodic charging at room temperature and gaseous charging at a temperature close to irradiation conditions. Considering this second charging technique, rather few studies address the hydrogen adsorption and diffusion in the metal. In this work we intend to provide a modeling of the gaseous charging of zirconium specimen as well as an experimental validation.

Modeling

Markowitz experiments evidenced that hydrogen diffuses towards the colder regions [START_REF] Markowitz | Hydrogen redistribution in thin plates of zirconium under large thermal gradients[END_REF] showing that a driving force for solid solution diffusion is temperature gradient. The diffusion is also driven by concentration gradient and to a lower extent by stress gradient as shown by Marino and Kammenzind [START_REF] Marino | A numerical calculation of the redistribution of an interstitial solute in a thermal gradient[END_REF][START_REF] Marino | Hydiz a 2 dimensional computer program for migration of interstitial solutes of finite solubility in a thermal gradient[END_REF][START_REF] Bf Kammenzind | Hydrogen pickup and redistribution in alpha-annealed Zircaloy-4[END_REF]. In the following, the modeling of thermal diffusion only neglects the stress gradient influence as suggested by Kammenzind [START_REF] Bf Kammenzind | Hydrogen pickup and redistribution in alpha-annealed Zircaloy-4[END_REF].

The total hydrogen content (C T ) is thus considered to be the sum of two components, the precipitated hydrogen content (C p ) and the dissolved hydrogen content in the -phase (C d ) according to Eq. 2, with the concentrations expressed in weight particles per million (wppm).

(Eq. 2)

During thermal diffusion, under time-increasing total hydrogen content, the dissolved hydrogen concentration is limited, assuming equilibrium conditions, by the hydrogen Terminal Solid Solubility for Precipitation (TSSP), C TSSP . During this precipitation phase, the dissolved hydrogen content in the -phase remains close, but not systematically equal, to the C TSSP limit [START_REF] Westerman | Charging Zircaloy-2 with hydrogen beyond the solubility limit[END_REF]. As a simple summary, during precipitation (Eq. 3), the solid-solution content is associated to equation 4.

(Eq. 3) (Eq. 4)

C supercharging corresponds to a few wppm dissolved hydrogen excess compared to the TSSP. This excess was evidenced by Marino [START_REF] Marino | Hydrogen supercharging in Zircaloy[END_REF] and is usually called "hydrogen supercharging".

The hydride precipitation kinetics is rather slow [8-9] and is assumed to be proportional to the level of supercharging as proposed by Marino and Kammenzind [5;[START_REF] Marino | Hydrogen supercharging in Zircaloy[END_REF]. This hydride precipitation kinetics is described in equation 5, where <X> is the positive value of X.

〈 〉 (Eq. 5)

The proportionality coefficient is given by equation 6 [START_REF] Westerman | Charging Zircaloy-2 with hydrogen beyond the solubility limit[END_REF],

(Eq. 6)
where  is a parameter related to the hydride precipitation kinetics in s

-1/2
, and T is the temperature in Kelvin.

During heating, the local hydrogen content is governed by hydride dissolution. Under such conditions, the dissolved hydrogen content is limited by the Terminal Solid Solubility for Dissolution (TSSD). Additionally, the dissolution kinetics is a fast process. During this phase, the dissolved hydrogen content is simply assessed by a maximum value being equal to C TSSD . The dissolution phase is defined by equation 7 and 8.

and

(Eq. 7)

(Eq. 8)

At a given temperature, the C TSSD value is significantly higher than C TSSP [START_REF] Une | Dissolution and precipitation behavior of hydrides in Zircaloy-2 and high Fe Zircaloy[END_REF][START_REF] Wh Erickson | The influence of alloying elements on the terminal solubility of hydrogen in alpha-zirconium[END_REF][START_REF] Slattery | The terminal solubility oh hydrogen in Zirconium alloys between 30 and 400°C[END_REF][START_REF] Jovanovic | Thermal diffusion of hydrogen and hydride precipitation in Zr-Nb pressure tube alloys[END_REF]. The dissolution and precipitation solubility limits are given by equations 9 and 10. TSSD is frequently considered as the solubility limit at equilibrium [START_REF] Do Northwood | Hydrides and delayed hydrogen cracking in Zirconium and its alloys[END_REF]. However, there is some controversy on this point.

( ) ( ) (Eq. 9)

( ) ( )
(Eq. 10)

The hydrogen diffusion within hydrides is slow and neglected [15]. The local hydride fraction is considered as negligible in order to simplify the diffusion equation. The hydrogen flux, without intense hydrostatic stress, derives from equation 11,

⃗ ( ⃗ ⃗⃗ ⃗ ⃗⃗ ) (Eq. 11)
where J is the hydrogen flux (mm/s), Q* is the heat of transport, equal to 24000 J/mol [16-17], T is the absolute temperature (Kelvin), R is the gas constant (J/mol/K) and D is the diffusion coefficient of hydrogen in the -phase (mm 2 /s). D can be described from an Arrhenius law according to equation 12 [18]. in mm 2 /s (Eq. 12)

In the presence of precipitated hydrides, the -phase diffusion is described by equations 13 and 14.

⃗⃗⃗⃗ 〈 〉 (Eq. 13)

⃗ ⃗⃗ (Eq. 14)

The literature was intensively analyzed in order to identify reference experiments for code validation [16-17;19-21] considering thermal diffusion only.

The modeling has been extended to account for solid-gas exchange. In the literature, it has been essentially addressed at elevated temperature. The hydrogen concentration in metal versus gaseous pressure at equilibrium was experimentally evidenced to follow a Sievert law [22-23] described in equation 15,

( ) (Eq. 15)
where k is a kinetic constant, p is the gaseous pressure and Q is the activation energy in the α-phase region.

The modeling relies on an explicit 1D integration scheme of diffusion equations in axisymmetric or cartesians coordinates. The model provides hydrogen radial concentration profile after any simulated exposure at elevated temperature. To model the solid/gas exchanges two possible assumptions are considered: a constant interface concentration (Eq 15 like conditions) or a constant absorbed flux depending on the environment gas pressure p as described by equation 16,

(Eq. 16)
where k is a kinetic constant.

Let's consider a thin plate of thickness, e, to clarify the meaning of this boundary condition. Solid/gas exchange at both faces of this plate following Equation 16 type boundary conditions and isothermal hydrogen charging are considered.

Solving diffusion equations provides a very simplified assessment of the diffusion behavior in such component. The average hydrogen concentration in the sample (C T ) is given by time integration of equation 17.

(Eq. 17)

In other words, at a given gaseous pressure, there is a constant increase rate of average hydrogen content in the sample. This kind of boundary condition can lead to unlimited hydrogen adsorption, whereas equation 15 type leads to equilibrium between hydrogen content and gaseous pressure.

To qualify the modeling, several gaseous charging experiments have been performed on SRA Zircaloy-4 cladding samples. In order to provide fully coupled experiments involving hydrogen motion in metal and surface exchange, complementary tests have been carried out using plain rod RXA samples.

Materials and methods

1-Materials

The specimens investigated in this study have been cut from low tin 17x17 SRA Zircaloy-4 fuel cladding tubes and RXA rods provided by AREVA-NP. The tube specimens are 70 mm long, 9.5 mm outer diameter and 570 µm thick respectively. Oversized grains located at the outer and/or inner diameter have been observed on as-received rod samples. Consequently, rods have been spark machined from 9,8 to 7 mm of outer diameter to remove the large grains zone. Chemical composition of these samples is detailed in Table 1

2-Hydrogen Gaseous charging

The hydrogen gaseous charging is performed at the École Centrale de Paris. Each Zircaloy-4 sample is inserted in a quartz tube under vacuum. The sample is then heated by a furnace at 420°C. Once the target temperature is reached, pure hydrogen is injected into the quartz tube. Indeed, the injected pressure is linked to the target hydrogen content. The charging temperature is chosen sufficiently low to avoid cladding recristallisation and sufficiently high for allowing charging high hydrogen concentrations in a reasonable time. The charging can be performed in several hydrogen injection steps depending on the desired final concentration. The sample is then cooled down without applied stresses at 20°C/h to obtain large and representative of as-irradiated circumferential hydrides (Figure 2). The hydrogen pressure decrease in the charging device provides a direct measurement of the average hydrogen content absorbed within the sample. The maximal allowable hydrogen concentration increase in the specimen is directly related to the initial pressure in the furnace.

3-Experiment

To describe only surface/gas exchange, gaseous charging experiments have been performed on SRA Zircaloy-4 cladding tubes. Due to the thin geometry of tubes, the hydrogen is close to uniformly distributed and the surface exposed to gaseous hydrogen is large, permitting an accurate measurement of adsorption kinetics.

RXA Zircaloy-4 rods have been also used to perform radial diffusion experiments. After gaseous charging, the cylindrical sample is first machined to obtain a cone and the outer volume is then truncated into many 1 mm long crowns with inner diameter depending on its axial location after sample machining (Figure 3). The hydrogen is measured on each crown sample which is allowed to determine radial hydrogen concentration assuming axial homogeneity. These complementary tests have been performed to provide fully coupled experiments involving hydrogen diffusion and surface/gas exchange.

Figure 3 : Conical machining and cutting after gaseous charging of the samples 4-Hydrogen content measurement

The hydrogen content of the samples is measured after gaseous charging using a Brücker ONH mat 286 by melting the sample at 1500°C in the presence of a carrier gas (argon). After calibration, the hydrogen content is obtained by comparing the conductivity of the mixture with the one of pure argon. For tube samples, four axial locations are systematically investigated. A minimum of four analyses for each area are performed to characterize the homogeneity of hydrogen content at each location. This method provides a measurement of the average hydrogen content and the radial hydrogen profile cannot be directly measured.

Results and discussion

1-Hydrogen diffusion in a thermal gradient

The experimental data of Sawatzky and Markowitz [16;[START_REF] Sawatzky | Hydrogen in Zircaloy 2; its distribution and heat of transport[END_REF] have first been considered for the validation of the diffusion modeling. The Sawatzky experiments have been performed on Zircaloy-2 cylindrical samples. The first test (called A) consisted in keeping during 34 days a linearly varying temperature of 130°C and 477°C at each ends of a 24 mm long sample containing initially 130 wppm hydrogen homogeneously distributed. In the second test (called B), a linearly varying temperature between 157°C and 454°C along the sample is kept during 41 days using a 24 mm long the sample containing initially 64 wppm hydrogen homogeneously distributed. The Markowitz experiments have also been performed on Zircaloy-2 cylindrical samples. In the first test, a thermal gradient of 73°C/cm during 30 days has been applied on the sample containing initially 300 wppm hydrogen homogeneously distributed. In the second test, a thermal gradient of 75°C/cm has been kept during 41.5 days on a sample containing initially 3360 wppm hydrogen uniformly distributed. Figure 4 compares the two test results of Markowitz with calculations and Figure 5 presents the Sawatzky results compared to simulations. A good agreement is obtained between experiments and simulation. There are, however, some discrepancies, mainly identified at the end of the samples. The two-phase, α/α+δhydride transition is acceptably simulated. However, in Markowitz first test (Figure 4A), this experimental transition is not as clear as in the simulation. Indeed, for high concentrations, hydride volume fraction should be taken into account. Actually, hydrogen diffusion within hydride is about three times slower than in the metal [15].

2-Surface/gas exchange

The present modeling has been compared to the experimental data of many samples charged by gaseous method. The experimental average hydrogen content within the sample is deduced from the hydrogen gaseous pressure decrease in the charging device. An example, where the charging has been performed in several hydrogen injection steps, is presented in Figure 6. First we can check that pressure decrease is really due to hydrogen absorption. The maximal allowable hydrogen concentration increase in the specimen is directly related to the initial pressure in the furnace by Equation 18,

(Eq. 18)
where is the initial hydrogen pressure of the current injection step (bar), m sample is the weight of the sample (g), is the maximal hydrogen concentration increase at each injection step i (wppm), V 0 is the hydrogen gas volume in the device, R is the gas constant (J/mol/K), T is the hydrogen charging temperature (Kelvin) and M H is the hydrogen molecular weight (1 g/mol).

Let

be the hydrogen content variation at each step (wppm) defined by equation 19,

(Eq. 19)
where is the initial hydrogen pressure and is the decrease in hydrogen pressure during absorption. The total average hydrogen content in the sample (

) is thus expressed by equation 20. In this equation, C 0 is the hydrogen content in the as received material.

∑

(Eq. 20)

The average hydrogen concentration of the studied tubes was approximately 100, 200, 400 and 500 wppm. The average hydrogen content calculated in the sample using eq. 15 for gas absorption can be compared to the measurements. The main results are summarized in Table 2. The good consistency between experimentally measured hydrogen content and the calculated one from pressure decrease shows that pressure can be used to measure the sample hydrogen content.

Table 2: Hydrogen contents after gaseous charging versus code predictions

Whereas hydrogen adsorption is clearly governed by a Sievert law at elevated temperature [22;23], the hydrogen charging experiments can lead to the absorption of the whole gaseous content of the furnace as illustrated in Figure 6. This trend is not compatible with Sievert like boundary conditions (Eq. 15) but with constant adsorption flux (Eq. 16 like condition). Analysis of the hydrogen absorption kinetics allows us to determine the pressure exponent of Eq. 16.

Hydrogen absorption kinetics can be fitted to a "Langmuir" type law [24;25]. The rate of average hydrogen content in the sample follows equation 21.

(Eq. 21)

where A is a kinetics constant and µ is a term proportional to .

The adsorption rate is proportional to , this means that the rate limiting phenomenon is the H 2 molecule adsorption at the sample surface rather than H 2 dissociation [25].

Considering that in our experiments , the combination of Eq. 17 and Eq. 21 results in Eq. 22, (Eq. 22)

where A*µ is in , k in , e in mm and D in mm 2 /s.

In the case of tube-type specimens and low hydrogen pressure, the solid/gas exchange is therefore described by Eq. 23,

(Eq. 23)
where is in and p in bar.

The agreement between experiment and theory is rather good as illustrated in Figure 7 in the whole range of hydrogen contents and charging conditions. The weight gain (g/cm 2 .s) is also used to characterize the hydrogen gaseous charging kinetics. The modeling provides assistance to laboratory gaseous charging experiments. The expected hydrogen profile and sample average contents can be anticipated from calculation. The modeling leads to a homogeneous distribution of the hydrogen which is consistent with the metallographic examinations (Figure 8).

3-Combined hydrogen diffusion and Surface/gas exchange

Complementary tests on RXA Zircaloy-4 rods have been performed to provide fully coupled experiments involving hydrogen distribution and surface/gas exchange. The radial profile of hydrogen content derived from measurements and the one derived from simulations are illustrated in Figure 9. An inverse technique based on destructive characterization of the sample has been used to deduce the hydrogen profile. This technique consists in postulating a hydrogen distribution and adjusting its shape to find average hydrogen content as close as possible to measurements on crown pellets machined on the rod sample (Figure 3). The hydrogen concentration varies between 60 wppm at the center of the rod up to 100 wppm at the outer diameter. For the simulations, the solid/gas exchange coefficient and/or the diffusion coefficient was adjusted to reproduce as closely as possible the experiment. The diffusion coefficient is potentially influenced by material texture and the charging kinetics appears to be faster than for tubes and some differences grain size and microstructure could also affect both, diffusion coefficient and solid/gas exchange. Kearns mentions that the diffusion coefficient in the c-axis direction of hcp Zircaloy-4 is about two times larger than in the a-direction. This would imply a lower diffusion coefficient to accurately simulate rods [18]. The measured hydrogen profile is acceptably reproduced by simulation. A sensitivity study to diffusion coefficient and adoption kinetics was done by adjusting the D 0 and k parameters in Eq. 24 and Eq. 25. It was not possible to find parameters leading to the measured profile in the entire range of radii. However two extreme sets of values are identified, the first one fitting the distribution at the center of the cylindrical sample and the second describing accurately the distribution at the outer diameter of the sample as illustrated in Figure 9. The results are different but remain consistent with the expected values obtained for tubes (Eq.12 and Eq. 23).

Conclusion

In this study, a quantitative modeling of hydrogen gaseous charging of Zircaloy-4 samples is proposed taking into account both the solid/gas exchange and the thermal diffusion of hydrogen in the metal.

Some gaseous charging experiments at 420°C of thin Zircaloy-4 fuel cladding tubes have been interpreted assuming a hydrogen flux at the solid/gas interface proportional to the hydrogen pressure. The modeling is consistent with measurements and can now be used to monitor hydrogen charging experiments. The expected hydrogen profile and sample average contents can be anticipated from calculations.

Experiments were designed on rod samples in order to obtain a radial hydrogen profile on RXA Zr-4 rods. The solid/gas exchange kinetics and diffusion coefficient had to be adjusted to reproduce the measured hydrogen profile in the sample. This is probably due to a strong sensitivity of the diffusion coefficients and solid/gas exchange to material microstructure.

Figure 1 :

 1 Figure 1: Cladding microstructure in PWR normal operating condition [1]

Figure 2 :

 2 Figure 2: Example of hydrides distribution in the cladding after gaseous charging at 420°C

Figure 4 :

 4 Figure 4: Hydrogen distribution in Markowitz samples [17] (A) 41.5 days anneal (B) 30 days anneal

Figure 5 :

 5 Figure 5: Hydrogen distribution in Sawatzky samples [16] (A) 34 days anneal (B) 41 days anneal

Figure 6 :

 6 Figure 6 : Evolution of hydrogen gaseous pressure versus time at 420°C of a sample charged at 500 wppm in 3 steps with an injection hydrogen pressure around 18 mbar

Figure 7 :

 7 Figure 7 : Hydrogen charging experiments of 8 different tubes at 420°C: experiments versus modeling (A) hydrogen content evolution in the sample (B) hydrogen charging kinetics

Figure 8 :

 8 Figure 8: Distribution of hydrides in a Zircaloy-4 tubes after gaseous charging with 100 wppm of hydrogen

Figure 9 :

 9 Figure 9: Hydrogen distribution within the sample: measurements versus modeling

Table 1 : Chemical composition of the Zy-4 samples used for the hydrogen gaseous charging

 1 .

	Sn (wt%)	Fe (wt%)	Cr (wt%)	O (wt%)	H (wppm)
	1.3	0.21	0.10	0.13	10


Michal Swierczewski

Plinio Maroni

Alexis Chenneviere

Mohammad Dadras

Lay-Theng Lee

Thomas Bürgi

Deposition of extended ordered ultrathin films of Au_38(SC_2H_4Ph)_24 nanocluster using Langmuir-Blodgett technique


Amal Zeaiter
email: zeaiteramal@gmail.com

Etienne Videcoq

Matthieu Fénot

Determination of Electric Motor Losses and Critical Temperatures through an Inverse Approach

Keywords: Electrical machines, Future time steps, Losses, Lumped Parameter Model, Thermal behavior, Regularization

In this study, a practical numerical method is proposed to estimate losses of high specific power density electric motors, using few simulated temperature data. In such electric motors, these losses generate high heat fluxes inside the motor components that can be critically sensitive to temperature. Electromagnetic and mechanical friction phenomena are behind the occurring of these thermal dissipations. For both phenomena, losses could be difficult to compute with electrical or mechanical approaches. However, thermal management of electric motors requires a precise knowledge of those losses, in particular for high-performance motors such as those considered in future hybrid planes. To determine electric motor losses in a Permanent Magnet Synchronous Motor (PMSM) in real time, an inverse method using a Lumped Parameter Thermal Model (LPTM) is elaborated. In the first step, the dynamic profile of losses is determined through the inverse method, based on temperature data at easy-access points of the motor. In a second step, the identified losses are used to find temperatures at critical non-accessible hot spot points of the motor through forward LPTM. The method is applied for three useful cases, from the simplest case scenario, where only one type of losses has to be identified, to the most complicated case where all losses are simultaneously estimated. A global strategy for the choice of the number of future time steps used for regularization of the ill-posed problem is also proposed. Results show that this method enables adequate real-time supervision of the critical motor temperatures, mainly rotor and winding core.

I. INTRODUCTION

Nowadays, with the rising interest of electrification in the transportation sector, the study and conception of very high-specific power density electric motors is essential. Indeed, weight and volume are known to be classical constraints for electrified vehicles. Moreover, the recent use of electric motors for the propulsion of the new generation of ground and air vehicles requires the development of series of more efficient and higher performance electric motors [START_REF] Bilgin | Making the Case for Electrified Transportation[END_REF]. This driving trend already exists in the automotive field for electric and hybrid-electric vehicles [START_REF] De Santiago | Electrical Motor Drivelines in Commercial All-Electric Vehicles: A Review[END_REF]. Recently, propulsion electrification is emerging into air vehicles design development. Numerous projects [START_REF] Siemens | Electric propulsion components with high power densities for aviation[END_REF]- [START_REF]Airbus[END_REF] are being conducted aiming to develop electrically propelled aircrafts and are targeting higher values of specific power density in electric motors. In general, these projects focus on hybrid aircrafts in which electric motors, coupled to fans, would propel the plane, and the turbo-reactor will be mainly used as a turbine for electricity generation. The increase of specific power density of electric motors could be limited by electrical and thermal constraints. The multiple electrical limits and the approaches allowing an increase of the motors specific power density are not discussed here, since the related issues are not the main purpose of the paper and are associated to loadability concepts. The interested reader may refer to [START_REF] Lefevre | Performance assessment tool based on loadability concepts[END_REF] for further comprehension of these issues. Thermal limits in electric motors, majorly related to maximum temperature allowed in electric motors windings depending on the insulation thermal class, are at least of the same order of challenge -if not more -than the electrical limits [START_REF] Popescu | Modern Heat Extraction Systems for Power Traction Machines-A Review[END_REF][START_REF] Nategh | A Review on Different Aspects of Traction Motor Design for Railway Applications[END_REF]. Excessive temperature increase may irreversibly demagnetize magnets [START_REF] De | Multiphysics Design Optimization of a Permanent Magnet Synchronous Generator[END_REF] or could lead to shortening motor lifetime due to overheated electric insulation material [START_REF] Farsane | Experimental study of the cooling of a closed type electric motor[END_REF]. The temperature increase in electric machines results from electromagnetic and mechanical losses, generating heat in motor components. Generally, these losses are empirically estimated. Besides being hardly evaluated with analytical and some computational methods, it is believed that even with a motor prototype, the process of determination of the resulting losses and distinguishing their types or origins is quite complicated. Generally, such information can be barely obtained with long iterative methods.

Nomenclature

Our study focuses on a Permanent Magnet Synchronous Machine (PMSM). In such machines, the temperatures are very sensitive to losses, specifically at hot spot zones. A possible solution to identify electric motor losses is to measure temperature inside the motor and then apply an inverse method. Such methods require a thermal model with quick solving process. Thus, a thermal nodal network could be a good choice as it is commonly used to carry out numerical simulations in the design process of high-performance electric machines.

Being mathematically ill-posed, the solution of the inverse problem may be unstable and not unique. Hence, several regularization techniques have been developed to ensure stable solutions [START_REF] Woodbury | Inverse Engineering Handbook[END_REF]. Examples of inverse problems for electric machines can be found in the literature with different objectives: geometric, thermal or magnetic parameter estimation, loss mapping, source strength determination from temperature measurements. In [START_REF] Khelifa | An inverse problem methodology for design and optimization of an interior permanent magnetic BLDC motor[END_REF], the authors propose an inverse problem in order to identify a set of 9 design parameters (geometry, magnetic parameters) for a permanent magnet brushless motor. The results show that the density of the magnetic flux in the air gap is increased, as well as the motor efficiency. Moreover, an inverse method for the identification of heat convection coefficients on coil end-windings from experimental data is presented by [START_REF] Davin | Motor Cooling Modeling: An Inverse Method for the Identification of Convection Coefficients[END_REF]. An interesting algorithm is developed and described for this nonlinear inverse heat conduction problem using a lumped parameter model. In [START_REF] Vansompel | An Inverse Thermal Modeling Approach for Thermal Parameter and Loss Identification in an Axial Flux Permanent Magnet Machine[END_REF], a low-order model is built in a first stage to estimate a set of thermal parameters. In a second stage, this model is used to identify thermal losses in the stator of an axial flux PM machine. In [START_REF] Nair | Sensitivity Analysis of Inverse Thermal Modelling to Determine Power Losses in Electrical Machines[END_REF], the power loss density is estimated using a 2D finite-element thermal model associated to a conjugate gradient method from noisy temperature measurements. A prediction of the strength of the heat source field in an induction motor is also carried out in [START_REF] Nair | Inverse Thermal Modeling to Determine Power Losses in Induction Motor[END_REF] using a steady state thermal model. Another example of reconstruction of thermal field sources can be found in [START_REF] Marinova | Inverse source problem for thermal fields[END_REF]. The aim is to estimate a 3D source distribution using infrared data. The least square method is used in order to solve the inverse problem. Concerning now the on-line identification of electromagnetic losses, a recursive parameter estimation method is implemented in [START_REF] Hey | On-line monitoring of electromagnetic losses in an electric motor indirectly through temperature measurement[END_REF]. Furthermore, a 3D inverse heat conduction problem is presented by [START_REF] Huang | A three-dimensional inverse problem in estimating the internal heat flux of housing for high speed motors[END_REF] to estimate the timedependent internal heat flux generated in rotor and stator for high speed electric motor.

This article deals with the temperature prediction in the electric machines at low accessibility locations through the identification of the heat sources and using an inverse technique. The originality of our numerical study lies in the fact that a low order lumped parameter thermal model that requires a very low computation time is used first to identify the unknown heat sources through an inverse procedure, and secondly to compute the temperature field in the electric machine at each time step in real time. The technique used to solve the mathematically ill-posed problem relies on Beck's contribution [START_REF] Scott | Analysis of order of the sequential regularization solutions of inverse heat conduction problems[END_REF], [START_REF] Osman | Numerical solution of the general two-dimensional inverse heat conduction problem[END_REF], which takes into account the lagging and damping effects due to the diffusion process by using a sequential method associated to a function specification procedure. This technique acts as a regularization procedure, as it will be shown in this study.

After presenting the electric machine and associated heat transfer modeling, the inverse method applied to identify the heat sources dissipated inside the motor is detailed. Finally, results of parametric studies and corresponding analysis are presented.

II. THERMAL MODEL

A. DESCRIPTION OF THE ELECTRIC MOTOR

The electric motor illustrated in Fig. 1 is a Surface Mounted-Permanent Magnets Synchronous Motor (SM-PMSM) used for hybrid electric aircraft propulsion with specifications found in [START_REF] Touhami | Original Optimization Procedures of Halbach Permanent Magnet Segmented Array[END_REF][START_REF] Touhami | Electro-thermal Models and Design Approach for High Specific Power Electric Motor for Hybrid Aircraft[END_REF]. It is a high efficiency electric motor (more than 96%) and maximum losses resulting in thermal sources are around 16 kW. These losses are divided into winding Joule losses, Iron losses and mechanical losses.

Geometrically, the motor consists in two parts: the rotor (rotating in the center) and the stator around it.

As shown in Fig. 1, the rotor consists of a shaft (brown) supported by two bearings (grey), surrounded by the rotor lamination (yellow) and the permanent magnets (orange). In the studied motor, the considered heat sources are due to mechanical losses and electric losses. Mechanical losses are mainly due to friction inside the bearings (PBearings), and to aerodynamic heating resulting from very high velocity gradient in the air gap between the rotor and the stator (PAirgap) and near the rotor end-caps (lateral surfaces of the rotor) termed PRotor.

The stator consists of the stator laminations (purple) and the winding (green) whose biggest part is embedded in the stator slots and the remaining part are the end-windings at the extremities of the stator, which are the complex parts in the cavities from geometrical and thermal viewpoints. Windings are heterogeneous elements made up of two types of materials: copper conductors and electric (and thermal) insulation material. The iron losses (PIron) and Joule losses (PJoule) are thermal sources distributed in the stator laminations and winding respectively rising from the resistivity of copper in winding conductors. Iron losses are produced by magnetic effects (hysteresis and Foucault's current) and the Joule losses in windings are due to Joule effect.

It is assumed in this study that the thermo-physical properties are temperature independent. The motor cooling system in this case is a combination of a stator water-jacket with water-glycol flowing around the stator, and a water-glycol axial flow at the center of the shaft for rotor cooling. Both types of liquid cooling are being used in electric motors for vehicular applications (examples can be found in [START_REF] Gai | Cooling of Automotive Traction Motors: Schemes, Examples, and Computation Methods[END_REF]). After flowing through the motor, the water-glycol is cooled with a simple heat exchanger: a channel in contact with ambient air through a flat plate. Moreover, high conductivity resin potting of end-windings is used to increase the heat transfer in those critical regions (due to their high Joule losses). This technique is investigated in [START_REF] Tighe | Assessment of cooling methods for increased power density in electrical machines[END_REF][START_REF] Nategh | A Study on Thermal Effects of Different Potting Strategies in Traction Motors[END_REF]. These techniques have dynamic impact on motor temperatures and thus influence the extraction of the generated heat resulting from the losses. They must be taken into account in any proposed model of the motors.

B. HEAT TRANSFER EQUATIONS AND MODELING

Temperature distribution within the motor is determined from solving the heat equation associated with initial and boundary conditions. The equation governing heat diffusion and convection in the domain Ω is written as:

(1)

Where T(M,t) is the temperature at a point M in the domain and is the velocity field, ρ the density, Cp the specific heat and k the thermal conductivity. Pj is the time varying strength of the j th volumetric heat source and Volj the corresponding volume. The characteristic function χj(M) is equal to 1 if M is located in source j, else χj(M) =0.

Heat exchanges between the motor and its environment are modeled through a global heat transfer coefficient h. Boundary conditions hence are written as follows:

(2) Moreover, the initial condition is:

(3)
If the velocity field is not known in a particular region of domain Ω, the transport term in (1) may be modeled by a convective heat transfer as in [START_REF] De Santiago | Electrical Motor Drivelines in Commercial All-Electric Vehicles: A Review[END_REF]. Then a convective coefficient h has to be evaluated from Nusselt correlations. Besides, the fluid temperature Text is required.

C. STATE SPACE REPRESENTATION

As analytical solutions for this 2D axisymmetric transient problem are not available, a Lumped Parameter Thermal Model (LPTM) is used to simulate motor thermal behavior according to (1)-( 3). It consists in dividing the studied system into N discrete volumes, which are connected together by thermal conductances. For brevity reasons, details of the network are not given in this paper but can be found in [START_REF] Zeaiter | Thermal Sensitivity Analysis of a High Power Density Electric Motor for Aeronautical Application[END_REF]. It leads to a system of ordinary differential equations with a low number of parameters. In our study, N=58. This set of equations can be written as:

(4)

Where T(t) (dimension N) is the vector of temperatures, function of time t, at the N discretization nodes, its derivative with respect to time t. Matrix A (dimension N,N) is the state matrix which connects temperatures at discretization nodes and contains diffusion and transport terms, as well as terms related to convective boundary conditions. The vector Bc is associated with these latter conditions and links corresponding nodes to environment nodes. Matrix Bp (dimension N,np) links discretization nodes to the internal heat sources gathered in vector function

P(t)=[PJoule PIron PRotor PAirgap PBearings] T (dimension np).

In what follows, let us separate vector P into two sets: unknown heat sources included in U (dimension nU) and known heat sources included in K (dimension (np-nU). The command matrix Bp is then also split into two matrices Bp=[BK BU] Defining the vector V=BcText+BKK, the state space representation takes the following form:

(5a) (5b)
Where an observation matrix C (dimension nq,N) allows to select nq temperatures in the whole temperature field T(t) and to store them in vector function Y(t).

In this purely numerical study, a maximum of 9 temperature locations will be used as outputs. These temperatures are depicted in Fig. 1, and are located in the stator lamination (Ts1, Ts2 and Ts3), in the end-winding (Tw) at the surface of the end-winding (Tsw), on the surface of the stator in the airgap (Ta), in the cavity between the rotor endcap and the frame (Tc), in the bearings (Tb), and in the magnets (Tm). These locations have been chosen given the heat sources but also considering plausible experimental measurement system positions. Temperature in the static part of the motor could be obtained using thermocouples (Ts1, Ts2, Ts3, Tsw, Ta, Tb and Tc). The rotor flange temperature could be obtained by infrared measurement. Tm and Tw are the critical temperatures difficult to obtain experimentally due to their locations, i.e. Tw inside the end-windings is confined and thought to be the hottest point and Tm requires a slip ring due to the rotating rotor.

III. INVERSE PROBLEM SOLUTION

The form of heat transfer model, developed in subsection II.C, is convenient to be used now in the inverse procedure. The inverse problem aims at finding an estimation of the input vector U from measured temperatures Y*. Note that, in this study, real temperature measurements are not available. The temperatures included in vector Y* derive from the direct problem with known heat sources. In order to simulate measurement errors, each temperature is altered with an additive Gaussian error, whose standard deviation is given by σ.

Subsequently, in order to have simplified and clear annotations, we will write Uk for U(tk), Tk for T(tk), Yk for Y(tk).

A sequential estimation method is adopted. Having the vector of noisy temperatures , and an estimated of , an estimation of is then computed. Using an Euler implicit scheme and Δt being the time step, (5a) becomes: [START_REF]Airbus[END_REF] Then: [START_REF] Lefevre | Performance assessment tool based on loadability concepts[END_REF] Hence, according to (5b): [START_REF] Popescu | Modern Heat Extraction Systems for Power Traction Machines-A Review[END_REF] Due to the lagging and damping effects of heat diffusion and convection, a variation of heat source strength does not immediately affect the sensors. In order to take into account these effects, Future Time Steps (FTS) [START_REF] Beck | Inverse heat conduction: ill-posed problems[END_REF]- [START_REF] Blanc | A guide for the use of the function specification method for 2D inverse heat conduction problems[END_REF] are used. It means that output data at further times tk+2, tk+3, … are used to correctly estimate . If nf is the number of FTS, then for 1≤f ≤nf, (8) is written for k+1+f instead of k+1. A temporary approximation of Uk+1+f is therefore needed to look for . In this study, a constant value is chosen:

for (9)
Hence, a global matrix formulation is obtained as in ( 10):

With:

(11) (12)
And:

(13)

The size of vector Y* is (nf+1)×nq and matrix S is of size ((nf+1)×nq,nU). The objective is to identify the pseudosolution of the inverse problem, such that Y*-Y≈0, where Y is the temperature vector computed by LPTM. As matrix S is not square due to the addition of future time steps, the least square method is used to solve [START_REF] De | Multiphysics Design Optimization of a Permanent Magnet Synchronous Generator[END_REF] and leads to the sequential solution: [START_REF] Davin | Motor Cooling Modeling: An Inverse Method for the Identification of Convection Coefficients[END_REF] In order to evaluate the accuracy of the estimations for the different studied cases, a mean quadratic discrepancy between exact heat sources (vector U exact ) and estimated ones (vector ) is defined as follows: [START_REF] Vansompel | An Inverse Thermal Modeling Approach for Thermal Parameter and Loss Identification in an Axial Flux Permanent Magnet Machine[END_REF] Where nt is the number of time steps in the inverse problem. It is worth to remind that, in a practical application, this quantity is not reachable.

In a similar way, for each output Yi, a mean quadratic error between measured temperatures (or artificially noisy temperatures) and those computed with the estimated set of heat sources strengths, is defined as follows: [START_REF] Nair | Sensitivity Analysis of Inverse Thermal Modelling to Determine Power Losses in Electrical Machines[END_REF] When the unknown heat sources are identified, they will be used as inputs in the direct problem in order to predict temperature in inaccessible parts of the motor, for example critical temperatures Tcrit in windings or magnets. The complete procedure is illustrated in the block diagram of Fig. 2. At each time-step, the model is solved in inverse problem to get motor losses. The inputs for inverse LPTM are the noisy temperature at accessible points (sensor temperature measurement on a real bench).

Then, forward LPTM is used to determine indicated critical temperatures.

IV. RESULTS AND ANALYSIS

A scenario of motor operation is chosen corresponding to a profile mission in an aircraft flight application. The corresponding profiles of heat sources generated in the motor are depicted in Fig. 3 during 3600 s. The different phases correspond to: taxi-out phase (0 to 450s), take off and climb (450 to 900s), cruise flight (900 to 2250s), descent and landing (2250s to 2700s) and taxi-in (2700s to 3600s). Rotor, Airgap, and Bearings losses remain relatively low and quite constant all along the mission whereas Joule and Iron losses change significantly. Those losses are particularly high during the take-off phase as they are linked to propulsion power. Airgap and Bearings losses depend principally on the rotor rotation speed, which varies slightly between take-off and cruise, and is very low during taxi phases.

Temperatures are then computed using LPTM and a time step equal to 1 s. Among these temperatures, a set of observables is chosen depending on the different case studies (see subsection II.C). A noise of standard deviation equal to σ=0.1°C is added to these observables so that they are used as temperature measurements in various inverse problems cases.

A. CASE 1: IDENTIFICATION OF IRON LOSSES

In this first case, among the five heat sources, four are supposed to be known (PJoule PRotor PAirgap PBearings) and only iron losses (PIron(t)) have to be identified. Fig. 4a shows the evolution of temperature Ts1 with respect to time. From this temperature, the inverse problem is solved in order to identify P(t)=PIron(t) knowing K(t)=[PJoule PRotor PAirgap PBearings] T . Hence, for this case, nq=nU=1.

The procedure is sequential, as mentioned before, expressed in ( 14), and depends on the number of future time steps. Results of identification are presented in Fig. 4b and Fig. 4c for respectively nf=2 and nf=10. It is remarkable that, without regularization (nf=0), the oscillations of PIron (t) are huge. Hence, we have chosen to draw PIron (t) for nf=2 in Fig. 4b.

It clearly appears that the number of future time steps filters the oscillations of PIron, what is confirmed by the quadratic criterion σU in Table 1. The main criterion, and the only one to be computed in a real case, is σY. Its evolution with respect to nf is quite interesting. For nf=0, even if the problem is ill-posed, large oscillations in allow to follow exactly the evolution of the noisy temperature. Hence, as indicated in Table 1, σY≈0°C. When nf increases, regularization tends to damp the time evolution of , leading to a bias in the temperature computation. We choose then the number nf for which σY≥σ to "get out" of the noise level. The identification of PIron(t) for nf=10 is depicted in Fig. 4c and shows a very good agreement with the exact evolution of these losses drawn in Fig. 3. Once PIron is identified, the aim is to predict temperatures Tw and Tm in inaccessible parts of the electric motor. They have been computed for the different values of nf and are shown in Fig. 4d for nf=0. Note that, for nf=0, temperature prediction occurs without any delay whereas it is executed with a 10s delay with nf=10 since the time step Δt is equal to 1s. For any value of nf, the temperature evolution in critical zones is the same, due to the fact that oscillations in are damped by the thermal path and inertia of the system. Same as in [START_REF] Vansompel | An Inverse Thermal Modeling Approach for Thermal Parameter and Loss Identification in an Axial Flux Permanent Magnet Machine[END_REF], critical temperatures Tw and Tm are compared to the exact temperatures computed with the exact inputs. For nf=0, a quadratic error σTcrit=0.003°C is obtained, which gives a very accurate prediction. Now, a new test is carried out using a second measurement point Ts2 located in the stator. Graphical results are summarized in Fig. 5.

It appears obviously that the oscillations of PIron(t) decreased compared to Fig. 4, which is confirmed by Table 2 where σU goes from 229 W to 200 W for nf=10. It is interesting to observe that in this case (nq=2, nU=1) it is not possible to obtain σY≈0°C for each of the two outputs from only one heat source PIron(t). The temperature prediction in windings and magnets is also very accurate with an error σTcrit=0.002°C for nf=0.

B. CASE 2: IDENTIFICATION OF JOULE AND IRON LOSSES

For the second case, Joule losses (PJoule) are also unknown. The evolutions of PJoule(t) and PIron(t) have to be identified. Firstly, outputs Ts1 and Tsw located in the stator and at the winding surface are used. It can be remarked in Fig. 6 that both heat sources are well identified when increasing the number of future time steps. In this case, nf=12 gives σY≥σ for each output, as shown in Table 3.

The condition number of the matrix S T S to invert is also indicated. It is noteworthy that this last one decreases when the number of future time steps increases, which confirms the regularization effect of the function specification, knowing that in case 1, S T S is a scalar, hence its condition number is equal to 1. Compared to case 1, the winding temperature is quite sensitive to the oscillations in Joule losses identification, hence the error increases: σTcrit=1.161°C for nf=0. In the other hand, by adding two other outputs Ts2 and Ts3, located in the stator, better inversion results are again obtained, as depicted in Fig. 7 and Table 4, since σU goes from 200 W to 180 W. The addition of outputs acts as a filter in the inversion procedure and the condition number of matrix S T S decreases also. The temperature prediction in windings and the magnets remains relatively accurate with an error σTcrit=1.158°C for nf=0.

C. CASE 3: IDENTIFICATION OF ALL LOSSES

In this third case, the five heat sources altogether are considered unknown and have to be identified simultaneously from five outputs Ts1, Tc, Tsw, Ta and Tb, whose temperature evolutions are given in Fig. 8a. These nodes are located in the stator, in the cavity, on the surface of the windings, in the air-gap and on the bearings respectively. It can be remarked that the temperature magnitude is very different according to the outputs locations.

As for previous cases, the quadratic errors and the condition number are given with respect to the number of future time steps. In this case, the condition number of the matrix S T S is very high, hence, the noise in temperature measurements may be greatly amplified by the inverse problem. As in previous cases, nf=0 corresponds to a square case (nq=nU=5), thus, even if the problem is ill-posed, large oscillations in allow to follow exactly the noisy temperature evolutions. Hence, as indicated in Table 5, σY≈0°C. An increase of nf leads to σY≥σ for all outputs for nf=12. Fig. 8c shows large oscillations in PRotor due to the fact that the outputs are located far from this heat source.

Concerning now the temperature prediction in windings and magnets (Fig. 8d), the error is completely acceptable with σTcrit=1.178°C for nf=0.

V. CONCLUSION

The study presented in this paper deals with the determination of electric motors' losses or heat sources, and prediction of critical temperature at non-accessible parts. A low-order lumped parameter thermal model is used to solve both inverse and direct problems. The method is sequential and includes regularization with the function specification technique. Different configurations are tested according to the number of heat sources to be identified. For each case, the influence of the required number of outputs to be used to solve the inverse problem and the number of future time steps are analyzed. The technique gives accurate prediction of the temperature in unreachable critical parts of the motor and is computationally efficient: only 10 -3 s at each time step is needed for identification of the heat sources and prediction of the motor temperatures.

FIGURE 1 .

 1 FIGURE 1. SM-PMSM schematic section with specific temperatures' locations.

FIGURE 2 .

 2 FIGURE 2. Block diagram: use of LPTM for both inverse and direct problems.

FIGURE 3 .

 3 FIGURE 3. PMSM heat sources evolution as function of time PMSM.

FIGURE 4 .

 4 FIGURE 4. Iron losses identification from node Ts1 with nf=2 and nf=10 and computation of critical temperatures Tw and Tm.

FIGURE 5 .

 5 FIGURE 5. Iron losses identification from nodes Ts1 and Ts2 with nf=2 and nf=10 and computation of critical temperatures Tw and Tm.

FIGURE 6 .

 6 FIGURE 6. Joule and iron losses identification and computation of critical temperatures from nodes Ts1 and Tsw with nf=2 and nf=12.

FIGURE 7 .

 7 FIGURE 7. Joule and iron losses identification and computation of critical temperatures from nodes Ts1, Tsw, Ts2 and Ts3 with nf=2 and nf=10.

FIGURE 8 .

 8 FIGURE 8. Identification of all the heat sources and computation of critical temperatures from nodes Ts1, Tc, Tsw, Ta and Tb with nf=4 and nf=12.

TABLE 2 . Inversion results for iron losses identification from Ts1 and Ts2; bold lines correspond to results shown in Fig. 5.

 2

	nf
	0	3713.2	0.0690	0.0700
	2	673.8	0.0921	0.0929
	4	336.0	0.0957	0.0960
	6	237.5	0.0976	0.0975
	8	205.1	0.0992	0.0987
	10	200.1	0.1011	0.1001
	12	205.0	0.1034	0.1020

TABLE 3 . Inversion results for Joule and iron losses identification from Ts1 and Tsw; bold lines correspond to results shown in Fig. 6.

 3

	nf				CN(S T S)
	0	4125.2	3e-14	6e-15	3.74
	2	741.2	0.0859	0.0869	3.64
	4	359.2	0.0926	0.0931	3.54
	6	246.4	0.0949	0.0971	3.45
	8	209.3	0.0970	0.1007	3.36

TABLE 4 . Inversion results for Joule and iron losses identification from Ts1, Tsw, Ts2 and Ts3; bold lines correspond to results shown in Fig. 7.

 4

	nf					CN(S T S)
	0	2849.0 0.0803	3e-5	0.0832 0.0816	1.26
	2	518.4	0.0941 0.0867 0.0974 0.0952	1.23
	4	262.6	0.0964 0.0937 0.0994 0.0974	1.21
	6	194.5	0.0976 0.0975 0.1004 0.0987	1.18
	8	179.2	0.0990 0.1012 0.1012 0.0999	1.16
	10	180.3	0.1007 0.1055 0.1025 0.1016	1.13
	12	188.3	0.1028 0.1109 0.1040 0.1038	1.11

TABLE 5 . Inversion results for all losses identification from Ts1, Tc, Tsw, Ta and Tb; bold lines correspond to results shown in Fig. 8.

 5

	nf

Declarations

This project has received funding from the European Union's Horizon 2020 (Cleansky 2 JTI) research and innovation program, 2014-2024 under grant agreement No 715483.