Stéphane De Chaisemartin
email: stephane.de_chaisemartin@em2c.ecp.fr

Lucie Fréret
email: lucie.freret@em2c.ecp.fr

Frédérique Laurent
email: laurent@em2c.ecp.fr

Marc Massot
email: massot@em2c.ecp.fr

Corine Lacour
email: corine.lacour@em2c.ecp.fr

Anne-Laure Birbaud

Sébastien Ducruix
email: ducruix@em2c.ecp.fr

Daniel Durox
email: durox@em2c.ecp.fr

Pulsated free jets with spray injection : Eulerian multi-fluid modelling and simulation versus experimental measurements

Keywords: Liquid Sprays, Polydispersity, Multi-fluid models, Spray Equation, Pulsated round free jets

The accurate simulation of the dynamics of polydisperse evaporating sprays in unsteady gaseous flow with large scale vortical structures is both a crucial issue for industrial combustion application and a challenge for modelling and scientific computing. Various appraoches have been developped in order to resolve the dispersed liquid phase starting from a "mesoscopic" level of description through the Williams equation. The usual Lagrangian approaches in polydisperse unsteady configurations lead to a very high computational cost and induce coupling difficulties due to the different type of description of the two phases involved. Several Eulerian models have been recently developed in order to avoid these two difficulties. However, the validation of the obtained models through detailed comparisons with measurements on well-instrumented experimental devices are seldom conducted in the litterature. In the present contribution, we consider the configuration of acoustically pulsated free jets with a polydisperse spray injection in a 2D axisymmetrical pulsated jet. It creates large vortical structure which are representative of the flow in more complex configurations and a strong interaction with the injected spray. In this context, we provide both a series of detailed experimental measurements through the coupling of laser diagnostics, as well as the multi-fluid modelling and associated numerical schemes in order to simulate accurately the dynamics of the spray. It allows us to study the physics of the coupling between the two phases and to validate the Eulerian multi-fluid model using the experimental measurements.

Introduction

In many applications, ranging from combustion chambers up to solid propellant rocket boosters, two-phase flows, where the carrier phase is a gaseous flow with large scales vortices and the second phase is a dispersed liquid one constituted of a cloud of droplets, play a significant role. Describing the dynamics of such flows is crucial for both the physical understanding and the development of predictive numerical tools. It relies on the ability to capture the dynamics and evaporation of droplets of various size.

In the present contribution we consider a mesoscopic point of view and describe the droplets as a cloud of point particles called a spray. This level of description is the only one for which numerical simulations at the scale of a combustion chamber or in a free jet can be conducted. The principal physical processes that must be accounted for are : transport in real space, droplet evaporation, acceleration of droplets due to drag conditionned by size, and eventually coalescence of droplets leading to polydispersity. The major challenge in numerical simulations is to account for the strong coupling between these processes. [START_REF] Williams | Spray combustion and atomization[END_REF] proposed a relatively simple transport equation based on kinetic theory that has proven to be useful for treating polydisperse, dilute and dense liquid sprays. In the literature, the Lagrangian Monte-Carlo approach is generally considered to be the most accurate for solving Williams equation. However, its computational cost is high, especially in unsteady configurations, and the method is difficult to couple accurately with Eulerian descriptions of the gas phase. There is thus considerable impetus to develop Eulerian methods for treating Williams equation. In this paper, we limit our attention to laminar flows and thus no turbulence models are required to close the spray equation.

In a recent paper Laurent et al. (2004a,b); [START_REF] Fox | Numerical simulation of polydisperse, dense liquid sprays in an eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF] have demonstrated the capability of an Eulerian multi-fluid model to capture the physics of polydisperse evaporating sprays. This approach relies on the derivation of a semikinetic model from the Williams equation under the assumption of mono-kinetic number density functions marginals at fixed droplet size. The idea is to keep size as a discretized internal variable since the dynamics and evaporation depend essentially on droplet size, thus yielding conservation equations for velocity moments of order zero and one, in fixed droplets size intervals, of the kinetic number density function. However, such models, even if they have proven their S4_Tue_C_23 6 th International Conference on Multiphase Flow, ICMF 2007, Leipzig, Germany, July 9 -13, 2007 efficiency in the framework of comparisons between numerical simulation with other types of solvers such as Lagrangian solvers [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF], have seldom been compared to experimental measurements in well-controlled configurations which can mimic the difficulties encountered in realistic devices.

Consequently, the purpose of the present study is two-fold. Based on a recent study of the dynamics of pulsated free jets [START_REF] Birbaud | Dynamics of free jets submitted to upstream acoustic modulations[END_REF]) and on experimental studies of reacting two-phase flows [START_REF] Pichard | Effects of acoustic modulations on a spray flame[END_REF], we consider the configuration of a polydisperse spray injected in a nonreacting pulsated laminar round jet. The first part of the study is the design of the associated experimental configuration and experimental set-up which are presented in section 1 and 2. We then present the coupled measurement techniques involved in this study (section 3) in order to be able to characterize the physics of such a phenomenon. The characteristics of the flow and spray are then analyzes in section 4 and the phase-locked experiments presented in section 5. The experimental results non only allow a study of the physics of the proposed configuration, but also provides the possibility of validating the Eulerian multi-fluid model by comparing numerical simulations to their measured counterparts. This second aspect is based on the recent developments of the numerical tools for the simulation of Eulerian multi-fluid models [START_REF] De Chaisemartin | Evaluation of Eulerian Multi-Fluid versus Lagrangian methods for the ejection of polydisperse evaporating sprays by vortices[END_REF]) and on an extension of such tools, in the present contribution, to 2D axisymmetrical configuration. The Eulerian multi-fluid model is thus presented in section 6 and the new dedicated 2D axisymmetrical numerical scheme is introduced in section 7, where its numerical validation is conducted in a simple test case. We are then able, in a last section, to compare the experimental measurements to the numerical simulations in this interesting and challenging configuration.

Configuration and modeling

There exists only a few comparative studies between numerical simulations and experiments in multiphase flows in perfectly controlled situations. An interesting configuration is that of a laminar, unsteady multiphase flow, with large stretches and high levels of vorticity. The dynamics of the droplets and the dynamics of the gas flow can then be clearly observed. The configuration which seems the simplest to set up is that of an excited jet, laden by droplets. One knows, since [START_REF] Crow | Orderly structure in jet turbulence[END_REF], that a circular jet has a natural instability at a frequency corresponding to the preferred mode of the jet. This instability can be characterized by the Strouhal number St , defined by:

St = f D U b (1)
based on the frequency f , on the jet diameter D and on the mean axial bulk velocity U b . According to [START_REF] Gutmark | Preferred modes and the spreading rate of jets[END_REF], the Strouhal number of the preferred mode ranges between 0.2 and 0.64, depending on the experimental conditions, but it is more generally in the range 0.3-0.5. Crow and Champagne obtained a Strouhal number of 0.45 for the preferred mode of the free jet, but they observed, by exciting the jet with an actuator placed upstream, that the jet maximum sensitivity was for St =0.3. Many studies were carried out on these excited jets, but generally with discrete frequencies of excitation, and by remaining close to the frequency corresponding to the preferred mode, or harmonics of this one [START_REF] Morrison | Axial wavenumber measurements in axisymetric jets[END_REF], [START_REF] Petersen | On the preferred mode of jet instability[END_REF], [START_REF] Hussain | The 'preferred mode' of the axisymmetric jet[END_REF], [START_REF] Chao | Visualisation and image processing of an acoustically excited jet flow[END_REF], [START_REF] Cho | Vortex pairing in an axisymmetric jet using two-frequency acousting forcing at low to moderate strouhal numbers[END_REF], [START_REF] Cerecedo | Changes in a coflowing jet structure caused by acoustic forcing[END_REF]). In order to better understand the dynamics of the excited jets at low Reynolds numbers, and to better define the response domain of the jets submitted to longitudinal acoustic excitations, [START_REF] Birbaud | Dynamics of free jets submitted to upstream acoustic modulations[END_REF] undertook a systematic study by varying the excitation frequencies in a large range. Reynolds numbers were 2000 and 4000. They varied the Strouhal number between 0.18 and 2.64. They highlighted that the jets strongly respond for Strouhal numbers close to 0.4. With this condition, the jets exhibit very strong and very reproducible coherent structures.

The vorticity can easily reach 800 s -1 inside the vortex core for velocity fluctuations of about 10 % of the mean velocity at the nozzle exit. An example of tomographic visualization is presented in Figure 1. For Strouhal numbers of about 0.2, the initial oscillation changes quickly to reveal the second harmonic of the main frequency, while leaving the fundamental one visible (first harmonic) at a distance of a few diameters downstream from the nozzle outlet. For a Strouhal number of 0.8, they observed pairing between vortices and the creation of a subharmonic in the jet dynamics. For Strouhal numbers higher than 2, the jet very quickly looses the sensitivity to the imposed excitation and it becomes unstable, revealing the preferred mode corresponding to St = 0.4. According to these results, it is thus possible to define the excitation frequency range for which very robust and very repetitive coherent structures can be obtained at low Reynolds numbers. For Strouhal numbers between 0.3 and 0.65, the structures are conserved for distances of some diameters downstream from the nozzle. From a situation of excited reactive jet laden by fuel droplets with a mean Sauter diameter D 32 of 80 µm and a mean diameter D 10 of 20 µm, [START_REF] Pichard | Effects of acoustic modulations on a spray flame[END_REF] studied the dynamics of a laminar spray flame submitted to an acoustic excitation. Their results show that the spray flames are extremely sensitive to low Strouhal numbers and that they present little variations of the heat release rate for Strouhal numbers higher than 1. This result is, of course, conditioned by the size of the flame, but it indicates that the relevant Strouhal numbers to modify the spray flames will be lower than one. This result is S4_Tue_C_23 ICMF 2007, Leipzig, Germany, July 9 -13, 2007 corroborated by the study of [START_REF] Wark | The structure of an acoustically forced reacting two-phase jet[END_REF], who studied a methanol spray flame, burning mainly in diffusion mode. The mean diameter D 10 of the methanol droplets was 15 µm and the D 32 diameter was of 48 µm at the jet exit. They had a Reynolds number of 5000 and they used a Strouhal number of 0.49 to strongly perturb the reactive jet. They observed that very strong structures were imposed to the flame, corresponding to the beginning of swirling rollings up of the jet. They highlighted clusters of droplets corresponding to the survival of the largest droplets to the crossing of the reactive front. They thus observed an increase in the mean droplet diameter after the flame and a variable concentration on the jet axis. Their study shows that the droplet distribution strongly undergoes the effect of the swirling action of the vortices and that the droplets have a dynamical behaviour, which is not that of the gas flow but obviously conditioned by the stretch rate, the vorticity, the droplets drag and the evaporation rate. Thus this excited multiphase jet seems to be an interesting situation for a specific study on the behavior of a spray in a flow exhibiting zones of high vorticity and for the validation of numerical codes. The experimental configuration retained in this study is similar to the one of [START_REF] Pichard | Effects of acoustic modulations on a spray flame[END_REF], in a non-reactive laden jet situation. Reynolds numbers selected will be of the order of 1000, in order to have a laminar gas flow at the nozzle outlet. Several liquids are usable in the configuration chosen, but results presented here are obtained with decane, in order to have a reduced evaporation rate in a first step. Other liquid fuels like n-heptane, methanol or ethanol, which have a higher evaporation rate, could also be used. In order to allow easy comparisons between numerical simulations and experimental results, values will often be presented in a non-dimensional form. The reference time is the acoustic time corresponding to the frequency imposed to pulse the jet :

6 th International Conference on Multiphase Flow,
τ f = 1 f . ( 2
)
The reference velocity U 0 is the mean axial velocity at the nozzle exit. A reference length can be deduced from these two quantities, which corresponds to a convective time. This length is characteristic of the coherent structures displacement [START_REF] Birbaud | Dynamics of free jets modulated by plane acoustic waves : Part 1. experiments[END_REF]). Another size used for the nondimensionning is the characteristic size of the droplets. In the present study, this dimension will be taken equal to 60 µm, which corresponds to the diameter of the largest droplets in the spray at the burner outlet. The Stokes number, or relaxation factor R, comparing the droplet dynamical time to the acoustic time will be defined as:

R = τ p (S) τ f , (3)
where the droplet dynamical response time reads :

τ p (d p ) = ρ l d 2 p 18 µ g , τ p (S) = 2 ρ l S 9 π µ g (4)
and ρ l is the liquid volumic mass, d p the droplet diameter and µ g the dynamic viscosity of the gas.

Experimental Set-up

An experimental facility has been set up for the study of interactions between acoustic waves and a spray. This set up is presented in figure 2. The injector consists of a converging nozzle, and a cylindrical tube 120 mm long, placed upstream from the nozzle, and containing various grids and honeycombs to produce a laminar flow. The end piece is a cylindrical tube 30 mm long with a diameter of 15 mm. The air jet flowing through this injector is acoustically forced to control the formation and evolution of vortices in the near field of the jet. Two loudspeakers are placed on the sides of the facility to generate longitudinal acoustic waves in the injector cylinder. A spray is obtained from an ultrasonic atomizer, which is fed by a syringe. A specific device is used to control and regulate the liquid flow through the syringe. Different liquids can be used such as heptane, decane, methanol and ethanol. The air flows through the cylinder into the facility to drag the droplets. In order to obtain a homogeneous flow at the nozzle exit, the distance between the atomizer tip and the nozzle exit is quite long (180 mm). This gives a droplet residence time of about 1s [START_REF] Pichard | Effects of acoustic modulations on a spray flame[END_REF]). To benefit from the positive effect of gravity on the droplet formation and motion, the facility is put in the downward direction. In the present study, only decane is used with a flow rate of 26 g.h -1 delivered by the syringe. The air flow corresponds to a mean bulk exit velocity U b =0.75 m.s -1 (Re = 715). The mean exit velocity on the axis is U 0 =1.15 m.s -1 as it will be shown later. This operating condition aims at preventing any turbulent fluctuations of the flow. Moreover, these conditions are chosen to yield a high vorticity level of the gaseous flow (up to 600 s -1 ). For that, the velocity modulation at the outlet was set at u ′ 0 =0.21 m.s -1 (rms value). The modulation frequency is f = 50 Hz so that the Strouhal number is St = 1.0 in this study. This Strouhal number corresponds to a high sensitivity of the jet and an important generation of vortices closer to the nozzle exit, than for a lower Strouhal number corresponding to the more amplified instability mode [START_REF] Birbaud | Dynamics of free jets modulated by plane acoustic waves : Part 1. experiments[END_REF]). With this high Strouhal number, more vorticies are generated in the vicinity of the nozzle exit. The S4_Tue_C_23 ICMF 2007, Leipzig, Germany, July 9 -13, 2007 Stokes number has been defined previously and is R=0.49 for these experimental conditions.

6 th International Conference on Multiphase Flow,

Measurement techniques

Particle Image Velocimetry (PIV)

The PIV technique is used to characterize the gas flow and droplet velocities. The set-up consists in a double pulsed Continuum Nd:Yag laser (energies of 50 mJ/pulse or 10 mJ/pulse were used, depending on the kind of measurements). A laser sheet formed by means of spherical convergent and cylindrical divergent lenses is focused on the central plane of the injector. The thickness of the laser sheet is about 500 µm. The particle Mie scattering is collected on a Dantec Hi-Sense CCD camera (1600×1186 pixels 2 -8 bits) equipped with a Nikkor 60 mm lens. An interferential filter centred on 532 nm (minus or plus 5 nm) is placed in front of the lens to filter the background emission. The air flow is seeded with small oil droplets of 2.5 µm mean diameter [START_REF] Durox | Flow seeding with an air nebulizer[END_REF]). The time delay between two successive laser shots is 400 µs. Velocity calculations are performed with an iterative cross-correlation algorithm (Dantec Flow-Manager). An iteratively decreasing interrogation window size is used with final value fixed to 16×32 pixels 2 with an overlap of 50 %.

Interferometry Particle Imaging (IPI)

The IPI technique aims at determining the size of transparent droplets, which flow through the light sheet and scatter light towards a camera. This technique has been developed recently, thanks to the progress in the field of CCD cameras and computer power. IPI measurements are not yet widely used as can be Phase Doppler Velocimetry. This optical diagnostic has been designed under different acronyms, like PMSI (Planar Mie Scattering Interferometry), ILIDS (Interferometric Laser Imaging for Droplet Sizing), PPIA, (Planar Particle Imaging Analysis). The first work that showed the possibility to measure the size of a droplet thanks to interference fringes is due to [START_REF] König | A new light scattering technique to measure the diameter of periodically generated moving droplets[END_REF]. Recent IPI methods are attributed to [START_REF] Glover | Interferometric laser imaging for droplet sizing: a method for droplet size measurement in sparse spray system[END_REF] which used for the first time a light plane and coupled PIV and IPI techniques.

In this paper, we only present the principles of this technique based on the interference between reflected and refracted light rays, travelling through a transparent droplet. The ray path difference between reflected and refracted light rays induces an interference fringe pattern for a coherent light source. The fringes are visualized in the far field by defocusing the diffraction image created on the camera lens. Each defocalised spot, corresponding to a droplet passing through the light sheet, has a circular shape on the image plane (figure 4(b)). Using the geometrical approach [START_REF] Glover | Interferometric laser imaging for droplet sizing: a method for droplet size measurement in sparse spray system[END_REF]), the difference between reflected and first order refracted light ray paths can be expressed. A simple linear relation is found between N f the number of fringes and d p the particle diameter:

N f = κd p (5)
where the geometrical factor κ is expressed by :

κ = α λ     cos φ 2 + m sin φ 2 m 2 -2m cos φ 2 + 1     (6)
λ is the light source wavelength, φ is the observation angle, relatively to the light source direction, m is the relative refraction index of the liquid related to the gas, α is the collection angle and can be expressed as:

α = 2 sin -1 d a 2z l (7)
where d a is the lens aperture diameter and z l the distance between the measurement plane and the camera lens plane. Standard liquid fuels lead to a refraction index between 1.2 and 1.5. For a laser wavelength of about 500 µm, and for a standard lens, the geometric coefficient κ is of the order of 10 5 m -1 . The IPI technique is generally used with an observation angle of 90 deg. For these conditions, it has been demonstrated that κ is roughly independent of the refraction index m for common liquids in air. The present experiment is carried out with two CCD cameras collecting the scattered light at 90 deg relatively to the laser light source (figure 3). The set up is constituted on the base of the PIV set up described in this paper, with two identical Dantec Hi-Sense CCD cameras equipped with a 60 mm Nikkor lens and an interferential filter centered on 532 nm (plus or minus 5 nm). The PIV camera is not only used for the droplet velocity measurement but for the droplet position detection too, as described by [START_REF] Damaschke | Multidimensional particle sizing technique[END_REF]. The second defocused camera is placed at 90 deg relatively to the first camera and collect the scattered light through a semi reflection mirror. The geometrical parameters for the IPI measurements are imposed by the set up configuration. The distance z l between the laser sheet and the camera lens is 320 mm and the camera aperture is f /# = 2.8. The aperture is thus about 21.4 mm in diameter. The corresponding collection angle is α = 3.8 deg and the magnification factor is 27 px.mm -1 . The laser power is strongly decreased to 10 mJ per pulse so that the collected signal does not saturate the CCD sensor of the camera.

To identify a fringe in the interference pattern, the fringe has to be 4 pixel wide at least. The maximum measurable droplet diameter, which can be measured, is thus determined by:

d max = N d i,max 4κ (8)
with N d i,max the number of pixels of the defocused droplet circle diameter. The defocused circle corresponds to the defocused image of a droplet on the image plane. The minimum diameter is defined by at least one fringe in this defocused circle. The minimum diameter, which can be measured is thus:

d min = 1 κ (9)
For the present experimental set up, the droplet diameters are measured in the range [9.3 µm -74.5 µm].

Characteristics of the flow and spray

Mean axial velocities measured at the nozzle exit (z/D=0.25) are presented in figure 5. The axial velocity profile shows a top-hat shape, the boundary layer development is thus limited. The axial velocity fluctuations are weak at the jet center and correspond to less than 1% of the axial velocity. The gas flow at the nozzle exit is laminar.

The diameter D 10 distribution is presented in figure 6 for a non excited jet. This mean diameter is obtained in windows of 3.2 mm 2 and from 40 instantaneous images. As visible in figure 6, the diameter spatial distribution is rather homogeneous in the jet, ranging from 20 to 30 µm. The size distribution at the burner exit is presented in figure 7. It has been checked that this distribution depends neither on time nor on radius. It also remains the same when the jet is submitted to an excitation.

The mean size distribution presented in figure 7 is obtained from 40 instantaneous images in a window of 150 mm 2 corresponding to the nozzle exit. The mean droplet size is D 10 = 26 µm and the Sauter diameter is D 32 = 36 µm. This decane spray is weakly polydispersed, with a weak droplet size dispersion. The diameter distribution is limited to a maximum droplet size of 60 µm. This value is representative of the maximum droplet diameter, given that only a few droplets are measured with diameters bigger than 60 µm. It should be noted that the processing algorithm generates systematic errors for strong overlapping rate. In the region of dense spray, the fringes detection may be confused with the bright circle of the defocused droplets. These errors induce virtual droplet size between 10 and 15 µm in diameter, which have been arbitrarily removed from the size distribution.

Phase locked experiments

The full cycle of oscillations is divided into 20 phases regularly spaced. A set of acquisition is composed of two focused tomographic images for the PIV computation and the particle position detection, and two defocused images for the droplet size measurement. In order to obtain correct statistics, the image acquisition has been phase-locked by synchronizing the laser and the camera with the loudspeaker excitation signal. For each of the 20 phases, 20 phase-locked pairs of images are acquired and treated. Results are then phase-lock averaged. The time evolution of the gas flow velocity at the nozzle outlet, on the central axis, is presented in figure 8. The gas flow velocity has been measured by seeding the air flow with small oil droplets. The sampling rate of the experiment is rather good and the experimental data are able to reproduce the sinusoidal evolution of the velocity. On the same plot, the decane droplet velocity evolution measured at the same point is shown. This velocity has been measured independently but with the same reference time given by the loudspeaker signal. The phase difference measured between gas flow and droplets velocities is probably due to the droplet inertia and the gravity effect which prevent them to perfectly follow the sinusoidal motion of the gas through the convergent nozzle. An amplitude difference is also observed between the gas velocity and the droplet velocity. It has been observed that, for the non excited jet case, the mean axial droplet velocity exceeds the mean axial gas velocity at the nozzle exit by 9%. This velocity difference can be attributed to the gravity effect, inducing larger velocities for heavier particles such as liquid droplets. A free fall velocity difference between the gas velocity (U g ) and the particle velocity (U p ) can be estimated by the product τ p (S) g which is equal to 0.081 m.s -1 for the biggest droplets of 60 µm diameter, with g the gravity.

The average velocity difference weighted by droplet surface which appears when we want to evaluate the mean free fall velocity of the spray thus seems to correspond quite well to the velocity differential measured between the gas and liquid phases.

In figure 8, the gas flow is modulated with a fluctuation of

u ′ 0 / U 0 =18 % (U 0 =1.15 m.s -1
).

An image of the decane spray is presented in figure 9(a) for an acoustic modulation at 50 Hz. This image represents an instantaneous view of the flow. The two phases, gas and liquid can be clearly identified on this image. The gas flow is traced with small oil droplets and can be visualized by light gray dots, while the decane droplets appear as bigger black dots. The maximum vorticity reaches 600 s -1 in the core of the vortex (figure 9(a)) and the maximum axial strain rate is of the order of κ s =50 s -1 at the head of the vortex on the central axis. This image shows that the interaction between the high vorticity gas flow and the spray leads to a non uniform distribution of droplets. The vortex cores present voids of decane droplet concentration whereas high droplet concentration appear at the leading edges (head) of the vortices. Droplets are ejected from the high vorticity core and follow the external vortex arms. In figure 9(b), the corresponding defocused image is presented. Circles of interference fringes appear for each single droplet. At the jet center, near the vortex head, the spray is dense. High droplet concentration zones induce a strong overlapping rate of the defocused circles. For these particular zones of the flow, the IPI validation reaches a maximum value of 40 %. Even if the results are biased due to the overlapping, we have chosen to keep the same experiments setting in order to have a good resolution on the whole domain of interest. Droplet concentration is measured from the tomographic images by considering the number of droplets in a given volume and is presented in figure 10(a). Because of a strong laser reflexion on the nozzle exit, the first measurement . High droplet concentration zones correspond to a combined effect of the gaseous axial velocity and radial velocity but are closely related to a maximum of the axial velocity gradient. At the head of the vortices, the gas flow is strongly decelerated followed by a new acceleration further downstream, inducing a peak of velocity gradient.The droplets slow down in this specific zone of the flow leading to a liquid concentration increase.

The gas velocity divergence is written for a cylindrical system (r, z):

dV z (r, z, t) dz + dV r (r, z, t) dr + V r (r, z, t) r = 0. (10)
For an axisymmetrical jet, at a given time, the radial velocity can be linearized in a neighborhood of the axis of symmetry which should correspond to zero radius, with a coefficient k.

V r (r, z) = k(r, z) r. (11)
The zero divergence condition should thus imposes that:

dV z (r, z) dz = -2 • dV r (r, z) dr . ( 12
)
It can be observed from the experimental measurements that the radial velocity field is effectively linear as a function of radial position in a given zone which should be related to the axis of symmetry of our 2D axisymmetrical configuration.

The radial velocity gradient is represented in figure 10(b). It has been multiplied by minus two, in order to check the zero-divergence condition along the axis. For the first two structures (z/D<1.6) the zero-divergence condition is satisfied, since the two curves are superimposed. For larger axial distances, the flow structure is perturbed by the surrounding air flow fluctuations. In such a case the 2D axisymmetrical assumption of the experimental configuration slightly deteriorates at larger distance from the nozzle exit.

Eulerian model

We point out some characteristics of the Eulerian multi-fluid model and complete details can be found in [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF]; de [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF]. The purpose of the present section is to introduce the multif-fluid model for the description of the dynamics of polydisperse sprays.

Mesoscopic level description: Williams equation

The (2007), from the ideas of [START_REF] Greenberg | An opposed jet quasi-monodisperse spray diffusion flame[END_REF] allows to describe polydisperse evaporating sprays and the associated size-conditioned dynamics. This approach relies on the derivation of a semi-kinetic model from the Williams equation using a moment method for velocity, but keeping the continuous size distribution function [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]).

The evolution of the spray can be described at what can be called a mesoscopic point of view by a simplified Williams transport equation :

∂ t f + ∇ x • (uf ) + ∇ u • (Ff ) = 0,
where ∇ x is the gradient vector in two-dimensional physical space and ∇ u is the gradient vector in two dimensional velocity space. The Stokes drag force F = F(t, x, u, S) is due to the velocity difference with the gaseous phase which depend on the local gas properties and this dependence is written in the (t, x) dependence. Many more phenomena can be described by a Williams equation, such as evaporation, heating of the droplets, coalescence and break-up of droplets as shown in Laurent et al. (2004b,a); [START_REF] Massot | Multiphase Reacting Flows: Modelling and Simulation[END_REF]. However, for the present study, we start with a non-evaporating configuration and the proposed equation is sufficent for our purposes. The drag force is given by a Stokes law :

F(t, x, u, S) = U g (t, x) -u τ p (S) + g with τ p (S) = ρ l S 18πµ g ,
where U g is the gas velocity, µ g , represents its dynamics viscosity and ρ l is the liquid density. The dynamical viscosity µ g and liquid density ρ l are assumed to be constant. A non-dimensional form of the kinetic equation is given by

∂ t * f * + ∇ x * • (u * f * ) + ∂ u * • U * g -u * St(S * ) + g * f * = 0.
Characteristic quantities of the gazeous flow are used in the dimensionless variables: ICMF 2007, Leipzig, Germany, July 9 -13, 2007 and where the non-dimensional gravity reads g * = gτ g /U 0 . The characteristic value of the velocity is taken to be the mean axial velocity U 0 and the characteristic time the acoustic time τ f . The characteristic length L 0 is such that L 0 = τ f U 0 . We take a maximum droplet surface S 0 which was related in the experimental conditions to a droplet diameter of 60µm (see section 1). Since the description of the whole dynamics in the velocity, space and size phase space is not required, we will reduce the dimension of this phase space in the next subsection. From now on, asterisk of dimensionless variables will be omitted.

u * = u U 0 , U g * = U g U 0 , x * = x L 0 , S * = S S 0 , t * = t τ g , S4_Tue_C_23 6 th International Conference on Multiphase Flow,

Semi kinetic model

The key idea is to reduce the size of the phase space and to consider only the moments of order zero and one in the velocity variable at a given time, a given position and for a given droplet size. The obtained conservation equations, called the semi-kinetic model for two fields n(t, x, ν) = f du and u(t, x, S) = f udu/n(t, x, S) are only in a closed form under a precise assumption on the support of the original NDF in the whole phase space: the velocity distribution at a given time t, given location x and for a given droplet size S in a Dirac delta function [START_REF] Laurent | Multi-fluid modeling of laminar poly-dispersed spray flames: origin, assumptions and comparison of the sectional and sampling methods[END_REF]. Assumptions are done on the spray distribution function:

• For a given size, at a given point (t, x), there is only one characteristic average velocity u(t, x, S).

• The velocity dispersion around the average velocity u(t, x, S) is zero in each direction, whatever the point (t, x, S).

Assumptions define the structure of f : f (t, x, S, u) = n(t, x, S)δ(uu(t, x, S)), and the semi-kinetic model is given by two partial differential equations in the variables n(t, x, S) and u(t, x, S) which express respectively the conservation of the number density of droplets and their momentum at a given location x and for a given size S:

∂ t n + ∇ x • (n u) = 0, ∂ t (n u) + ∇ x • (n u ⊗ u) = F.
(13)

Multi-fluid model

We choose a discretization 0 = S 0 < S 1 < ... < S N for a droplet size phase space and average the obtained system of conservation laws over each fixed intervals [S k , S k+1 [. Two assumptions are then introduced:

• In ine section, the characteristic averaged velocity do not depend on the size of the droplets,

• The form of n as a function of S is supposed to be independent of t and x in a given section, thus decoupling the evolution of the mass concerntation of droplets in a section from the repartition in terms of sizes.

These assumtions are equivalent to assume the NDF in size and in velocity inside each section:

∀S ∈ [S k , S k+1 [ , n(t, x, S) = m (k+1) (t, x)κ (k+1) (S), u(t, x, S) = u (k) (t, x).
where m (k) is the mass concentration of droplets in the k th section. The conservation equations for the k th section then read:

∂ t m (k) + ∇ x • m (k) u (k) = 0, ∂ t m (k) u (k) + ∇ x • m (k) u (k) ⊗ u (k) = F
The multi-fluid model's system thus possesses a similar structure as the pressureless gas dynamics.

7 Numerical resolution

General scheme

In our computations, transport and drag processes are treated separately through a Strang splitting algorithm that has the following structure:

• Drag force during ∆t/2,

• Transport during ∆t,

• Drag force during ∆t/2.

This approach has the great advantage to preserve the property of the different schemes we use for the different contributions as for example maximum principle or positivity. This splitting algorithm is second-order accurate in time and allows if necessary highly optimized algorithms de [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF].

With the generic notations ρ for the mass density of one section, u for its radial velocity and v for its axial velocity, the system to solve for the transport is :

               ∂ t ρ + 1 r ∂ r (rρu) + ∂ z (ρv) = 0, ∂ t (ρu) + 1 r ∂ r (rρu 2 ) + ∂ z (ρvu) = 0, ∂ t (ρv) + 1 r ∂ r (rρuv) + ∂ z (ρv 2 ) = 0. ( 14
)
For the 2-D axisymmetric case we are dealing with on structured grids, we can further use a dimensional splitting of the 2D axisymmetrical transport scheme such that we successively resolve :

• Transport in the radial direction of a step ∆r/2,

• Transport in the axial direction of a step ∆z,

• Transport in the radial direction of a step ∆r/2, which has the great advantage to be second order accurate in space and time once the "building blocks" are second ordre accurate in time. The transport in z-direction has already been treated in [START_REF] Massot | Eulerian Multi-fluid method for the numerical simulation of evaporating polydisperse sprays : modelling and numerical issues in multi-dimensional configurations[END_REF]; [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] and the one in r-direction need a special treatment which is detailed in the next section.

The drag force is characterized by a relaxation of the liquid velocity to the gas velocity with a characteristic Stokes number R(S). This leads, at the semi-kinetic level, to a system of equation in time parametrized by the spatial coordinate which can be solved exaclty for a known unstationary gas velocity given by the experiment. ICMF 2007, Leipzig, Germany, July 9 -13, 2007 The obtained scheme for the multi-fluid reads, for a stationary gaseous flow field without gravity :

m (k) (t n+1 ) = m (k) (t n ), u (k) (t n+1 ) = U g + u (k) (t n ) -U g exp -∆t St(S k u )
.

This numerical scheme ensures a good accuracy in time but it is only first-order accurate in droplet surface because of the use of a constant velocity in each section. In the case of an unsteady gaseous flow field and with the gravity field, a θ-scheme has to be used.

Transport scheme for the radial direction

The radial part of the system ( 14) reduces to the following system :

               ∂ t ρ + 1 r ∂ r (rρu) = 0, ∂ t (ρu) + 1 r ∂ r (rρu 2 ) = 0, ∂ t (ρv) + 1 r ∂ r (rρuv) = 0. (15)
The first two equations constitute the pressureless gas equations for rρ and u, as described in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], except that r is positive, and the last equation is just a transport equation for v at the velocity u.

The same technics as in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] can then be used in order to derive a second order scheme which preserves the maximum principle on the velocities and the positivity of ρ. It is a finite volume scheme based on the equivalence between a macroscopic and a microscopic level of description for the pressureless gas equations, that is to say between the system (15) and the following kinetic equation :

∂ t f + ξ r ∂ r (r f ) = 0, (t, r, ξ, ζ) ∈ R + ×R + ×R×R (16)
with :

f (t, x, ξ, ζ) = ρ(t, x)δ(ξ -u(t, r))δ(ζ -v(t, r)). ( 17
)
The values of ρ, u and v are then recovered from f by the formula :

  ρ u v   (t, x) = R 2   1 ξ ζ   f (t, x, ξ, ζ)dξdζ. (18)
One defines the discretization 0 = r -1 2 < r 1 2 < . . . < r i+ 1 2 < ... over R + and sets

∆r i = r i+ 1 2 -r i-1 2 ∀ i ∈ {0, 1, ...}.
The scheme is a finite volume method giving approximations ρ n j , u n j and v n j of averaged values on each cell [r j-1/2 , r j+1/2 ] of ρ, u and v at each discrete time t n :

  ρ n j ρ n j u n j ρ n j v n j   ≃ 1 r i ∆r i r i+ 1 2 r i-1 2 r   ρ(t n , r) ρ(t n , r)u(t n , r) ρ(t n , r)v(t n , r)   dr. (19)
The principle of the derivation of the scheme is the following : first, at time t = t n , a distribution function f n (r, ξ, ζ) is reconstructed from the averaged values ρ n j , u n j and v n j . This comes from Eq. (17) and, for example, a piecewise linear reconstruction of ρ(t n , r), u(t n , r) and v(t n , r) with adequate slope limiters. Second, the kinetic equation is solved analytically between t n and t n+1 :

f (t, r, ξ, ζ) = f n (r -ξ(t -t n ), ξ, ζ). (20)
Finally, a projection of f (t n+1-, r, ξ, ζ) is done to find ρ n+1 j , u n+1 j and v n+1 j , which corresponds to the average on each cell of (18) at t = t n+1-. This leads to the following scheme :

  ρ n+1 i ρ n+1 i u n+1 i ρ n+1 i v n+1 i   =   ρ n i ρ n i u n i ρ n i v n i   - ∆t r i ∆r i F i+ 1 2 -F i-1 2 , (21)
with the fluxes

F i+ 1 2 = r i+ 1 2 ∆t tn+1 tn R 2   1 ξ ζ   ξf (t, r i+ 1 2 , ξ)dξ dζ dt. (22) Note that F -1 2 = 0.
The obtained fluxes rely, through (17), on the reconstructions of ρ n (r) = ρ(t n , r), u n (r) = u(t n , r) and v n (r) = v(t n , r) from the discrete values ρ n j and u n j . Different type of reconstructions are proposed in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. We choose the one which gives good results without being too complex : a piecewise linear reconstruction. The functions ρ n (r), v n (r) and v n (r) are then written, for r between r j-1/2 and r j+1/2 :

         ρ n (r) = ρ i + Dρ i (r -r i ), u n (r) = u i + Du i (r -r i ), v n (r) = v i + Dv i (r -r i ), , (23)
with

ρ i = ρ n i -Dρ i ∆r 2 i 12r i , (24)
u i = u n i - ∆r 2 i 12r i Du i 1 + Dρ i ρ n i r i - ∆r 2 i 12r i , (25)
v i = v n i - ∆r 2 i 12r i Dv i 1 + Dρ i ρ n i r i - ∆r 2 i 12r i , (26)
in order to have the conservation property on the momentum, i.e. in order that the property ( 19) is exact for the reconstructed functions.

The slopes Dρ n j , Du n j and Dv n j are chosen (like in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] for the velocity) in such a way that the maximum principle property is preserved. We then choose, for i > 0 :

Dρ i = 1 2 [sgn(ρ n i+1 -ρ n i ) + sgn(ρ n i -ρ n i-1 )] × min    |ρ n i+1 -ρ n i | ∆r i 1 -∆ri 6ri , |ρ n i -ρ n i-1 | ∆r i 1 + ∆ri 6ri    , (27)
Du i = 1 2 sgn(u n i+1 -u n i ) + sgn(u n i -u n i-1 ) × min |u n i+1 -u n i | ∆r i (1 -β i ) , |u n i -u n i-1 | ∆r i (1 + β i ) , (28)
with

β i = ∆r i 6r i 1 + Dρ i ρ n i r i - ∆r 2 i 12r i (29)
and

Dv i = 1 2 sgn(v n i+1 -v n i ) + sgn(v n i -v n i-1 ) × min |v n i+1 -v n i | ∆r i (1 -β i ) , |v n i -v n i-1 | ∆r i (1 + β i ) , (30)
For i = 0, the axial symmetry (at r = 0) imposes that Dρ 0 = 0, Dv 0 = 0 and the formulas (28) can be used for Du 0 if we define u n -1 = -u n 0 . Let us denotes ρ I j , ρ I I j , u I j , u I I j , v I j , v I I j the corresponding values of ρ n (r), u n (r) and v n (r) at the bounds r j-1/2 and r j+1/2 of the j th cell. The flux defined by ( 22) are then calculated, like in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF], using the kinetic solution (20), the shape (17) of this function and the previously defined reconstruction. Under the CFL condition

∆t sup|u n (r)| < ∆r, (31)
and the piecewise nonovertaking condition ∆tDu n j > -1, we get

F i+ 1 2 = F + i+ 1 2 + F - i+ 1 2
, with

F + i+ 1 2 = ρ I I i α I I i ×    r i+ 1 2 -∆t 2 α I I i α I I i r i+ 1 2 1 + ∆t 2 Du i -∆t 2 α I I i 1 + ∆t 3 Du i r i+ 1 2 v I I i -∆t 2 α I I i Dv i -∆t 2 α I I i v I I i -2∆t 3 α I I i Dv i    - ∆t 6 Dρ i (α I I i ) 2 ×    3r i+ 1 2 -2∆tα I I i α I I i r i+ 1 2 (3 + ∆tDu i ) -2∆tα I I i 1 + ∆t 4 Du i r i+ 1 2 3v I I i -2∆tα I I i Dv i -2∆tα I I i v I I i -3∆t 4 α I I i Dv i    with α I I i = max{u I I i , 0} 1 + ∆tDu i and
F - i-1 2 = -ρ I i α I i ×    r i-1 2 + ∆t 2 α I i -α I i r i-1 2 1 + ∆t 2 Du i + ∆t 2 α I i 1 + ∆t 3 Du i r i-1 2 v I i + ∆t 2 α I i Dv i + ∆t 2 α I i v I i + 2∆t 3 α I i Dv i    - ∆t 6 Dρ i (α I i ) 2 ×    3r i-1 2 + 2∆tα I i -α I i r i-1 2 3 + ∆t 2 Du i + 2∆tα I i 1 + ∆t 4 Du i r i-1 2 3v I i + 2∆tα I i Dv i + 2∆tα I i v I i + 3∆t 4 α I i Dv i    with α I i = -min{u I i , 0} 1 + ∆tDu i .

One-dimensional validation case of radial transport

We conducted different numerical tests to validate the numerical scheme we have derived. A Rieman problem at a point r > 0 has for example be tested inducing the development of a delta-shock that the numerical scheme capture with the correct propagation speed. Since it is not very different to what is done in [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF] (replacing ρ by rρ), this test is not presented here. We merely focus on an other test case which emphasizes the particular properties of our system. Moreover, we presents only results for the ρ and u since, for v, it is just a transport equation at the velocity u and it does not lead to any difficulty. The computational domain is the one-dimensional segment [0,1] discretized with 0 = r 0 < r 1 < ... < r N = 1 such that r j+1 -r j = ∆r for j = 1, ..., N -1 and

Numerical set up

The numerical computation requires experimental gas velocity fields and initial droplet size and velocity. The experimental velocity fields are used to determine the unsteady gas phase velocity in the computed domain and are interpolated on the computational grid. A phase-locked mean gas velocity field measured by PIV, as described previously in the experimental section, is obtained from 40 instantaneous fields and for each of the 20 time steps chosen in the experiments in order to sample one phase. The velocity map corresponds to one half of the jet domain and cover part of the external region of jet. In order to measure velocities of the air flow surrounding the jet, small oil droplets are injected through a ring-shaped perforated pipe, placed around the nozzle. Low injection velocities and homogeneous seeding are used, in order to avoid any perturbation of the central jet. Velocity calculations are performed with an iterative cross-correlation algorithm (Dantec FlowManager) and interrogation window size with final value fixed to 16×32 pixels 2 with an overlap of 50 %.

For the numerical simulations, the 20 phase-locked velocity fields, after being interpolated in space on the computational grid, are interpolated in time during the computation from the 20 phases available at each needed time. The initial conditions are defined at a domain boundary roughly corresponding to the nozzle exit (z/D=0.2mm). For each time step and radial position, initial conditions are defined by a gas velocity (as described previously), a droplet velocity and a droplet density value. Droplet exit velocities are measured experimentally by PIV and the 20 phase-locked mean exit velocity profiles are used as initial conditions for the numerical computation. The droplet density is measured at the nozzle exit by IPI (see experimental section). The size distribution presented in figure 7 is used to define the droplet density for the computation initialization. For each droplet size section (3 µm in diameter) a droplet density (number of droplets in a given volume) is determined and imposed as initial conditons. The same droplet density is imposed for each radial position corresponding to the nozzle width as it has been observed experimentally that the size distribution depends neither on time nor on radius. The numerical discretization used in our simulation is a 200 × 200 domain thus offering sufficient resolution for the purpose of the comparisons with experimental measurements.

Comparison between experimental and numerical results

The experimental mean diameter D 10 map is presented in figure 12(a). The gas flow velocities are superimposed on the mean diameter map. To ease the analysis, only one half of the velocity field is presented in this figure. The vortex propagation velocity, corresponding roughly to U 0 /2, has been substracted to axial gas flow velocities, so that the gas flow velocity is presented in the vortex system of reference. Voids of droplets correspond to the core of the vorticies, and the spray is confined between these vortices. The droplet diameter D 10 ranges from 20 µm to 40 µm in the excited spray, the size dispersion is weak and is similar to the case of the nonexcitated spray (figure 6). Indeed, the size dispersion of the injected spray is weak as seen on the droplet number density function at the injector exit, presented in figure 7. Moreover, the decane evaporation rate is low (less than 0.01 mm 3 .s -1 ).

The mean diameter distribution is not homogeneous in the excited flow. Biggest droplets (30-40 µm dia.) are present along the central axis and in the leading edge of the vorticies where there is a high droplet concentration. In the arms of the vortices, smaller droplet diameters are measured (20-30 µm dia.). Numerical results for the mean diameter D 10 are presented in figure 12(b). Numerical results are able to reproduce the experimental observation, since the diameter spatial distribution is rather homogeneous in both cases and D 10 ranges between 20 µm and 40 µm. Both results present bigger droplets around the core of the vorticies. Arithmetic mean droplet velocities are deduced from a cross correlation algorithm applied to the tomographic focused spray images. The interrogation window size is fixed to 32×32 pixels 2 with an overlap of 50 %. Averaged ve- inside the vortex cores and because of a poor droplet density in the vortex arms. Indeed, the droplets are injected in the central part of the jet and follow the inner side of vortex as shown in figure 9(b). The droplet radial velocities present successive zones of positive and negative velocities, corresponding to the roll up of the vorticies. The successive zones of positive and negative droplet velocities are localized for similar axial distance in the numerical velocity field and in the experimental velocity field. Numerical results are thus able to reproduce the velocity motion of the droplets.

Concerning the numerical mean droplet concentration presented in Figure 15-top, high concentration zones are located at the periphery of the vortices. Whereas experimentally, the concentration shown in figure 15-bottom, is higher in the head of each vortex. Some differences between experimental and numerical results still remain. However, the model seems to well reproduce the mean concentration related to

Conclusions

In the present paper, we have considered the configuration of acoustically pulsated free jets with a polydisperse spray injection in a 2D axisymmetrical pulsated jet. The large vortical structure created by the pulsation leads to well controlled configuration for phase-locked experimental measurements which are representative of the flow in more complex configurations.

In this context, we have provided both a series of detailed experimental measurements through the coupling of laser diagnostics, as well as comprehensive numerical simulations of the polydisperse spray in the gaseous flow field extracted from the experiments. Moreover, we have provided the necessary numerical ingredients in order to be able to simulate such a 2D axisymmetrical configuration with the potential associated singularities. Finally we have compared the numerical simulations and the experimental measurements for four averaged quantities of the spray : mean size, axial and radial velocities, as well as droplet number density. This yields very good comparisons which can be improved following two directions. First, we need an improve level of statistical information in the experimental measurements in order to tackle the characterization of size-conditioned dynamics and eventually vary the structure of the vortical velocity fields at the exit of the nozzle by changing the way the flow is pulsed. Second, we will focus on another fuel for which the evaporation characteristic time is much lower in order to be able to characterize the dynamics of polydisperse evaporating sprays.

Figure 1 :

 1 Figure 1: Tomographic images of an excited jet. Flow exit velocity on the axis: 2.5 m.s -1 , diameter: 22 mm, frequency: 46 Hz. St = 0.41 and Re = 3660. Each image is shifted of π/2 (Birbaud et al. (2007))

Figure 2 :

 2 Figure 2: Lateral view of the experimental setup.

Figure 3 :

 3 Figure 3: PIV and IPI experimental set-up.

  (a) Focused tomographic image. (b) Defocused image.

Figure 4 :

 4 Figure 4: Image pair for IPI processing showing the same zone of the spray. The gray scale has been reversed.

S4_Tue_C_236

  th International Conference on Multiphase Flow,ICMF 2007, Leipzig, Germany, July 9 -13, 2007

Figure 5 :

 5 Figure 5: Gas flow characteristics of the non excited jet. Measurements are made at z/D=0.25.

Figure 6 :

 6 Figure 6: Mean droplet diameter D10 (µm) in the spray without frequency modulation.

S4_Tue_C_236

  th International Conference on Multiphase Flow,ICMF 2007, Leipzig, Germany, July 9 -13, 2007

Figure 7 :

 7 Figure 7: Number density function (droplet number) as a function of the droplet diameter at the injector exit. D10 = 26 µm. D32 = 36 µm.

Figure 8 :Figure 9 :

 89 Figure 8: Exit axial velocity at z/D=0.25 for two complete modulation cycles -circle: gas flow velocity -square: droplets velocity. The line corresponds to a sinusoidal fit for the gas flow velocity

S4_Tue_C_236Figure 10 :

 10 Figure 10: (a) Droplet concentration given in droplet number per mm 3 for the modulated spray. (b) Corresponding axial and radial velocity gradient of the gas flow.

  Eulerian Multi-Fluid model, extended by Laurent and Massot (2001); Massot et al. (2007); de Chaisemartin et al.

S4_Tue_C_23 6 th

 6 International Conference on Multiphase Flow,

S4_Tue_C_23 6 th

 6 International Conference on Multiphase Flow,ICMF 2007, Leipzig, Germany, July 9 -13, 2007

Figure 11 :

 11 Figure 11: Comparison between the numerical (circle), the analytical (square) and the initial (triangle) solution with ∆t/∆r = 1 and Deltar = 0.004. (top) density fields with δ-shock formation at r = 0 (middle) zoom on the density fields witout taking into account the behavior around r = 0 (bottom) velocity field.

S4_Tue_C_236Figure 12 :

 12 Figure 12: Mean droplet size D10 for the modulated spray at phase φ = 0. (a) experimental mean droplet size D10 measured by IPI from 40 instantaneous images and averaged in windows of 3.2 mm 2 . Gas flow velocities (vortex system of reference) superimposed. (b) the numerical solution at time t = 4π on the grid 200 × 200.

Figure 13 :

 13 Figure 13: Mean droplet axial velocities for the phase φ = 0. (a) experimental mean droplet axial velocities obtained by PIV with interrogation windows size of 32×32 pixels 2 with an overlap of 50 %. Averaged velocities are obtained from 40 instantaneous image. (b) the numerical solution at time t = 4π on the grid 200 × 200.

S4_Tue_C_236Figure 14 :

 14 Figure 14: Mean droplet radal velocities for the phase φ = 0. (a) experimental mean droplet radial velocities obtained by PIV with interrogation windows size of 32×32 pixels 2 with an overlap of 50 %. Averaged velocities are obtained from 40 instantaneous image. (b) the numerical solution at time t = 4π on the grid 200 × 200.

Figure 15 :

 15 Figure 15: Mean droplet concentration obtained for the phase φ = 0. (a) experimental mean droplet concentration measured by IPI from 40 instantaneous images and averaged in windows of 3.2 mm 2 . (b) the numerical solution at time t = 4π on the grid 200 × 200.

. This ensures the symmetry around the axis r = 0. We take the initial data to be

The exact solution is obtained using the characteristics for the variables rρ and u. At time t = 0.5, it is given by rρ(t, r) = f (t)δ 0 (r) + µ(t, r) with f (t) = t 0 rρ(0, r)dr = 5/48 and

Concerning the velocity, the exact solution is drawn on figure 11(bottom). In this test, the initial velocity jump to a higher value at r = 0.35, which leads to a vacuum state, followed by a linearly decreasing part, where the mass accumulates and causes the density to increase. Figures 11(bottom) and 11(top-middle) are respectively the comparison of numerical and analytical solution of density and velocity. It can be seen that the scheme gives good results with, however, some numerical diffusion. Moreover, the initial negative velocity near r = 0 leads to an accumulation of mass density at r = 0, which is particular to our case representing the r part of a 2D axi-symetric problem. It is in fact a singularity for ρ but the averaged value is all the same defined, for the first cell :

For ∆r = 0.004, the exact value of this quantity is 52133.48 whereas we obtain 52586.2 showing a good behavior of the scheme.

S4_Tue_C_23 ICMF 2007, Leipzig, Germany, July 9 -13, 2007 11 Acknowledgements

The present research was done thanks to a Young Investigator Award from the French Ministry of Research (New Interfaces of Mathematics -M. Massot, 2003[START_REF] Massot | Multiphase Reacting Flows: Modelling and Simulation[END_REF], the support of an ANR (National Research Agency -France) Young Investigator Award (M. [START_REF] Massot | Multiphase Reacting Flows: Modelling and Simulation[END_REF]Massot, -2009)). The authors also acknowledge the support from IDRIS-CNRS (Institut de Developpement et de Ressources en Informatique Scientifique, Centre National de la Recherche Scientifique) where some of the computations were performed.


Minh Son Phan

Étienne Baudrier

Loïc Mazo

Mohamed Tajine

ESTIMATION OF ANGULAR DIFFERENCE BETWEEN TOMOGRAPHIC PROJECTIONS TAKEN AT UNKNOWN DIRECTIONS IN 3D

Keywords: tomography, unknown direction, angular difference, Euclidean distance, invariant moments

This paper deals with the estimation of angular difference between two tomographic projections belonging to a set of projections taken at unknown directions. The proposed method extends our former work from 2D to 3D. The method is potential for many applications such as projection refinement or projection classification, which are important in the process of tomographic reconstruction. Unlike to common line based angular estimation, the proposed method does not need reference projections. Our method relies on the selection of projection neighbors with local adaptive thresholds, the calculus of the angular difference for neighboring projections by using properties of moments. The accuracy and the robustness of our method are shown on a test database including 50 3D gray-level images at different resolutions and with different levels of noise.

INTRODUCTION

Tomographic reconstruction for unknown direction projections is a well-known problem in various domains such as cryo electron microscopy or medical imaging. During the process of reconstruction, the distance between two projection vectors is measured and the Euclidean distance is often used. For example, in [START_REF] Coifman | Graph laplacian tomography from unknown random projections[END_REF][START_REF] Fang | sLLE: Spherical locally linear embedding with applications to tomography[END_REF], the authors compute the Euclidean distances between the neighboring projections in order to classify the projections and estimate their directions. In [START_REF] Van Heel | Multivariate statistical classification of noisy images (randomly oriented biological macromolecules)[END_REF], after generating an initial reconstruction model, Van Heel computes the Euclidean distance between the projections of the reconstruction model and those acquired from data for the purpose of refining the reconstruction model. The Euclidean distance is also used in [START_REF] Ben Cheikh | A tomographical reconstruction method from unknown direction projections for 2D gray-level images[END_REF] for simultaneously estimating the projection directions and reconstructing the object.

However, the calculus of the Euclidean distance between two projections has a wide variation even if the two projections are close. This may lead to errors when classifying or refining projections. On the other hand, our previous work in 2D presented a method that can estimate the angular difference between two projections by using the properties of mo-ments [START_REF] Phan | Angular difference measure between tomographic projections taken at unknown directions in 2d[END_REF]. The variation of the angular difference estimated from this method is very smaller than the Euclidean distance. This enables us to extend our work to 3D for estimating the angular difference between the neighboring projections. The organisation of the paper is as follows. We first provide some background notions in Sec. 2. Then, the estimation problem of the angular difference is presented in Sec. 3. The performance of our method is carried out in Sec. 4. Finally, we conclude in Sec. 5.

BACKGROUND NOTIONS

Rotation in 3D

We first present the rotation convention used in this paper. The rotation of a vector in 3D is defined as the rotation of this vector around the three axes x, y and z by the three angles ϕ, θ and ψ. It can be represented as the multiplication of the three rotation matrices R x ϕ , R y θ and R z ψ . For example, the product

R ϕ,θ,ψ = R x ϕ R y θ R z ψ (1)
represents a rotation of a vector around the z, y and x axis by the angles ψ, θ and ϕ, respectively. We also note that

R -1 ϕ,θ,ψ = R z -ψ R y -θ R x -ϕ = R -(ψ,θ,ϕ) .

Projection

The projection of an object can be obtained by first rotating the coordinates system using the rotation matrix R ϕ,θ,ψ , then project the object onto the new (xOy) plane. This procedure is similar to the rotation of the object using the rotation matrix R -(ψ,θ,ϕ) and project the object onto the (xOy) plane. We define the projection as follows Definition 1 (Projection). Let f : R 3 → [0, 1] be a measurable function with compact support and let ϕ, θ, ψ ∈ R.

The projection P f (ϕ, θ, ψ) of f in the direction (ϕ, θ, ψ) is defined by

P f (ϕ, θ, ψ)(x, y) = R f (ρ -(ψ,θ,ϕ) (x, y, z))dz, (2)
where

ρ ψ,θ,ϕ (x, y, z) =   R ψ,θ,ϕ   x y z     T .
Note that P f (ϕ, θ, ψ) = P f (ϕ + π, θ + π, ψ + π). Thus we can assume that ϕ, θ, ψ ∈ [-π/2, π/2].

We also present the direction vector of the projection, that is the unit vector in the new z axis direction as follows v P (ϕ, θ) = R ϕ,θ,ψ (0, 0, 1) T = (sin θ, -sin ϕ cos θ, cos ϕ cos θ) T .

(3)

Note that the direction vector v P (ϕ, θ) does not depend on the angle ψ.

Moment of projection

We recall that for a function g : R 2 → R, the moment of g with order (c, d), where c, d ∈ N is given by

µ c,d (g) = R R
x c y d g(x, y) dx dy.

For a function f :

R 3 → R and ϕ, θ, ψ ∈ R, we note µ f,c,d (ϕ, θ, ψ) = µ c,d (P f (ϕ, θ, ψ)) a moment with order (c, d) of the projection P f (ϕ, θ, ψ).
The moments of projection in 2D are trigonometric polynomials as a function of projection direction as shown in [START_REF] Natterer | The Mathematics of Computerized Tomography[END_REF]. In 3D, we obtain the following result Proposition 1 (Trigonometric polynomial of moment). Let f : R 3 → [0, 1] be a measurable function with compact support, let ϕ, θ, ψ ∈ R and c, d ∈ N. The moment of the projection P f (ϕ, θ, ψ) of the function f , with the order (c, d), in the direction (ϕ, θ, ψ) is a trigonometric polynomial as a function of projection direction.

In the following, for simplicity, we write P and µ c,d instead of P f and µ f,c,d .

ESTIMATION OF THE ANGULAR DIFFERENCE

Problem of estimating the angular difference

We consider a function f whose centroid is at (0, 0, 0), a set of unknown directions

Θ = {(ϕ i , θ i , ψ i ) | 1 ≤ i ≤ n} ⊆ [-π/2, π/2[ 3 , n ∈ N. We set Π = {P(ϕ i , θ i , ψ i ) | 1 ≤ i ≤ n} a set of projections associated to Θ and V = {v P (ϕ i , θ i ) | 1 ≤ i ≤ n} a set of direction vectors of Π.
Our purpose is to estimate the angular difference between two neighboring projections from the set of projections Π. The following formula calculates the angular difference between two projections P and P , that is the inner product between two direction vectors v P and v P d ang (P, P ) = arccos (v P . v P ) .

(4)

In order to calculate the angular difference from the equation (4), the direction vector must be calculated. There are two well-known approaches to the direction estimation problem. The first approach is based on the common line technique as shown in [START_REF] Van Heel | Angular reconstitution: a posteriori assignment of projection directions for 3D reconstruction[END_REF]. The main idea is that two projections share one line in common which is perpendicular to the direction vectors of the two projections. Thus, using the two common lines between three projections enables us to find the direction vectors of these projections. However, the common line technique is time consuming and may lead to accumulative errors when detecting the false common lines. We therefore employ the second approach that estimates the direction vector of projection by analysing the relationship between the moments of projections and the moments of the 3D function [START_REF] Salzman | A method of general moments for orienting 2D projections of unknown 3D objects[END_REF][START_REF] Goncharov | Determination of mutual orientation of identical particles from their projections by the moments method[END_REF]. The method of moments is fast and does not lead to the accumulation of errors as in the first approach. We have developed from [START_REF] Salzman | A method of general moments for orienting 2D projections of unknown 3D objects[END_REF] the formula of direction vector as follows

v P (ϕ, θ) = v 1 P , v 2 P , v 3 P T , (5)
where :

v 1 P = ± 1 a (µmax -µ 2,0 (ϕ, θ, ψ)) (µmax -µ 0,2 (ϕ, θ, ψ)) 1/2
, with a = (µmax -µ med )(µmax -µ min ),

v 2 P = ± 1 b (µ med -µ 2,0 (ϕ, θ, ψ)) (µ med -µ 0,2 (ϕ, θ, ψ)) 1/2
, with b = -(µmax -µ med )(µ med -µ min )

v 3 P = 1 c (µ min -µ 2,0 (ϕ, θ, ψ)) (µ min -µ 0,2 (ϕ, θ, ψ)) 1/2
, with c = (µmax -µ min )(µ med -µ min ) and µmax = max

(ϕ,θ,ψ)∈R max(µ 2,0 (ϕ, θ, ψ), µ 0,2 (ψ, θ, ϕ)), µ med = min (ϕ,θ,ψ)∈R max(µ 2,0 (ϕ, θ, ψ), µ 0,2 (ψ, θ, ϕ)), µ min = min (ϕ,θ,ψ)∈R min(µ 2,0 (ϕ, θ, ψ), µ 0,2 (ψ, θ, ϕ)).
There are four possible direction vectors in (5) due to the unknown signs. These four direction vectors are on four different octants of the coordinates system. However, we can obtain a unique solution in (5) if we make the assumption that both of the projections are in the same octant. This is the case for most neighboring projection pairs if there is a large number of projections and if the distribution of projections is uniform. Thus, we can forget the problem of signs in (5) by simply setting v P = (|v 1 P |, |v 2 P |, v 3 P ) T . Then the angular difference formula between two neighboring projections is modified as follows d ang (P, P ) = arccos

3 i=1 |v i P v i P | . (6)
The problem of searching the neighbors of each projection is presented in the following section.

Searching for the neighbors of projection

In [START_REF] Phan | Angular difference measure between tomographic projections taken at unknown directions in 2d[END_REF], we have developed a method for finding the neighbors of each projection in 2D . The principle of this method is that, for each moment of a given projection, with a given order, we compute a threshold that allows us to find at least one moment corresponding to a neighbor of the projection. The method can also be applied in 3D. However, the moments are not used directly since the projection can be rotated around its direction vector by an arbitrary angle. This changes the value of moments and leads to errors for neighboring searching. Instead, we use Hu moments [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF], for their invariance to the rotation of projection. These invariant moments can be calculated easily from of the (geometric) moments. Note that the (geometric) moments are trigonometric polynomials (Proposition 1) and it is known that the set of trigonometric polynomials is a ring. Thus, these invariant moments are also trigonometric polynomials.

Hu moments are now applied for selecting the neighbors of a given projection. Since, these moments are nonmonotonic, close moments are not synonymous with close direction vectors of projections. Thus, for each Hu moment M(ϕ, θ) of a given projection, we find an interval M(ϕ, θ) ± ε p (ϕ, θ) in which can be found, with some chosen probability p ∈ [0, 1], at least one moment corresponding to a neighbor of the projection. Assuming that the number of projections is sufficient large and uniformly distributed on demi-sphere. We first find for each angle λ ∈ {ϕ, θ}, an interval λ ± δ λ in which there is at least one angle λ . Since the distribution of projections is uniform, we can obtain the formula of δ λ as follows

δ λ =      π 2 (1 -(1 -p) 1 √ n-1 ), if λ = ϕ 1 (v 2 P ) 2 + (v 3 P ) 2 (1 -(1 -p) 1 √ n-1 ), if λ = θ .
We also note that Hu moments M(ϕ, θ) are trigonometric polynomials. Then, by using Bernstein's inequality [START_REF] Bernstein | Sur l'ordre de la meilleure approximation des fonctions continues par des polynômes de degré donné[END_REF] for trigonometric polynomials, we have

ε p (ϕ, θ) ≈ k (δ ϕ + δ θ ) max (ϕ,θ)∈R |M(ϕ, θ)|, ( 7
)
where k is the order of M(ϕ, θ).

In particular, when M = µ 2,0 + µ 0,2 , a better value of ε p can be obtained as follows

ε p (ϕ, θ) ≈ 2 (ξ ϕ δ ϕ + ξ θ δ θ ), (8)
where

ξϕ = (M(ϕ, θ) -µ med )(µ max + µ med -M(ϕ, θ)) 1/2 , ξ θ = (M(ϕ, θ) -µ max )(µ max + µ med + µ min -M(ϕ, θ)) 1/2 .
In conclusion, ( 7) and ( 8) allow us to find sets of projection neighbors for distinct moment orders. Then, the final result is obtained by intersecting these set of projection neighbors.

Remark 1. In Section 3, the exact values of µ max , µ med and µ min used in ( 5), ( 7) and ( 8) are not known. We can only estimate µ max , µ med and µ min from the finite set of projections Π. However, if the angular difference d ang between two neighboring projections is small in probability, then we can assume that the errors on µ max , µ med and µ min less than (d ang ) 2 and do not modify the formulae obtained above.

RESULTS

We apply our method on a set of 50 phantom images of different resolutions (32 3 , 64 3 , 128 3 and 256 3 voxels) as shown in Fig. 1. We first generate randomly a set of 400 triplets (ϕ, θ, ψ), where ϕ, θ, ψ ∈ [-π/2, π/2[. Then, a set of projections P(ϕ, θ, ψ) is computed for each image. We next find the neighbors of each projection P(ϕ, θ, ψ) by using the Eq. ( 7) and [START_REF] Salzman | A method of general moments for orienting 2D projections of unknown 3D objects[END_REF]. Finally, the angular differences d ang between the neighboring projections are estimated by the Eq. ( 6). The Euclidean distance between these neighboring projections are also computed. Our objective is to measure the dispersion of our estimated angular difference d ang (resp. the Euclidean distance) with respect to the exact angular difference d ang . We use the variance-to-mean ratio (VMR) presented in [START_REF] Phan | Angular difference measure between tomographic projections taken at unknown directions in 2d[END_REF] for measuring the dispersion. The VMR should be lower than 0.5 for a small dispersion. As seen in Table . 1, the VMR of our estimated angular difference d ang is very smaller than 0.5 for all image resolutions, whereas VMR of the Euclidean distance is higher than 1. An example for the dispersion of the Euclidean distance and the d ang is illustrated in Fig. 2. We see that the dispersion of the Euclidean distance (x axis) with respect to d ang (y axis) is greater than the d ang . The robustness of our method against noise is also evaluated. A white noise with a given SNR is added to the projections as seen in Fig. 3. The results of estimating d ang are We further applied our method to a real image TAF7 of 160 3 voxels as shown in Fig. 5. 400 noisy projections are generated from the TAF7 and the Wiener filter is used to denoise the projections. The positive results are obtained at the SNR ≥ 10 dB with the VMR ≤ 0.5.

CONCLUSION

This paper extends our new method to the 3D case. We deal with the problem of estimating the angular difference between two neighboring projections taken at unknown directions. We first find the neighbors of each projection based on local adaptive thresholds for the moments of projections. Then formula of estimating the angular difference between two neighboring projections is given. Our method has been evaluated on a set of 3D images at different resolutions and with different levels of noise. The results demonstrated that our estimated angular difference performs better than the Euclidean distance and is robust at moderate levels of noise. Our method can be applied for many applications such as projection refinement, projection classification, etc. Finding the neighboring projections and estimating the angular differences between them are important for the following step of our work. This allows us to calculate the angular difference between any two projections by constructing a neighborhood graph whose edges connect the neighboring projections, weighted by the estimated angular differences and then using the shortest path algorithm on the graph. This work will be completed in the future.

Fig. 1 :

 1 Fig. 1: Examples from the 3D database

Fig. 2 :

 2 Fig. 2: Dispersion of the Euclidean distance (a) and our estimated angular difference (b), from an image of size 64 3 voxels.

Fig. 3 :

 3 Fig. 3: Examples of noisy projections.

Fig. 4 :

 4 Fig. 4: Mean and standard deviation of the VMR for d ang against white noise, from the images of 256 3 voxels.

Fig. 5 :

 5 Fig. 5: Real 3D image TAF7.

Table 1 :

 1 Mean and standard deviation of the VMR for the Euclidean distance and the d ang , according to the image resolutions.

	Resolution Euclidean distance	d ang
	32 3 voxels	1.03 ± 0.40	0.06 ± 0.03
	64 3 voxels	0.98 ± 0.32	0.05 ± 0.01
	128 3 voxels	1.09 ± 0.34	0.05 ± 0.02
	256 3 voxels	0.83 ± 0.37	0.05 ± 0.03