Stéphane Mouty
email: stephane.mouty-ext@converteam.com

Abdollah Mirzaian

Frederic Gustin

Alain Berthon

Daniel Depernet

Christophe Espanet

Permanent Magnet Generator Design Solutions for Wind Turbines

Introduction

Converteam is a worldwide specialist in power conversion and in offshore solutions. Wind energy needs power converters for energy conversion and transport. Converteam has realized wind turbine generators and their converters for several years. The company has decided to propose only direct drive solutions for future orders and is developing new concepts of permanent magnet (PM) generators in collaboration with the Energy department of FEMTO-ST Institute.

In order to be competitive in the market and to respect the environmental constraints, the international standards have to be respected and economic solutions are researched. Thus, the aim is to obtain a machine which has a low cost and a high efficiency. The full power converter, which is necessary to connect the machine to the grid, needs to be optimized too.

Firstly, a description of the common supply chains met in high power wind turbine applications is given. Secondly, advantages and drawbacks of the direct drive permanent magnet generators will be detailed in comparison with the other topologies. Then, the study will focus on the permanent magnet generators and improvements to increase their interest. A comparison between various permanent magnet generator designs focusing on active parts weight and efficiency is proposed.

High Power Wind Turbine Conversion Chains

In order to optimise energy production, the use of variable speed generators is required for large wind turbine, i.e. power production above 1 MW.

Four main topologies are met in the market for high power wind turbine systems.

The nominal generator speed is different according to the considered topologies.

1. Doubly Fed Induction Generator (DFIG)

The most popular system is induction machine (see Fig. 1). There are two reasons to explain this interest, firstly induction machine cost is lower than synchronous machine and secondly power electronics converter doesn't need to convert all the power but only around thirty per cent. These converters are used to modify the rotor properties and allow variable speed behaviour around the rated point.

To keep a good efficiency, an induction machine will be used at high speed (>1500 rpm).

To obtain this speed, the only way is to use a gearbox, the turbine speed is around 15 rpm for high power wind application, so that some gearbox's stage are necessary. Nevertheless with a three stages gearbox the power is limited at 5 MW.

The main suppliers are: Vestas, Gamesa, GE Wind, Nordex, Alstom wind …

Synchronous Generator with Wound Rotor (SGWR)

With the experience in high power machines (oil, hydro or nuclear) where synchronous generator with wound rotor (or direct current machine in old installations) is used, it is not surprising to find this type of technology in wind turbine generator. This topology is met for a direct drive application (without gearbox). In this case, the working speed is the same that the turbine speed; thus at a given power the torque is more important than one obtained with the DFIG. Machines will have a big diameter and a short length; due to their form these machines may be called torque or ring generator.

This topology leads to expensive active parts and lots of auxiliaries are necessary to create the rotate field as we can see Fig. 2.

Two topologies of excitation exist: firstly with a brush to commutator contact, secondly with a rotating transformer and a diode bridge. This second solution is preferred because there are none mechanical contacts between rotated and fixed parts. To supply the rotating transformer, converters are necessary but their power is lower than the generator one. To work with variable speed, a full power converter is necessary to convert the electrical power provided by the stator of the machine. With such a power chain, the machine is able to produce energy for almost all wind speeds.

The main supplier is Enercon.

3. Permanent Magnet Generator

In order to increase power density and to simplify field creation of the direct drive generator, a solution is to use permanent magnets instead of wound rotor. In this case the power chain is the same, as shown in Fig. 3.

The main suppliers are GE and Siemens.

4. Hybrid Concept

A hybrid concept, with permanent magnet machine and a gearbox having a low transformation ratio, leads to rotational speeds between 100 and 300 rpm and enables the reduction of the generator's size.

The aim is to have a smaller generator than the direct drive options and a better efficiency than DFIG. The supply converter is designed for the full power (see Fig. 4).

The main supplier is Multibrid, electrical active parts are realized by Converteam.

Direct Drive Permanent Magnet Generator (DDPMG)

Some studies can be found in the literature concerning comparison between the different wind turbine topologies. For the four previous cases we can recall the works of H. Polinder. In [START_REF] Bang | Promising direct-drive generator system for large wind turbines[END_REF], a comparison of cost and losses are given. The results are summarized in Fig. 5.

1. Magnets rather than Wound Rotor

The fact that no supply is necessary to create the rotate field is the main advantage of the permanent magnets. Most of the time, permanent magnets synchronous machines are called "brushless" to highlight this aspect.

With magnets the rotor weight and rotor losses are reduced.

In industry, most of permanent magnet machines have their magnets mounted at the surface of the rotor because there is an experience background, a mounting process which is under control and seems to be the easiest solution.

But, when magnets are mounted at the surface of the rotor, defluxing workings are more difficult to obtain because poles are not salient. These modes can be interesting when over speeds or over loads happen in wind turbine applications.

2. Direct Drive rather than Gearbox

Even if gearbox is not the main cause of downtime, risk of failures is still important [START_REF] Puigcorde | Wind turbine gearbox reliability, 2010, renewable energy world[END_REF]. The strongest constraint is the necessary time for the replacement of the gearbox when it is out of order. And this time may be long, especially in offshore applications.

Losses for a permanent magnet generator are less important than for an induction one and machine is able to work on a larger wind speed range. Moreover, choosing permanent magnet generator, instead of DFIG, enables the reduction of the entire nacelle weight [START_REF] Fairley | Wind turbines shed their gears both siemens and GE bet on direct-drive generators[END_REF], limits the number of elements in the conversion chain and so reduces the risk of failures.

Consequence

Focusing on the technology, solution with permanent magnets seems to be the best, because it leads to lower losses, lower total weight, lower risk of failures, more produced energy.

Therefore, direct drive permanent magnet generator is the topology which has been selected by the company to be suggested to the customers.

The main permanent magnets drawback is their cost, which is linked to the increase of the demand of rare earth magnets. Some risks must be taken in consideration too: demagnetisation, short circuit torque, bonding with other magnetic parts during the manufacturing, the use of non-magnetic tools.

Improvements

In spite of its technology benefits, we see that the cost of permanent magnet direct drive active parts is high. Thus, even if rules are not the same for active parts than for global turbine (i.e. to minimize the active parts cost doesn't necessary lead to minimize the wind turbine cost), we are going to propose some ways to optimise the active parts weight of the generator.

For the active parts design of a wind turbine generator, the main points which have to be taken in consideration are: power density optimization and losses reduction. These two aspects do not match because the weight increases when the losses decrease, see Fig. 6, then a compromise must be found.

Some solutions can increase the power density of the surface mounted permanent magnet machine and/or limit the losses:

• To reduce the copper with a concentrated winding and a fractional number of slot per pole and per phase [START_REF] El-Refaie | Fractionnal-slot concentratedwinding synchronous permanent magnet machines: opportunities and challenges[END_REF]; • To modify the magnet's form and adapt the currents waveforms [START_REF] Jahns | Pulsating torque minimization techniques for permanent magnet AC motor drives-A review[END_REF]; • To have an axial flux machine topologies instead of a radial [START_REF] Balkan Sirnsir | A comparison of torque capability of axial flux and radial flux type of brushless DC (BLDC) drives for wide speed range applications[END_REF]; • To use an outer rotor topology instead of an inner rotor one, if the diameter is constrained; • To increase the phases number. As some cases are ever been described in the literature we do not focus on these solutions. In the rest of the paper only influence of the phase number and the location of the rotor compare to the stator are presented. Nevertheless, in order to clarify the understanding it is necessary to describe the choice of the air gap flux density and currents waveforms.

Air gap flux density

In permanent magnet machines, two shapes can be considered for the air gap flux density waveform (Fig. 7) and the back electromotive forces (back-EMF) have similar shapes. For each of these cases, in order to obtain a constant torque, the current waveform in the slots will be different: the current will be sinusoidal when the flux density is sinusoidal whereas the current will be trapezoidal (nearly rectangular) when the flux density is trapezoidal. As the waveforms are similar to the classical alternative machines, when the waves are sinusoidal, the machine is called BLAC (brushless AC).

In the others case, as the principle is closed to the brushed DC machine, these machines are called BLDC (brushless DC).

Two control strategies can be used with rectangular currents to keep the torque constant [START_REF] Qiang | Comparison of brushless DC motor drives with 180/120-degree inverter systems[END_REF]. In the first case, each phase is open circuited during a part of the period. With 3 phases, a phase must be supplied with a constant current during 120 electrical degrees, the switches command signals are adapted, as shown in Fig. 8.

In the second case, the phases are always connected so that they are fed with a constant current during 180 degrees as shown in Fig. 8.

In order to have an air gap flux density as close as a square wave we need to have magnet width close to the pole pitch and a number of slots per pole and per phase equals to one. In this case, the air-gap flux density shape is given Fig. 9. To obtain a constant torque, the EMF waveforms should be constant during 180 electrical degrees, which is only a theoretical possibility. Indeed, due to flux leakages between two magnets, the flux density (and as a consequence the back EMFs) cannot be constant under the pole transition. Because of this, performances, in case of a 180-Degree, are not so interesting and it should be corrected.

Influence of Phase Number

To make a comparison between various machines for equivalent losses, the back-EMF peak values must be the same to have the same iron losses, and the current RMS values must be the same to have same copper losses. A comparison between the cases previously described is presented in Table 1. Real gain is obtained when the difference between back-EMF peak and RMS value is taken into account, this the correction mentioned in the previous section.

A comparison at equivalent losses can be interpreted as a comparison at equivalent active parts weight. Indeed, iron and copper volume are preserved and variations of magnet volume don't have a big influence on the total weight. In detail, to keep the same EMF peak value, in all cases, some geometric parameters can be preserved, as magnet thickness, air gap length and rotor diameter but the magnet width can vary. For the currents, with the same RMS value, keep the current density leads to a constant copper volume.

As a conclusion, BLDC machines have a better power density than the BLAC.

A drawback of BLDC motor is the difficulty to obtain a constant torque. With a 120-degrees converter it is necessary that the currents are perfectly rectangular and in phase with the back EMFs which must be constant on 120 electric degrees. With 180-degrees this becomes nearly unfeasible, as explained below.

The back EMFs are calculated by using finite elements simulations (FLUX2D) and they are plotted in Fig 10 for the example of a three-phased machine.. Torque is determined considering the currents as perfectly rectangular on 180 electric degrees; the waveforms are given in Fig. 11 for 3 phases.

We study now what it happens when the number of phases (i.e. the number of slots per pole in this case) increases. The winding scheme will have to be adapted according to the phase number (see Fig. 12). For a phase number of 7 the back-EMF is given in Fig. 13 and the torque in Fig. 14. It appears that the modification of the phase number has an influence on the torque ripple which does not seem to be negligible. On the other hand, when the number of phases rises, its influence becomes small as shown in Fig. 15.

The torque ripple is characterized using [START_REF] Bang | Promising direct-drive generator system for large wind turbines[END_REF].

MEAN MIN MAX Γ Γ - Γ = ∆Γ (1)
Moreover the phase number can be limited by the teeth size, which can become too small and unfeasible when the pole number is high.

Rotor Location

Once again in the aim of increasing the torque density, two structures of surface mounted permanent magnet machine can be studied. These structures are given in Fig. 16. The difference between these two structures concerns the rotating part (rotor) which is placed with regard to the fixed part (stator). In the first case corresponding to the outer rotor topology, the stator is placed in the center of the rotor. For the inner rotor topology, the rotor is in the center. A comparison between these solutions is made in the next paragraph.

Comparison of permanent magnet generators designs

Two machines respecting the specifications given in Table 2 are designed. The specifications correspond to a wind turbine application as one described in [START_REF] Polinder | Comparison of direct-drive and geared generator concepts for wind turbines[END_REF]. Previous results show that with a 180-degree converter and a trapezoidal flux density, the machine will have the best power density, so that this topology will be used for the designs.

To obtain the slightest heavy machine a minimization of the weight is made. To achieve this optimization, an analytical model is necessary to estimate the back-EMFs for this type of machine. This model is described below.

Analytical prediction of back-EMF

Radial component of air gap flux density can be determined with (2) and using (3)-( 5) for inner rotor machines. These equations can also be adapted to analyze flux density of a machine with outer rotor as in [START_REF] Espanet | Modélisation et conception optimale de moteurs sans balais à structure inverse application au moteur-roue[END_REF]. )

⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ = ∑ - - - p n r b r a p n r B odd n p n p n (2) with: ( ) c R R p n R R p n R R R R a p n p n p n p n . 2 . 1 ) . 1 ( 1 . 2 1 . . 2 2 1 . 3 2 . 3 2                 - - +         + ⋅         = + (3) ( ) c R R p n R R p n R b p n p n p n . 2 . 1 ) . 1 ( 1 . 2 1 . . 2 2 1 1 . 2                 - - +         + ⋅ = + + (4) ( ) ( ) rc a p n B K n p n p n R R p c       -                 - = 2 . sin . 1 . . 1 1 2 2 . . 2 3 1 π π (5) cos( ) ( ) , ( 1 . 1 . θ θ
where R1 is the radius at the bottom of permanent magnets, R2 is the radius at the top of permanent magnets, R3 is the stator interior radius, Ka is the factor between magnet arc and pole shoe and Brc is the corrected remanent flux density expressed in ( 6)

r R R rc B B µ µ 2 1+ = (6)
Br is the remanent flux density of the permanent magnet and μR is the relative permeability of permanent magnet.

It is then possible to predict the flux density in a slotless machine, as in Fig. 17. But to take into account the slot effect it is necessary to correct the flux density with the relative air gap permeance, as in [START_REF] Qiang | Comparison of brushless DC motor drives with 180/120-degree inverter systems[END_REF] ) , ( ) , ( ) , (

θ λ θ θ r r B r B G = (7)
where the relative permeance λ is given by ( 8) and ( 9

+       + + - = R M h g w S (9)
Expression of υ is given [START_REF] Zhu | Instantaneous magnetic field distribution in brushless permanent magnet DC motors, Part III: effect of stator slotting[END_REF], wS represent the slot opening width, NS is the slot number, g is the air gap length and hM is the magnet thickness.

The flux density waveforms taking into account the slot effects are given in Fig. 18.

The back-EMF rms value is obtained with ( 10)

Ω = G RMS G RMS LR B E ( 10
)
where L is the length of active parts, RG is the air gap radius and Ω is the angular speed in rad/s. Table 3 gives a comparison between the flux densities obtained with 2D finite-element simulations and the analytical model. Even if an error is made on the RMS flux density value when the slots are taking into account, this method gives a good estimation of the back-EMF.

Designs comparison

The stator external diameter is limited by the size of the impregnating tank, the characteristics of the two designed machines are presented in Table 4.

The torque ripple is always calculated assuming that the currents are perfectly rectangular.

Considering the design constraints of these machines, it appears that the active parts total weight decreases when the rotor is placed outer of the stator. This can be explained by the fact that the air gap diameter can be bigger when the rotor is out. Therefore only this case will be considered.

As it was foreseeable, amplitude of torque ripple is high; a design with 7 phases is made in order to limit the torque ripple. The main parameters of this design are given in Table 5. By keeping the same efficiency, increasing the phase number leads to a decrease of the torque ripple; it also appears that the weight decreases. Therefore this increase seems to correspond to an increase of the machine power density too.

All designs made are quite under the specified weight, which can be found in [START_REF] Polinder | Comparison of direct-drive and geared generator concepts for wind turbines[END_REF], for the same efficiency. In the paper mentioned, the generator is a conventional BLAC 3 phases machine. Difference on the weights agrees with the power gain given in Table 1, nevertheless this gap must be weighted by the fact that the diameter, the airgap are different.

Comments

When the number of phases increases, using a 180degrees or a (q-1)/q*180-degrees converter leads to fulfill the machine performances. But without a 180degrees converter, torque ripple amplitude becomes greater. In fact, without other improvement, a 120degrees converter permits to limit the torque ripple with 3 phases but if the number of phases is greater, a 180 degrees converter is preferable as shown in Fig. 19.

Conclusion

For a high power wind turbine application, direct drive permanent magnet generator seems to be the best technology: lower weight, lower losses and less failure risks. Nevertheless, it has some drawbacks, generator cost and size. Few improvements can be done concerning the size. In fact, less the speed is, higher the diameter should be; thus some restrictions can appear as the transport capabilities, but solutions are possible to increase the power density and limit the weight, and so the cost, of the active parts.

In the aim of increasing the machine power density, some solutions are proposed, and two points have been investigated. The first is the phase number; indeed in some cases the choice of the phase number can bring some advantages. But, this modification leads to increase the difficulty of the machine manufacturing and control. To achieve the control, models with fictitious machine could be use as proposed in [START_REF] Scuiller | Développement d'outils de conception de machine polyphasées à aimants utilisant l'approche multimachine[END_REF].

The second point deals with the rotor position respect to the stator. When the overall volume is limited, the use of an outer rotor is a solution to reduce the weight of the machine's active parts. But, the impact on the non active parts has to be studied carefully. Figure 19: Evolution of the ripple torque with the phase number with a 180-degree converter (dotted line) and a (q-1)/q*180-degree converter.

Other solutions can be found in the literature, with multi generators [START_REF]The Liberty 2.5 MW Wind Turbine[END_REF] or where magnets are not mounted at rotor surface [START_REF] Zhang | Design and comparison of a novel stator interior permanent magnet generator for direct-drive wind turbines[END_REF]. Because of the manufacturing constraints, which are given previously, these solutions are not studied by the company.

On the other hand, for the future wind turbine generators, Converteam develops a High Temperature Superconductors (HTS) machine [START_REF] Lewis | A Direct Drive Wind Turbine HTS Generator[END_REF].

Figure 1 :

 1 Figure 1: Nacelle scheme of a DFIG and description of the power chain.

Figure 2 :

 2 Figure 2: Power chain of the DDWR and picture of a machine.

Figure 3 :

 3 Figure 3: Nacelle scheme of DDPMG and description of it power conversion chain.

Figure 4 :

 4 Figure 4: Picture of the multibrid generator and scheme of the power chain.

Figure 5 :

 5 Comparison of wind turbine topologies in terms of cost (a) and energy (b)

Figure 6 :

 6 Figure 6: Evolution of the weight in function of efficiency.

Figure 8 :

 8 Machine converter (a) scheme, (b) 120-degree control and (c) 180-degree control.

Figure 9 :

 9 Figure 9: Air gap flux density when the magnet width is equal to pole shoe and its rms value (dotted lines).

Figure 12 :Figure 14 :

 1214 Figure 10: Back-EMFs in the 3 phases machine

Figure 15 :

 15 Figure 15: Evolution of the torque ripple with the phase number

Figure 16 :

 16 Figure 16: Scheme on a pole pair of outer rotor topology (left) and inner rotor topology (right)

Figure 17 :

 17 Figure 17: Air gap flux density without slot effects:Comparison between Finite-element (cross) and analytical

Figure 18 :

 18 Figure 18: Air gap flux density: Comparison between Finite-element (cross) and analytical when the phase number is 3 (a), and 7(b)

Table 2 :

 2 Specification

	Power	3 MW
	Speed	15 rpm
	Pole pair	80
	Outer stator diameter	<5 m
	Efficiency	>94.5 %
	Air gap	>5 mm
	Active material weight	<24 T

Table 3 :

 3 Comparison between FE and analytical calculations

	Method	FE	Analytical	Error (%)
	BRMS (T)	0.772	0.758	1.8
		With 3 phases
	BG RMS (T)	0.693	0.638	8
	ERMS (V)	1563	1522	2.6
		With 7 phases
	BG RMS (T)	0.731	0.69	5.6
	ERMS (V)	1647	1646	0.1

Table 4 :

 4 Comparison between the two designs

	Parameter	Inner rotor	Outer rotor
	Inner diameter (m)	4.72	4.78
	Outer diameter (m)	4.98	5.06
	Air gap diameter (m)	4.78	5
	Air gap length (mm)	7	7
	Magnet weight (T)	1.3	1.3
	Iron weight (T)	11.6	11
	Copper weight (T)	5.4	5.2
	Total active parts weight	18.3	17.5
	(T)
	Efficiency (%)	94.53	94.72
	Torque ripple (%)	37	38.9

Table 5 :

 5 Design with 7 phases and an outer rotor

	Phase number	7
	Inner diameter (m)	4.8
	Outer diameter (m)	5.07
	Air gap diameter (m)	5
	Air gap (mm)	7
	Magnet weight (T)	1.3
	Iron weight (T)	12.1
	Copper weight (T)	3.8
	Total active parts weight (T)	17.2
	Efficiency (%)	94.72
	Torque ripple (%)	4.5


R Conseil
email: robin.conseil@hotmail.fr

D Ramel

J Lebedinsky

Y Bailly

P Hervé

DETERMINATION OF HOT GAS TEMPERATURE PROFILES BY IR SPECTROSCOPY

Keywords: Main subjects: radiative transfer, inverse methods Fluid: exhaust gas Visualization method(s): absorption/emission recording, spectroscopy Other keywords: infrared spectroscopy

Determination of exhaust gas temperature and concentration is a critical point for the study and the optimization of engines combustion. A non-invasive method has been developed based on the analysis of the CO2 infrared spectrum, which is particularly sensitive to the temperature between 4.16 µm and 4.20 µm [1]. Inversion of the radiative transfer equation ensures the retrieval of temperature and concentration profiles, from transmission or radiative intensity measurement and complementary information about the profiles shape. The inverse problem is solved by minimizing the quadratic difference between the measured spectrum and the calculated one from the spectroscopic data listed in . This method has been validated on a flowing CO 2 /N 2 gas mixture heated up to 300°C. Identified temperature has shown good agreement with the one recorded by thermocouple: difference was less than 2.5%. This optical technique, without medium disturbance, may be used in various applications involving combustion, such as the study of Diesel exhaust pipe.

INTRODUCTION.

Temperature is a critical information on thermodynamic systems. In the specific case of combustion, the energetic or running optimization of systems is strongly linked to this parameter since their efficiency is directly proportional to the temperature. Environmental standards regularly impose stronger reduction of pollutants or greenhouse gas emission to the automotive industry. An additional solution to all the pollutant filtering or transforming systems such as particles filters or catalytic converter is to improve engine efficiency. In order to make such improvements, the solutions must be experimented while designers have to minimize trials and prototypes. Currently, some numerical tools are developed and their validation must be achieved by comparisons with experimental values. In this paper, a new method to determine temperature and gas concentration profiles in thermodynamic systems is proposed. This original non-invasive method has been developed as part of a research project on the engines air-loop for automotive industry: SIMBA. The resulting numerical model concerns the whole gas circulation in the engine from the inlet to the exhaust lines. Inlet pipes, turbocharger, coolers, cylinders, exhaust and EGR lines have been taken into account for the Diesel engine that has been studied. The presented method is based on the analysis of the CO 2 infrared spectrum, which is particularly sensitive to the temperature between 4.16 µm and 4.20 µm [1]. The inversion of the radiative transfer equation from transmission or radiative intensity measurements ensures the retrieval of temperature and concentration profiles. Finally, the method is applied to the exhaust of a DV6-TED4 PSA engine.

RADIATIVE PROPERTIES OF GAS AND PRINCIPLE OF THE METHOD.

The properties of the infrared radiations emitted by a hot gas are entirely different from those emitted by a hot solid surface. According to the Planck's law, the radiative intensity emitted by a blackbody is a function of temperature and wavelength (here the wavenumber ):

, . where T is the temperature of the body (K), the wavenumber (cm -1 ), 1.191 (W.cm 2 .sr -1 ) and 1.438 (K.cm). Usually, the wavelength λ of the radiation is given in µm but for IR spectroscopy, the conventional unit of the wavenumber is cm -1 . The wavenumber is given by the following formula: ! " 10000 $!μ " & .

Blackbody's law is represented on the following figure at 300K, 500K, 700K and 900K in the spectral range 0 cm -1 to 6000 cm -1 .

Fig. 1 : Energy emitted by a black body for various temperatures

On the contrary to a hot solid medium which emits radiations on a continuous spectrum, a hot gas emits radiations only on discrete wavelengths (and conversely, a cold gas absorbs radiations only on discrete wavelengths and not on a continuous spectrum). For a gas, the radiative emission or absorption is directly linked to the molecular energy level which can be divided in four types: kinetic energy, electronic energy, vibrational energy and rotational energy. Besides, these discrete wavelengths are gathered in wavelength stripes, whose width increases as the temperature, the concentration and the pressure of the gas rise. The infrared lines of carbon dioxide spectrum present an important dependence to temperature around 2400 cm -1 . Therefore, this molecule is the main gas of this study. When torque and rpm increase in an engine, the temperature and concentration of carbon dioxide increase in the exhaust gas. A Fourier Transform Infrared spectrometer (FTIR) can give really accurate information about this kind of lines and the following figure illustrates this point with carbon dioxide absorption lines measured on the exhaust gazes of a DV6 Diesel engine: As expected, temperature and concentration of CO 2 increase into the exhaust gases. The shapes of the spectra are modified and the absorption lines are getting larger. The shape and intensity of gas absorption or emission lines are deduced from the gas temperature and concentration data along the optical axis "x". At the wavenumber " " the transmission is given by: ' (

) * +, , .,+
where,is the monochromatic absorption coefficient of the gas. For a gas volume, this expression can be integrated from x =0 to x= d with "d" the thickness of the system as shown on the following figure:

Then, the radiation emitted by this specific volume is given by the transfer equation:

where, is the monochromatic transmission Since and finally:

With this last equation, the radiative heat transfer equation and nonisothermal gas volume can be divided in small elements. For each of them, the properties are homogeneous and the temperature constant. Thus and containing carbon dioxide, is here considered Here, k = 0 indicates the origin of the radiation in an emissivity "ε(σ)" and its temperature "Tp". The energy emitted by this element is defined as discretization of the temperature and concentration profiles for element, the temperature and concentration of gas are constant values attributed at the center of the element and calculated by: and 15 th International Symposium on Flow Visualization June 25-28 Then, the radiation emitted by this specific volume is given by the following form of the is the monochromatic transmission of the gas from "x" to "d" given by: the radiation can be defined by: With this last equation, the radiative heat transfer equation may be solved numerically assuming that the heterogeneous gas volume can be divided in small elements. For each of them, the properties are the temperature constant. Thus, a system divided into "N" elements of constant thickness "e = d/N" is here considered as shown below:

eterogeneous system divided into N homogeneous and isotherm elements.

Here, k = 0 indicates the origin of the radiation in the system. It may be null or bounded by a hot plate characterized by )" and its temperature "Tp". The energy emitted by this element is defined as of the temperature and concentration profiles for this system is necessary. It is element, the temperature and concentration of gas are constant values attributed at the center of the element and International Symposium on Flow Visualization 28, 2012, Minsk, Belarus following form of the integrated radiative heat assuming that the heterogeneous gas volume can be divided in small elements. For each of them, the properties are considered constant thickness "e = d/N" eterogeneous system divided into N homogeneous and isotherm elements.

the system. It may be null or bounded by a hot plate characterized by )" and its temperature "Tp". The energy emitted by this element is defined as . Of course, the assumed here that for each element, the temperature and concentration of gas are constant values attributed at the center of the element and Using the temperature and concentration calculated for each element, the corresponding value of the monochromatic absorption coefficient may be computed by an adequate model: element "k", "

" may be calculated with this coefficient and the emitted radiation from the previous element And finally:

With "N" iterations (k=1,…,N) the radiative intensity emitted at the end of the gas volume calculated. The initialization can be done in the case of a hot plate as shown on .

Note that the total transmission factor of the whole volume is

If spectral measurement are proceeded on a real system, it is possible to find the best parameters for "Tk" and "Ck" to fit the recorded data with the synthetic values. The final result has to be monochromatic radiative intensity calculated at the end of the system the measured one for the same wavelength. Assuming some properties on the shape of the temperature and concentration distribution discards non carbon dioxide lines in the infrared present important temperature dependence around 2400 cm an inverse problem is solved by minimizing the quadratic diff one from the spectroscopic data listed in HITRAN gives the best solution for:

where, "m" stands for T and C parameters vector. Of course, some parameters are necessary to find the best solution in good conditions. The total pressure of the gas must be known, the boundary conditions (hot plate properties, air thickness between the gas volume and the detection system, spectral response of the optical system), and an assumption about the shape's distribution of both temperature and concentration must be included as input data. following diagram describes the principle of the The validity of the monochromatic values calculated

INVERSION OF THE GAS VOLUME TRANSMISSION

The method has been validated on a synthetic gas bench where the concentration and the temperature of gas were controlled values. A FTIR spectrometer was following figure illustrates the experimental setup:

15 th International Symposium on Flow Visualization June 25-28

Using the temperature and concentration calculated for each element, the corresponding value of the monochromatic e computed by an adequate model: " may be calculated with this coefficient and the emitted radiation from the previous element With "N" iterations (k=1,…,N) the radiative intensity emitted at the end of the gas volume

The initialization can be done in the case of a hot plate as shown on Fig. 3 ctor of the whole volume is spectral measurement are proceeded on a real system, it is possible to find the best parameters for "Tk" and "Ck" to fit the recorded data with the synthetic values. The final result has to be calculated at the end of the system for the wavenumber the same wavelength. Assuming some properties on the shape of the temperature and non-physical results and these profiles can be recovered. As previously indicated, carbon dioxide lines in the infrared present important temperature dependence around 2400 cm inverse problem is solved by minimizing the quadratic difference between the measured spectrum and the calculated ic data listed in HITRAN-08 [2] [3] [4] with the developed model.

ere, "m" stands for T and C parameters vector. Of course, some parameters are necessary to find the best solution in he total pressure of the gas must be known, the boundary conditions (hot plate properties, air volume and the detection system, spectral response of the optical system), and an assumption about the shape's distribution of both temperature and concentration must be included as input data. following diagram describes the principle of the optimization code: values calculated has been verified with a specific gas characterization cell.

INVERSION OF THE GAS VOLUME TRANSMISSION SPECTRUM.

The method has been validated on a synthetic gas bench where the concentration and the temperature of gas were spectrometer was used to realize transmission measurements on the gas volume. The s the experimental setup:

International Symposium on Flow Visualization 28, 2012, Minsk, Belarus

Using the temperature and concentration calculated for each element, the corresponding value of the monochromatic [START_REF] Conseil | Spectrométrie Infrarouge de gaz de combustion : application à la détermination de profils de température par inversion de l'équation du transfert radiatif[END_REF]. Then, at the end of an " may be calculated with this coefficient and the emitted radiation from the previous element :

With "N" iterations (k=1,…,N) the radiative intensity emitted at the end of the gas volume may be with the following value:

. spectral measurement are proceeded on a real system, it is possible to find the best parameters for "Tk" and "Ck" to with the the wavenumber (direct problem), and the same wavelength. Assuming some properties on the shape of the temperature and ical results and these profiles can be recovered. As previously indicated, carbon dioxide lines in the infrared present important temperature dependence around 2400 cm -1 . In this spectral range, erence between the measured spectrum and the calculated 08 [2] [3] [4] with the developed model. The minimization code ere, "m" stands for T and C parameters vector. Of course, some parameters are necessary to find the best solution in he total pressure of the gas must be known, the boundary conditions (hot plate properties, air volume and the detection system, spectral response of the optical system), and an assumption about the shape's distribution of both temperature and concentration must be included as input data. Finally, the been verified with a specific gas characterization cell.

The method has been validated on a synthetic gas bench where the concentration and the temperature of gas were o realize transmission measurements on the gas volume. The Fig. 6 : Experiment setup on synthetic gas bench

The gas mixture flows through a pipe form the premix cell to the measurement cell. A heater controls the gas temperature. Firstly, nitrogen, a neutral gas interest) is flowing through the system. U spectrometer) crosses the gas flow and is collected on an MCT detector after modulation in the Michelson interferometer. The result is a reference spectrum for the radiative sour nitrogen in the premix chamber and flows through a pipe to the measurement cell. A new spectrum is acquired and gives . In addition, a "K type" thermocouple is used to give a compar minimization code and the real temperature distribution. A Kulite WCT 312 pressure transducer is also used as shown on Fig. 6 to measure the total pressure of the visualized mixture. The last step gives the spectral transmission of the gas mixture by the calculation of the ratio between the reference spectrum and the spectrum with carbon dioxide: This transmission measurement doesn't require a global calibration because all the source information and the environment influences are integrated into the reference measurement perturbations are avoided since the spectrum is determined with the modulated light coming from the Michelson interferometer. Finally, the energy from the infrared source is higher than the from surrounding because the temperature of the internal source of the spectrometer is around 1400K. radiated energy calculated by Planck's la energy of a real body is always lower than that of a black source of the spectrometer is close to the black the signal to noise ratio is really good in this study and the noticed noise level is about 1%. The data inversion requires a high resolution spectrum ( 4 minutes for a 0.125cm -1 spectral resolution. Moreover, the steady state has to be reached before acquiring data. As an example, for a 5%CO 2 and 95%N 2 mixture during spectra acquisition: setup on synthetic gas bench for measurement of the transmission spectrum flows through a pipe form the premix cell to the measurement cell. A heater controls the gas temperature. Firstly, nitrogen, a neutral gas within the studied spectral range (no absorption in the spectral range of interest) is flowing through the system. Using some optics, the light from the radiative source (internal source of the spectrometer) crosses the gas flow and is collected on an MCT detector after modulation in the Michelson interferometer. The result is a reference spectrum for the radiative source . Then, nitrogen in the premix chamber and flows through a pipe to the measurement cell. A new spectrum is acquired and , a "K type" thermocouple is used to give a comparison point between the result of the minimization code and the real temperature distribution. A Kulite WCT 312 pressure transducer is also used as shown to measure the total pressure of the visualized mixture. The last step gives the spectral transmission of the gas mixture by the calculation of the ratio between the reference spectrum and the spectrum with carbon dioxide: This transmission measurement doesn't require a global calibration because all the source information and the environment influences are integrated into the reference measurement . Moreover all the environment perturbations are avoided since the spectrum is determined with the modulated light coming from the Michelson interferometer. Finally, the energy from the infrared source is higher than the self-emission from surrounding because the temperature of the internal source of the spectrometer is around 1400K. calculated by Planck's law is about 33 times more important at 1400K than at 600K at 4µm. ways lower than that of a blackbody at the same temperature, considering that the infrared to the blackbody assumptions, this ratio is a minimum value. For all these reasons, e ratio is really good in this study and the noticed noise level is about 1%. The data inversion requires a high resolution spectrum (Fig. 2). Therefore, the acquisition time is quite important: about spectral resolution. Moreover, the steady state has to be reached before acquiring data. As an mixture, the monitoring of temperature and pressure gi emperature and pressure for a 5%CO 2 and 95%N 2 mixture Temperature stability is correct and around K: the temperature variation is about 0.6% in the standard Kelvin units. reached and good. The average value is ±1.98%.

For each gas mixture, the average measured temperature value thermocouple's position in the measurement cell. Moreover, as the exact concentration of gas injecte known, the result of the inversion for concentration will be compared to this exact value. identification of temperature and concentration of carbon dioxide for three different cases are presented: 95%N 2 under 1bar then 5%CO 2 / 95%N tests, the initialization values are T = 300 K, calculated and measured values for temperature and concentration:

Molar fraction of CO2 (%) 5,0

Ttck (K) 517,3 Toptical (K) 523,8 1,3

Table 1 : Differences between calculated parameters and real parameters for inversion of transmission measurement

These results confirm both the measurement method and the data processing of the inversion code.

APPLICATION TO CAR ENGINE:

Some measurements have been done on the exhaust pipe of a car engine (DV6 measurement cell used for the validation of the method on the synthetic gas bench is connected to the exhaust pipe behind a particle filter to prevent obstruction and the real position of the cell of the exhaust pipe is shown on the picture:

Fig. 8 : Instrumentation set-up for radiative measurement

In this experiment, the temperature distribution into the pipe is assumed to fit with a Gaussian curve. Thus the temperature can be parameterized with the center position of the Gaussian curve into the pipe, is the temperature difference between the asymptotic value and the center value. The center value is assumed at the center of the pipe . The concentration of carbon dioxide is assumed as constant through the section of the pipe and equal to C x . So the problem needs the identification of four parameters:

The presented results are acquired from thermocouple is 590K and the inversion result at the same position gives 601K less than 2%. The molar fraction has been determined at 7.8% an

th International Symposium on Flow Visualization

June 25-28

Temperature stability is correct and around ±3.6 K which as to be compared to the average K: the temperature variation is about 0.6% in the standard Kelvin units. Concerning the pressure, the stability is also value is 1.01 bar, the absolute variation is ±0.02 bar and the pressure variation is about measured temperature value is compared to the retrieved temperature value at the thermocouple's position in the measurement cell. Moreover, as the exact concentration of gas injecte known, the result of the inversion for concentration will be compared to this exact value. The identification of temperature and concentration of carbon dioxide for three different cases are presented: / 95%N 2 under 1.01 bar and finally 10%CO 2 / 90%N 2 under 0.99 bars. For all these = 300 K, X CO2 = 1%. The following tables give the relative differences between values for temperature and concentration:

ifferences between calculated parameters and real parameters for inversion of transmission measurement

both the measurement method and the data processing of the inversion code.

MEASUREMENTS AND RESULTS

Some measurements have been done on the exhaust pipe of a car engine (DV6-TED4; 1.6 HDI; 110 hp). The measurement cell used for the validation of the method on the synthetic gas bench is connected to the exhaust pipe truction on windows. The experimental setup is described on the following figure and the real position of the cell of the exhaust pipe is shown on the picture:

up for radiative measurement of exhaust gas for a DV6-TED4 Diesel engine

In this experiment, the temperature distribution into the pipe is assumed to fit with a Gaussian curve. Thus the temperature can be parameterized with the equation with center position of the Gaussian curve into the pipe, the width of the curve and the temperature asymptotic value. is the temperature difference between the asymptotic value and the center value. The center value is assumed at the . The concentration of carbon dioxide is assumed as constant through the section of the pipe . So the problem needs the identification of four parameters: , , and C x om a 2500rpm engine speed with 40Nm torque. The temperature value from the thermocouple is 590K and the inversion result at the same position gives 601K: the difference between the two values is less than 2%. The molar fraction has been determined at 7.8% and has been compared to the value measured with an International Symposium on Flow Visualization 28, 2012, Minsk, Belarus value calculated at 600.25 the pressure, the stability is also the pressure variation is about compared to the retrieved temperature value at the thermocouple's position in the measurement cell. Moreover, as the exact concentration of gas injected in the pipe is

The results obtained for the identification of temperature and concentration of carbon dioxide for three different cases are presented: first 5%CO 2 / under 0.99 bars. For all these %. The following tables give the relative differences between ifferences between calculated parameters and real parameters for inversion of transmission measurement.

both the measurement method and the data processing of the inversion code.

TED4; 1.6 HDI; 110 hp). The measurement cell used for the validation of the method on the synthetic gas bench is connected to the exhaust pipe . The experimental setup is described on the following figure

TED4 Diesel engine

In this experiment, the temperature distribution into the pipe is assumed to fit with a Gaussian curve. Thus the . Here is the the temperature asymptotic value. is the temperature difference between the asymptotic value and the center value. The center value is assumed at the . The concentration of carbon dioxide is assumed as constant through the section of the pipe

x .

a 2500rpm engine speed with 40Nm torque. The temperature value from the the difference between the two values is d has been compared to the value measured with an industrial FTIR analyzer equals to 7.5%. For this case, the difference on the concentration obtained between an intrusive measurement and the result of the optical method value is less than 3%. On the left shows the measured spectrum, the inversion result and the difference between these values. On the right side, the calculated temperature distribution and the intrusive value for the engine running point is presented.

Fig. 9 : Inversion results on the transmission measurement for a DV6

The results of the method for this real case are similar to the intrusive data. This study shows that the method is actually relevant for industrial applications. However it has some negative points particularly on the acquisition time and some improvements have to be considered.

CONCLUSIONS AND PROSPECTS

A non-intrusive method has been developed to determine temperature and The selected species is the carbon dioxide infrared bandwidth can be reduced between analysis of spectral data is achieved thanks to a numerical inversion code code is based on an algorithm which uses local minimization of the parameters T and C. It ensures rapid convergence to parameters is limited and it is necessary to limit the number of parameters, a Gaussian parameterization of the searching the value for T end C at each point of the grid. T bench by inversion of transmission spectra. 2.5% with the thermocouple values. The method is applied to carbon dioxide for temperatures under 1000K and near atmospheric pressure but the validity of this method Moreover, the monochromatic absorption coefficient of CO Thus, the method is not limited to the exhaust of car engines and could be appl

th International Symposium on Flow Visualization

June 25-28 industrial FTIR analyzer equals to 7.5%. For this case, the difference on the concentration obtained between an intrusive measurement and the result of the optical method value is less than 3%. On the left shows the measured spectrum, the inversion result and the difference between these values. On the right side, the calculated temperature distribution and the intrusive value for the engine running point is presented.

nversion results on the transmission measurement for a DV6-TED4 car engine at 2500rpm and 40Nm torque

The results of the method for this real case are similar to the intrusive data. This study shows that the method is actually elevant for industrial applications. However it has some negative points particularly on the acquisition time and some

CONCLUSIONS AND PROSPECTS

intrusive method has been developed to determine temperature and concentration distributions of a gas volume. carbon dioxide for which the high resolution spectrum of its between 2380cm -1 to 2400cm -1 where an important thermometric is achieved thanks to a numerical inversion code especially developed for this application. This is based on an algorithm which uses local minimization of the radiative intensity gradient rapid convergence to the best solution. For this method, t it is necessary to assume some particularities about the shape of the a Gaussian parameterization of the T and C distribution has been assumed instead of t each point of the grid. The validation of the method has been done on a synthetic gas ssion spectra. The inversion of these measures showed a temperature difference values. The method is applied to carbon dioxide for temperatures under 1000K and near he validity of this method and of the direct model may be extended to other molecules. he monochromatic absorption coefficient of CO 2 is done with partition functions given Thus, the method is not limited to the exhaust of car engines and could be applied to many other systems. International Symposium on Flow Visualization 28, 2012, Minsk, Belarus industrial FTIR analyzer equals to 7.5%. For this case, the difference on the concentration obtained between an intrusive measurement and the result of the optical method value is less than 3%. On the left side, the following figure shows the measured spectrum, the inversion result and the difference between these values. On the right side, the calculated temperature distribution and the intrusive value for the engine running point is presented.

TED4 car engine at 2500rpm and 40Nm torque

The results of the method for this real case are similar to the intrusive data. This study shows that the method is actually elevant for industrial applications. However it has some negative points particularly on the acquisition time and some concentration distributions of a gas volume. for which the high resolution spectrum of its radiative properties in the where an important thermometric effect exists. The specially developed for this application. This gradient value according to the For this method, the number of identifiable the shape of the distributions. In order to distribution has been assumed instead of he validation of the method has been done on a synthetic gas these measures showed a temperature difference reduced at values. The method is applied to carbon dioxide for temperatures under 1000K and near be extended to other molecules. is done with partition functions given valid up to 3000 K. ied to many other systems.

Fig. 2 :

 2 Fig. 2 : Typical gas absorption lines for carbon dioxide in the infrared.

Fig. 4 :

 4 Fig. 4 : Heterogeneous system divided into N homogeneous and isotherm elements.

Fig. 3 :

 3 Fig. 3 : Typical experimental configuration.

Fig. 5 :

 5 Fig. 5 : Principle of the optimization code.

Fig. 7 :

 7 Fig. 7 : Temperature and pressure for a 5%CO

  ACKNOWLEDGEMENTS.This work is supported by the research project SIMBA funded by "l'Etat "Communauté d'Agglomération du Pays de Montbéliard ». facility CEMBA.REFERENCES1. Al Khoury P, Chavent G, Clément F, Hervé P, and Legras O. transfert radiatif des milieux semi-transparents. Inversion de données spectroscopiques pour le CO RR INRIA 4693, 2003. 2. Rothman LS, Gordon IE, Barbe A, Benner DC, Bernath PF, Birk M, et al. spectroscopic database. Journal of Quantitative Spectroscopy and Radiative Transfer. juin; 3. Simeckov M, Jacquemart D, Rothman LS, Gamache RR, Goldman A. for molecular absorption transitions in the HITRAN database 4. Fischer J, Gamache RR, Goldman A, Rothman LS, Perrin A. the 2000 edition of the HITRAN database 5. Pascale Chelin. Etude des propriétés radiatives de la vapeur d'eau à haute température et diagnostic optique de la combustion H2/O2/N2. PhD thesis, Université 6. P. Hervé. Procédé et système de détermination de la distribution spatiale d'un paramètre milieu semi transparent, Brevet déposé en Mai 2009.

  This work is supported by the research project SIMBA funded by "l'Etat Français" (Contract No.072906524) and by the du Pays de Montbéliard ». Experiments were achieved in the R&D Al Khoury P, Chavent G, Clément F, Hervé P, and Legras O. Etude numérique du comportement de l'équation de transparents. Inversion de données spectroscopiques pour le CO Rothman LS, Gordon IE, Barbe A, Benner DC, Bernath PF, Birk M, et al. The HITRAN 2008 molecular . Journal of Quantitative Spectroscopy and Radiative Transfer. juin; 110(9 kov M, Jacquemart D, Rothman LS, Gamache RR, Goldman A. Einstein A-coefficients and statistical weights for molecular absorption transitions in the HITRAN database. JQSRT, 2006 mars;98(1):130-Fischer J, Gamache RR, Goldman A, Rothman LS, Perrin A. Total internal partition sums for molecular species in the 2000 edition of the HITRAN database. JQSRT, 15 November;82(1-4):401-12. Pascale Chelin. Etude des propriétés radiatives de la vapeur d'eau à haute température et tique de la combustion H2/O2/N2. PhD thesis, Université Paris X, 2003. P. Hervé. Procédé et système de détermination de la distribution spatiale d'un paramètre milieu semi transparent, Brevet déposé en Mai 2009.

Français" (

 ( Contract No.072906524) and by the Experiments were achieved in the R&D Moteur's test Etude numérique du comportement de l'équation de transparents. Inversion de données spectroscopiques pour le CO 2 . Technical Report The HITRAN 2008 molecular 110(9-10): 533-72. coefficients and statistical weights -55. Total internal partition sums for molecular species in Pascale Chelin. Etude des propriétés radiatives de la vapeur d'eau à haute température et haute pression par P. Hervé. Procédé et système de détermination de la distribution spatiale d'un paramètre thermodynamique d'un


X Q Huang

Y Bailly

F Guermeur

V Lepiller

SIMULTANEOUS DETERMINATION OF VELOCITY FIELD AND SIZE DISTRIBUTION OF WATER DROPLETS GENERATED BY AN ATOMIZER FROM IMAGE ANALYSIS

Keywords: Main subjects: particle sizing, velocity Fluid: water spray Visualization method(s): Particle Image Velocimetry Other keywords: polarization ratio, light scattering

A particle sizing method based on the analysis of polarization status of light scattered by targeted particles is introduced. It is combined with Interferometry Laser Imaging for Droplet Sizing method (ILIDS) to determine the size distribution of water droplets in a mist generated by an atomizer at different humidity level under different temperature. The optical setup is combined with Particle Image Velocimetry (PIV) for the simultaneous determination of particle velocity.

Introduction

Particle sizing techniques have been widely investigated and developed in a variety of industries and research laboratories, which aims at characterizing the form and the size distribution of particles ranging from several nanometers to several millimeters in a targeted group. The past several decades has seen the advancement of optical methods attribute to the application of laser techniques. Optical methods, which are based on the analysis of lightparticle interaction, evidently show their advantages by providing non-intrusive and rapid measurements as well as the possibility of combination with velocity measuring devices [START_REF] Black | Laser-based techniques for particle-size measurement: A review of sizing methods and their industrial applications[END_REF] . Fig. 1 illustrates light scattering by a spherical particle. The spherical polar coordinate system is centered on the particle. Theory of Lorenz-Mie provides a solution of Maxwell electromagnetic equations, which elaborates the absorption and scattering of incident light by a sphere in a linear, isotropic and homogeneous medium 2 . A remarkable phenomenon of light scattering is the polarization of scattered light which can be characterized by four experimental parameters including the incident wavelength (λ), the angle of observation (θ), the relative refractive index of the particle and the medium (m), and the particle diameter (d). Thus when the incident light is linearly polarized at 45°, polarization ratio (P) of scattered light is calculated by the following equation:

// // // // ) , , , , ( i i E E i i E E d m P i i s s ⊥ ⊥ ⊥ ⊥ = ⋅ = = θ λ Eq.1
where s E and i E represent respectively the electric field of scattered light and that of incident light. ⊥ i is the scattered irradiance per unit incident irradiance given the incident light is polarized perpendicular to the scattering plane, while

//
i is its counterpart when the incident light is polarized parallel. Based on theory of Lorenz-Mie, the polarization ratio as well as the scattered irradiance per unit incident irradiance of scattered light by water droplets of different diameters is calculated by a program developed in Mathematica, as shown in Fig. 2. The water droplets are dispersed in air, with a relative refractive index (m) of 1.33+0i. The incident laser light is double-pulsed laser Nd:Yag at a wavelength (λ) of 0.532µm. The measurements are performed at an observation angle (θ) of 66°.

Fig.2 Polarization ratio and scattered irradiance of water droplets in air

On the other hand, the two polarized components of scattered light can be determined experimentally from the analysis of particle images respectively polarized parallel and perpendicular taken by a high resolution CCD camera. Thus the particle diameter under investigation can be obtained from inverse calculation. From Fig. 2, fluctuation of the theoretical values of polarization ratio for the particles bigger than 0.5µm is observed, which causes uncertainty for the inverse calculation. Meanwhile the scattered irradiance presents an increasing tendency generally monotone with particle diameter, which provides a good reference for the inverse calculation of particle size.

Since the particles investigated seeding the flow, the experimental setup becomes a PIV device which enables the simultaneous acquisition of velocity profile when a double-pulsed Nd:Yag laser is used as incident light source. The two images illuminated by two successive laser pulses are correlated to render particle velocity in the observed field. Theoretical work and a more detailed introduction of the method of polarization ratio can be found in the authors' previous publication [START_REF] Huang | Particle sizing and velocity measurement of microspheres from the analysis of polarization of the scattered light[END_REF] .

In the current study, method of polarization ratio is combined with Interferometry Laser Imaging for Droplet Sizing method (ILIDS) to determine the size distribution of water droplets in a mist generated by an atomizer. Interferometry Laser Imaging for Droplet Sizing method (ILIDS) ILIDS is a developed optical method for the determination of droplet size bigger than 10µm. According to Van de Hulst [START_REF] Van De Hulst | Glare points[END_REF] , when a water drop is illuminated by a wide beam, glare points can be viewed and imaged from a certain direction, respectively corresponding to the intensity maxima of the rays reflected and refracted in the forward-scatter region, around 20° to 70° 5 . Glantschnig and Chen [START_REF] Glantschnig | Light scattering from water droplets in the geometrical optics approximation[END_REF] showed that regular fringes arise from interference between reflected and refracted rays of orders p=0 and p=1, as shown in Fig. 3. The angular frequency of fringes is to droplet diameter, as illustrated by the following formula 8 :

            - + + ∆ = 2 cos 2 1 2 sin 2 cos 2 2 θ θ θ θ λ m m m N d Eq.2
where N is the fringe number, and λ is the wavelength of the light incident. θ is the scattering angle and ∆θ is the collecting angle, while m is the refractive index of the droplet. fringe-intensity within a circle has the form of a periodic function that resembles a sine wave, the number of peaks, i.e. the fringe number, was obtained by fitting the intensity information to a parametric model. Thus the particle diameter can be calculated through Eq.2.

Experimental setup An optical test rig has been set up to perform the acquisition of particle images, as illustrated by the sketch in Fig. 5. Light scattered by the particles investigated is collected at a certain angle θ, then the two components polarized are separated by a cube polarizer and arrive at the camera after a couple of reflections on the mirrors. Particle images polarized respectively parallel and perpendicular are taken by a CCD camera (TSI 630062: 4008×2672 pixels). According to the basic principle of ILIDS, defocused images have been taken for the sizing of water droplets. 6 shows a defocused image of water droplets generated by an atomizer, with the image perpendicular polarized at left and that parallel polarized at right. For the experimental results presented, measurements were carried out at a scattering angle of 66°. Double-pulsed Nd:Yag laser was employed at a wavelength of 0.532µm, which was inclined at 45° to render an incident beam linearly polarized at 45°. With a 16-bit grayscale CCD camera of 4008×2672 pixels, the size of each pixel was 9×9 µm 2 . The size of measuring area was 4×4 cm 2 and the magnification ratio the optical system is 1:10. Thus, the diameter of a droplet producing exactly one fringe across the observed circle would be 15.7µm at f2.8. To determine the particle size, the fringe number of each detected particle was calculated. Those with a fringe number less than 1 were taken out to perform the calculation of polarization ratio, based on which their diameters were obtained. Measurements were carried out to determine the size distribution of water droplets generated by an atomizer under different heating temperature at different humidity level. Experimental results shown in Fig. 7 were drawn based on 30 images successively taken by CCD camera. Particle size decreased evidently when the water droplets were heated under a higher temperature, whereas larger particles became more important in the mist when humidity level was increased.

Fig. 7 Size distributions of water droplets in a mist generated by atomizer under different conditions

In the current study, particle velocity was obtained by Particle Image Velocimetry from the analysis of two successive images of particles illuminated by double-pulsed laser sheet. Here we take the image polarized parallel (left half of Fig. 6) for the determination velocity field. Le flow direction of the water droplets ascending naturally is influenced by the hot air flow blowing downward, as shown in Fig. 8.

Conclusion

The optical method of polarization ratio provides a non-intrusive particle sizing technique for spherical particles in twophase flow under low concentration. The combination with velocity measurements (PIV) makes it possible to determine the particle velocity simultaneously. The optical setup also allows the combination with other particle sizing methods for the application to particles of a larger size range.

Fig. 1

 1 Fig.1 Light scattering by a spherical particle

  Symposium on Flow Visualization June 25-28, 2012, Minsk, Belarus ISFV15 -Minsk / Belarus -2012

Fig. 3

 3 Fig.3 Principle of Interferometry Laser Imaging for Droplet Sizing method (Albrecht 7 )

Fig. 4

 4 Fig.4(a) presents a defocused snapshot of water droplets generated by an atomizer. The fringe number of each particle corresponds to different particle diameter. A circular Hough transform was applied to locate circles in the image from edge-strength information, which rendered the size and position of each circle detected, as shown in Fig.4(b). Since the

Fig. 5

 5 Fig.5 Sketch of optical setup

Fig. 6

 6 Fig.6 Defocused images of water droplets parallel and perpendicular polarized

Fig. 8

 8 Fig.8 Velocity field of water droplets under investigation

Acknowledgement This work is supported by the research project SIMBA funded by "l'Etat Français" (Contract No.072906524) and by the "Communauté d'Agglomération du Pays de Montbéliard ».


S Villalonga

C Thomas

C Nony

F Thiebaud

M Geli

A Lucas

K Knobloch

C Maugy

F Nony

K Kremer- Knobloch

Applications of full thermoplastic composite for type IV 70 MPa high pressure vessels

Keywords: Thermoplastic Composite, Type IV High Pressure Vessel, Hydrogen Storage, Fuel Cell Vehicles, Hydrogen Energy, New Energy Technologies

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

General Introduction

The fossil resources shortage and climate change both encourage an evolution of our energy supplying system. For more than a decade, hydrogen as an alternative to traditional energy sources has been focus of research and development efforts. Hydrogen, combined with efficient converters such as fuel cells, represents a very promising way to realize a more sustainable energy system. Nowadays, for automotive applications, even if several interesting prototypes have been displayed, no storage system meets every international standards and car manufacturers' requirements. However, developments on 700 bar type IV vessels with Thermoplastic PA6 liner [START_REF] Nony | Type IV 700 bar-vessel for compressed gaseous hydrogen storage: Material research and performance achievements[END_REF] [START_REF] Villalonga | Composite 700 bar-vessel for on-board compressed gaseous hydrogen storage[END_REF] have demonstrated very encouraging results (cycling resistance, burst pressure, H 2 tightness, gravimetric and volumetric storage capacities, H 2 barrier properties). Over the last decade, the CEA (French Atomic Energy and Alternative Energies Commission) has developed, with several industrial and academic partners, valuable knowledge and competences focused on reactive rotational molding of high performances polymeric liners. The CEA has improved design, calculation, manufacturing and testing of structural composite vessel shells too. This paper reviews the most recent technical and scientific achievements in comparison with the stateof-the-art of CGH2 storage systems, especially results from the French collaborative research Hype project. Vessel composite shells are manufactured by filament winding using Carbosatmp TM (Toray T700 fiber, PA6 thermoplastic matrix). This paper presents the material used, its properties, composite and vessel design and calculation, manufacturing aspects, testing and automotive application.

Material

Considering gravimetric and volumetric densities, as well as safety requirements, type-IV vessels are the most promising solutions for on-board hydrogen storage. They display no fatigue issues, the H 2 leakage is far below standards, and they show decent weight performances. CEA PA6 liner technology [3][4] has demonstrated abilities to withstand 15,000 cycles from 20 to 875 bar, as required by European or ISO standards for a nominal operating pressure of 700 bar. Moreover, the CEA is involved in the development of innovative materials such as materials for composite filament winding (thermoplastic and thermoset matrices) and materials for different types of liner (based either on thermoplastics and thermosets). In the scope of Hype project, the composite shell material, Carbostamp TM 'Fig. 1', is a polyamide matrix (polyamide 6) reinforced with T700 Carbon Fibers. Its fiber volume fraction is about 50 %. The material used for the liner is a PA6 thermoplastic material made by reactive rotomolding process so that the vessel is a PA6 monopolymer vessel with a very good consolidation between the liner wall and the composite shell. Thermoplastic matrix composites can be considered as relevant solutions since they present many advantages. They show numerous potential interests in term of storage, shell fiber volume fraction quality, impact resistance, durability and recyclability. Indeed, recyclability of CF/PA6 composite is easier than CF/Epoxy composite, which is very attractive for automotive industry. CF and PA6 can be separated by considering PA6 fusion temperature and re-used to make short carbon fibers for example. Moreover, CF/PA6 process could be less expensive and faster. Compared to epoxy prepreg which has to be stored at -20°C and used before 6 months, there is no special constraint to store CF/PA6 before the filament winding. Compared to wet filament winding, there is no constraint of impregnation during the winding, no time and cost wasted to prepare resin, no pot-life time, no impregnation process parameters which need to be adjusted, no volumic fiber rate variation of the composite in the ply or in the laminate of the vessel, no time and cost wasted to clean the filament winding machine, no safety constraint (equipments and workers) as using volatile resin, no time and cost wasted due to the need of thermal curing after the filament winding process. The main difficulty of the CF/PA6 process (impregnation / consolidation / curing) is the consolidation step.

Material Properties

The design of reinforced laminated composites for type IV high pressure vessels application takes into account the initial mechanical properties of the material including safety coefficients [START_REF] Thomas | Study of damages on thermoplastic composite structures: application for hydrogen storage vessels under high pressure[END_REF].

Carbostamp TM PA 6 is a polyamide matrix (polyamide 6) reinforced with T 700 Carbon Fibres from Toray SOFICAR. The fibre volume fraction, determined by pyrolysis at 500°C, is about 50 %. The matrix is semi crystalline and presents a melting point at 215°C. The density, measured according to Archimede's principle by weighing the samples in air and a liquid of known density (water), is about 1.47 kg/m3 for Carbostamp TM PA6.

The first characterizations are dedicated to the mechanical behavior of the material. Tensile tests have been conducted on unidirectional specimen with three different orientations 0°, 90° and 45° (orientations of fibres compared to the tensile load direction) and cross-plied laminates [±45]s. The test samples are manufactured by hot compression moulding. The materials and the mould are heated to reach a temperature above the melting temperature of the matrix. Then a pressure is applied. Finally the mould is cooled still under pressure to a temperature below the glass transition temperature before that the consolidated specimen is extracted.

The next mechanical experiments deal with different damages occurring in carbon fibers / polyamide matrix composite structures in the aim of understanding their influence on the mechanical properties. Transverse cracking and delamination have been studied 'Fig. 2'. The results are described in the following reference: [START_REF] Thomas | Research and achievements on carbon fiber reinforced thermoplastic composites for high pressure storage[END_REF].

Composite and Vessel Design and Calculation

The design of a composite vessel is extremely specific as it should take into account boss, liner and composite materials, liner manufacturing process specifications as well as constraints for the filament winding process. The composite architecture has to be optimized in order to approach the predefined target as closely as possible. Moreover, the design should take into account the service and test pressures, the external stresses, which are specific to the use (and dependent on, for instance, impact, aggressive media, temperature etc…). For the specific case of type-IV vessels, the main composite design test is the burst pressure.

Design of the dome shape region is one of the main difficulties for a high pressure vessel. The shape dome has to be compatible with the filament winding process. During the 18 th century, the French mathematician Alexis Claude Clairaut (1713-1765) 'Fig. 3.' studied the geodesic shape on the earth surface [START_REF] Clairaut | Theorie de la figure de la terre[END_REF] and developed a theorie. Today, this theorie is used to design and manufacture high pressure vessel by filament winding 'Fig. 4.".

The best estimate calculation code for design evaluation has to be enough qualified to minimize result uncertainties and optimized High Pressure Vessels. In each case, HPV should be optimized to hold a maximum of fluid with a limited volume at the best manufacturing costs, load bearing capability and lifespan. This previous step is vital before enlarging design parameters and finally answering industrial and economic stakes.

Current finite element calculation limitations are related to different models and behavior uncertainties. We can quote for instance the following aspects: -

The lay-up design about the winding evolution strategies on end-closures revolution areas (thickness evolutions, angles evolutions with geodesic winding trajectories or with non-geodesic winding techniques, laminate superposition, slippage parameters)

-

The finite element models about the choice of the elements types. Shell elements are mainly suitable for analyzing thin to moderately-thick shell structures whereas solid elements are more appropriate for calculation precision in spite of their more laborious use for winding mesh procedures, which could appears to be really time consuming in an optimization procedure where altering lay-ups is predominant. Another limitation concerns the accuracy of axisymmetrical models to faithfully characterize the composite cyclic periodic behavior.

-

The damage models. An optimization procedure has to take into account every damage mechanism (fatigue, failures…) occurring in the structure. The fatigue behavior for composite materials relies on the nature of the constituents, the process parameters and mainly on the application. -Numerous damage models for composites structures are available. For example, the most sophisticated one are based on micro-models [START_REF] Mishnaevsky | Micromechanical modeling of damage and fracture of unidirectional fiber reinforced composites: A review[END_REF] but are limited by their complexity and costs to relatively basic structures. Moreover, macro models [START_REF] Spottswood | Progressive failure analysis of a composite shell[END_REF], certainly more appropriate for more complicated and large structures, exist but still remain to be tested on pressure vessels.

Today, the CEA is working on the qualification step of the FE calculation code for Type IV HPV design. For Hype vessel design, the laminate has been defined by using Abaqus Wound Composite plugin. The stress rupture analysis is compute by Abaqus code 'Fig. 5'.

Vessel Manufacturing

The liner is manufactured by reactive rotomolding process. The CEA developed, with it partner AIR LIQUIDE, a polyamide-6 (PA-6)-based polymer. A specific formulation was chosen based on the best compromise for permeation, mechanical strength, elasticity and processability. Since this polymer was sensitive to thermal oxidation, the CEA created an innovative 1-step reactive rotational molding process [3][4] starting directly from liquid monomers or precursors. The mold was considered to be a chemical synthesis reactor in which the polymerization and the shaping/molding occurred simultaneously. The polymerization took place between 150 and 180°C and lasted from 30 to 90 min depending on the size of the pieces 'Fig. 6' and eventually on the number of layers.

The composite shell is obtained by filament winding process. The filament winding process developed in this project presents a preheater (infrared), which raises the temperature just below the melting point, and a heater (infrared), which raises it above the melting point near the mandrel. The consolidation is ensured by a compaction roller. An adjustable tension is magnetically applied to T700/PA6 tape spools. Other heating systems have been evaluated as hot air system, laser system. Several PA6 liners with internal liner volume of 2 H 2 liters, 37 H 2 liters 'Fig. 7' have been winded with Carbostamp TM PA6 with a good composite consolidation. 105 H 2 liters will be winded during july 2011.

Vessel Testing

One of the main issues of project development is the characterization of vessels. Each CEA designs are validated in full-scale 700-bar H 2 cylinders. Through liner polymer syntheses and formulation knowledge, reactive roto-molding and filament process capabilities, design and simulation methodologies, the CEA and its French partners have developed and fabricated a full PA6 thermoplastic type IV 700-bar vessel able of meeting the Hype project test criteria with regard to burst (>1,400 bar), permeation (<< 1 Ncm 3 /h/L) and cycling (>15,000 cycles). H 2 standards are changing and safety factors are decreasing. Today, for the next ISO standard, considering automotive application, the burst pressure ratio could be 2.and, for the permeation, the maximum rate could be 6 Ncm 3 /h/L. At this stage of the project, MaHyTec partner is in charge of burst and cycling tests. The 37 L vessels will be tested in june 2011 and 105 L vessels in august 2011.

Automotive Application

High pressure hydrogen storage (700 bar) is one of the different hydrogen storage ways studied by PSA Peugeot Citroën, for its Fuel Cell vehicle. The FISYPAC Fuel Cell hybrid vehicle developed by PSA in 2009 'Fig. 8' was already equipped with high pressure 700 bar H2 type III tanks. One of the main challenges for high pressure H2 storage system is to increase system's volumetric and gravimetric capacities. HYPE tanks are aimed at answering those requirements.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 12345 Fig.1. Toray-Soficar Carbostamp TM .

Fig. 8 .

 8 Fig.8. FISYPAC Fuel Cell hybrid vehicle from PSA

Acknowledgment : This work has been supported by French Research National Agency (ANR) through "Plan d'Action National sur l'Hydrogène et les piles à combustible" program (ANR-07-PANH-006-02).


Ali Zeaiter
email: ali.zeaiter@femto-st.fr

David Chapelle

Philippe Nardin Femto-St

Searching out the hydrogen absorption/desorption limiting reaction factors: strategies allowing to increase kinetics

Keywords: Hydrogen, absorption, desorption, hydride materials, thermodynamic model, kinetics

Hydrogen gas has become one of the most promising energy carriers. Main breakthrough concerns hydrogen solid storage, specially based on intermetallic material use. Regarding the raw material abundance and cost, the AB type alloy FeTi is an auspicious candidate to store hydrogen. Its absorption/desorption kinetics is a basic hindrance to common use, compared with more usual hydrides. First, discussions based on literature help us identifying the successive steps leading to metal hydriding, and allow to introduce the physical parameters which drive or limit the reaction. This analysis leads us to suggest strategies in order to increase absorption/desorption kinetics. Attention is then paid to a thermodynamics model, allowing to describe a macroscopic solid storage reactor. Thus, we can achieve a simulation which describes the overall reaction inside the hydrogen reactor and, by varying the sub-mentioned parameters (thermal conductivity, the powder granularity, environment heat exchange…), we attempt to hierarchy the reaction limiting factors. These simulations are correlated to absorption/desorption experiments for which pressure, temperature and hydrogen flow are recorded.

1-Introduction

Solid storage coupled with intermetallic is one of the most promising technologies that can provide an efficient and secure holding of the hydrogen gas comparing with the traditional ways such as gas and liquid storage. The solid storage method can be accomplished with a variety of structural type of intermetallic like (AB, AB 2 , AB 5 , AB 3 and A 2 B 7 ). By inspecting deeply the different structural types, the FeTi AB intermetallic has received a special attention because of the plenty in raw material and the reversible characteristic of the hydride reaction, in other word: absorption and desorption processes. However, in order to ensure that the hydride reaction meets engineering requirements, this alloy still need to be characterized to overcome all challenges such as a harsh activation behavior or a slow reaction kinetics. A pure FeTi needs an applied pressure in the range of 30 bar and a temperature in the range of 670 K to be activated, and it was reported that kinetics are very slow (an absorption stage of a fully activated sample can take over than 10 hours to reach the maximum of storage capacity ~ 1.8 wt.%). Focusing on the hydride reaction of FeTi alloy, this reaction is described according to the following chemical equation: nFeTi + H 2 ↔ nFeTiH 2/n ∓ ∆ , where ∆ is the heat generated or needed during absorption and desorption respectively, because as well known, while the absorption stage the chemical reaction is exothermic, contrariwise the reaction is endothermic for the desorption stage.

All physics phenomena should be described mathematically by using concrete equations which allow to find the system response based on initial conditions, i.e. input variables and material properties for the hydrogen absorption/desorption process. A thermodynamic model is developed based on several physics principles: it describes the reaction from a macro-scale point of view. One of the most relevant issues on the reaction rate is the thermal management and specially the overall heat transfer coefficient. It has been reported that the thermal conductivity of the powders decrease in correlation with the number of repeating cycle to reach a very low value in the range of 0.1 W/mK. The thermal conductivity can be mended by several ways such as powders compacting or/and the intercalation of some materials that have an important thermal conductivity. Those metals can be copper wire, porous graphite, aluminum foam, etc. K wang et al. [1] used the porous graphite with powders compacting to enhance thermal conductivity and have reported that λ is increased to reach 3W/mK. Also several works have been reported on material doping, as an example P. Mandal et al. [2] have worked on the FeTi (Mn) which have 46% Fe, 50% Ti, and 4% Mn, the alloy is mixed with Mg according to the following composition 40% (FeTi Mn) and 60% Mg. They have shown that 80% of the total amount of hydrogen is desorbed within 10 min. In the other hand, fabrication methods can be listed as one of the pre-treatment ways to increase kinetics. H. Aoyagi et al. [3] have monitored the behavior of a FeTi sample, and as a result from the performed mechanical treatment, it was seen a formation of new activated surfaces and an improvement in the absorption/desorption reaction rate. Another example of fabrication method is the SICS (self-ignition combustion synthesis). It was proposed by NaotoYasuda et al. [4] for the following interests: purifying the product, hydriding an alloy directly and reduction of the hydrogenation time. It was seen also from literature that powder morphology or grain size is one of the major parameters that can be monitored in order to increase kinetics.

2-Hydriding modeling and experimental comparison

The present section is devoted to the thermodynamic model description. The model aims to describe the (de)hydriding reactions occurring inside a reactor containing a defined amount of hydride powders, fully activated and ready to start an absorption/desorption cycle. The hydride tank has a cylindrical shape which is coupled with a heat exchange (water flow on the outside surface). Based on the principles of energy balance, mass balance and ideal gas law equation, equations for both the absorption and desorption stage are proposed by B.A.Talaganis et al. [5]. These equations are listed in Table 1. The simulation is carried out assuming a set of assumptions: constant applied pressure during absorption and desorption, constant cooling or heating fluid temperature, uniform distribution of temperature inside the tank. The finite difference method is employed, i.e. all equations are discretized in time only, alloying to determine according to time: tank temperature, reaction rate, and hydride mass formed during reaction.

Absorption Desorption

Energy balance + = ( ( - + ( --∆ (1)
+ = ( !" ( - + ( -+ ∆ (6) Rate reaction = # $% & '( ( ) ) #* (+ - (2) = # $% & ( ) , #* #* ( (7)
Ideal gas law ) = & .. ( 3
) ) = & .. ( 8
)
Equilibrium pressure #* = # ( ∆ & , ∆/ & ) 0 (4) #* = # (, ∆ & 1 ∆/ & ) 0 (9)
Hydride mass

= (5) = (10)

Table 2-Input variables and their values

In this study, attention is paid on the behavior of the FeTi-X when 100g of this intermetallic powder is filled into the cylindrical tank and subjected to 10 bar of constant hydrogen pressure. The external heat exchange is natural convection. Figure 1 shows hydride reaction rate and temperature response during absorption. Same profile is observed for both temperature and absorption rate evolution. Simulations predict an extremely fast response of the system. Maximum temperature and absorption are reached before 100s, then both quantities decrease rapidly and reaction ends after 10 minutes. The pick of absorption rate occurs a bit earlier, after less than 1 minute.

Having a simulating tool to predict numerically the FeTi-X alloy response during hydriding, turns out to be really advantageous, in order to understand the various steps of the reaction. However, in order to make it efficient, it is required to have a comparison with experimental results. To do so, a 100 g of FeTi-X is filled into a hydrogen reactor coupled with a thermocouple, and putted under the following thermal conditions: cooling and heating is done under room temperature (natural convection with a coefficient of convection of 10 W/mK). Fig- 2 present schematically all the components of the bench test that has been used to perform this investigation.

Figure. 2-Hydriding reactor and bench test

The FeTi-X sample is first subjected to a pressure of 20bar (constant according to time). It is observed a small increase in temperature after 3000 min. Then the powder is subjected to vacuum, and again a pressure of 60bar is applied in order to force the activation mechanism. In this stage, a higher amplitude temperature pick is observed but still requires 1350 min to occur. The powder is subjected once again to vacuum, and finally a pressure of 10bar is applied twice, allowing to guaranty the complete activation of the sample and to ensure repeatability. The two temperature curves beneath in Fig. 3 show a comparison between experimental result (on top right) and simulation result (on top left). Simulation result confirms the overall tendency, i.e. a fast increase of temperature with comparable maximum value, followed with a slow cooling till ambient temperature. Retrieving the accurate properties of FeTi-X in order to improve the good correlation between the model and the experimental results will require to describe more precisely physical phenomena and to lighten hypothesis.

3-Sensiblity analysis of the thermodynamic model

Regarding mathematical models, they are most of the time of complex form and involve a huge numbers of input parameters. It's so important to classify input parameters by order of their impact on the output data. The prioritization allows to define on what input parameters attention is to be paid. This classification is here obtained using a sensibility study based on screen design methods. Doing so, all input factors are ranked in order of influence, but from a qualitative point of view. In this paper, Morris design method is used to manage the sensibility analysis on the thermodynamic model. This method imposes variation of one input factor at a time. Sensitivity is related to two factor indicators, proposed by Morris method: µ that evaluates the overall effect of the input on the output and σ that gives the coupling effect between input factors. In order to perform a global sensitivity approach, elementary effect of each input must be identified as follow: (13) Where CD 3 is a small perturbation induced to an input in the direction of the unit D 3 vector.

Calculation of the sensibility analysis is done using dedicated in-house Matlab® software, called SIMBAD. Here after, all input and output parameters used in this study are listed:

1-Input parameters: T in , A, T wa , ∆H a , ∆S a , m s , Ea, Pa, C a , U, Cp H2 , and V.

2-Output parameter:

(absorption reaction rate) In Fig. 4, all input variables are ranked according their elementary effect mean value and standard deviation. This analysis leads to the conclusion that activation energy Ea has the most important nonlinear effect, i.e. it is coupled with other input parameters. Then, ∆H a , ∆S a come after and their effect is also nonlinear. As a general explanation, variables that have the greatest mean value are the most sensible parameters and depending on the standard deviation value, effect can be linear or nonlinear: the more Std(2 3 is great, the more non linearity effect increases. From Fig. 4, we can rank input parameters by order of importance as follow:

Ea> ∆S a >∆H a >T wa >C a and other inputs have negligible effects. Consequently, improving the accuracy of predicted hydride tank responses requires determining with care values of these parameters but also their evolution with the temperature.

4-Conclusion

The present paper aims to contribute at developing strategies in order to increase absorption/desorption of Fe-Ti based alloys. We present a macro scale thermodynamic model, that simulates the response of a FeTi-X hydride tank, and in parallel we introduce a bench test which is employed to establish a comparison between experimental and numerical results. It has been clearly seen that material properties have to be more precisely identified in so far to ensure more accurate fitting between experiments and simulation results. Finally, a sensibility analysis is performed to evidence the most influent input variables: i.e. variables that change widely output results, such as reaction rate. In order to get a robust model and excellent fitting between simulation and experiment, a model calibration must be performed and lighter hypothesis introduced.

Figure 1 -Figure 3 -

 13 Figure 1-Simulation results (reaction rate and temperature evolution according to time of FeTi-X)

Figure 4 -

 4 Figure 4-Morris design sensitivity indicator

Table 1 -

 1 Absorption/desorption equationsIn Table2, all input parameters are listed with their values.

	Variables	Definition	FeTi-X	LaNi 5
		Heat capacity of H 2	14.3	14.3
		Solid mass	100g	100g
		Heat capacity of solid	0.468 J/g.K	0.355 J/g.K
		Over all heat transfer coefficient	240 W/mK	240 W/mK
		H 2 inlet temperature	296 K	296 K
	∆	Cooling water temperature Formation enthalpy of absorption Pre-exponential factor	296 K -27.5 kJ/mol 16.4 1/s	296 K -30.478 kJ/mol 59.2 1/s
	% ) ∆/	Activation energy of absorption Absorption pressure Formation entropy of absorption Atmospheric pressure	23.8 kJ/mol 10 bar -104.32 J/mol.K 298 K	21.17 kJ/mol 10 bar -108 J/mol.K 298 K
	∆ sc	Heating water temperature Formation enthalpy of desorption Stoichiometric coefficient	296 K 27.6 kJ/mol 1	296 K 30.8 KJ/mol 2.76
		Pre-exponential factor	2.6 1/s	9.6 1/s
	% ∆/	Activation energy of desorption Formation entropy of desorption Molar mass of formed hydride	19.87 kJ/mol 103.04 kJ/mol 103.7 g/mol	16.42 kJ/mol 108 J/mol.K 432 g/mol
	Pd	Applied desorption pressure	0.068 bar	0.068 bar

Acknowledgements

Authors would like to thanks Scott Cogan, researcher at the Applied Mechanics Department of the FEMTO ST Institute, to offer them the opportunity to use the SIMBAD Matlab® software. It is dedicated to the sensitivity analysis but also to model optimization.