J Selmi

J P Costes

P Lorong

G Poulachon

P Carras

Spindle position influence on milling operation stability

Keywords: Spindle, Dynamic stiffness, Chatter, Stability, Criteria, Milling

Machine tool is the most used manufacturing means in automotive industry. This kind of manufacturing machine is very expensive and has immediate influence on produced workpiece quality.

During machining process design, a tool machine specification sheet is written for each machining operation depending on workpiece required quality, material hardness and cutting conditions. In this way, in order to optimize investments, new machining center should be able to manufacture the quality set by workpiece specification sheet. The main aim of the work presented in this paper is to define and test a new method for spindle dynamic stiffness comparing.

The used approach is firstly based on dynamic criteria attained through FRF (Frequency Response Function) measured experimentally on machining center spindle. Secondly, stability limits are calculated and then obtained experimentally. In this case the criterion of steady state limit is used to compare different spindle configurations regarding the tool.

Introduction:

Besides of machining-center power and torque limits, self-excited vibrations occurrence (chatter) constitutes a real productivity limitation. Otherwise, Peklenik and Gartner [START_REF] Peklenik | workpiece accuracy criterion for the dynamic machine tool acceptance test[END_REF] and Pruvot [START_REF] Pruvot | Conception et calcul des machines-outils[END_REF] have shown that stiffness lack in machining systems could have a negative effect on dimensional produced work piece accuracy.

Lobe stability diagram is a chart which represents stability and instability domains during machining and allows the choice of preferential cutting condition in order to increase productivity and avoid vibration. Tlusty [START_REF] Tlusty | Self-excited vibrations on machine tools[END_REF] and Tobias [START_REF] Tobias | Machine-Tool Vibration[END_REF] have invented the theory of regeneration; they demonstrated that the instability in milling is due to this phenomenon and appears when a depth of cut threshold is reached. They also pointed the theory of stability lobes and demonstrated that areas allowing big material removal rate are situated in high speed region. Usually working in those regions, especially for soft materials machining, activates mode related to the tool. However, hard material machining is associated to structural low frequency mode of machining center subset [START_REF] Ezugwu | Key improvements in the machining of difficult-to-cut aerospace super alloys[END_REF].

Many other studies have been carried on spindle characterization, Gurney [START_REF] Gurney | A graphical method for the determination of the dynamic stability of machine tools[END_REF], Sadek [START_REF] Sadek | Comparative dynamic acceptance tests for machine tools applied to horizontal milling machines[END_REF] and Yuce [START_REF] Yuce | Pulse excitation technique for determining frequency response of machine tools using an-line minicomputer and a non-contacting electromagnetic exciter[END_REF] have used the coefficient of Merit to compare spindles dynamic behavior by setting the structure into forced vibrations. In addition, Tobias [START_REF] Tobias | Dynamic acceptance tests for machine tools[END_REF] used stability chart as a comparative dynamic test certificate. Altintas [START_REF] Altintas | Chatter stability of metal cutting and Grinding[END_REF] has set a model for stability chart prediction. However machining center reliability from avoiding chatter point of view depends on how machining center will be used.

Otherwise, Abele [START_REF]Creating stability lobe diagrams during milling[END_REF], and Pruvot [START_REF] Pruvot | Les broches étude cinématique et statique[END_REF] showed that spindle stiffness depends on its design, especially geometrical dimension, and ball bearing specifications and positions. This paper presents a general approach for machining center qualification based on instability behavior. In the first section, a global approach of spindle qualification is shown using experimental and mathematical methods. In the second part some milling tests are made in order to confirm the proposed method.

Hammer tests:

Experimental method:

The aim of this section is to study the influence of spindle position in a three dimension space on its dynamical behavior in order to validate a method for spindle stiffness evaluation. All tests presented in this paper were made on the same machining center designed by DMG. The machine center used is three axes architecture, equipped with a spindle (KESSLER) and ball screw axis technology. The experimental setups are depicted in FIG. 1. For each axis movement, spindle was positioned at the extreme positions, and then the spindle nose FRF measurements were made in each one. However, the spindle shaft is not easily reachable to make hammer test because of spindle headstock protection. In order to get FRF on spindle shaft nose, the spindle was equipped with a tool holder, which was excited by hammer. Hammer tests were made following instructions and details described below:

-A low mass bandwidth accelerometer [Sensitivity: 100mV/g] was located on the tool holder and connected to computer through a high frequency data acquisition card.

-A modal hammer [Sensitivity: 2,27mV/N] was connected to the acquisition card.

-Impulse hammer shocks are conducted on the tool holder using soft and hard tips in order to enhance signal quality and get a larger signal representation for high and low frequencies.

Results and interpretation:

The different tests were made for all positions described FIG. 2. Spindle positions on (X, Y) map have no influence on spindle nose stiffness. The spindle Y axis positions FRF shown in FIG. 2-a) are approximately superposed.

However, FRFs show that the spindle stiffness is affected by its position on Z axis. The spindle modes amplitude are lower when the spindle is in the highest position, which means that spindle has a better dynamic behavior in this position. Otherwise, when the spindle is set in the lowest position, flexion modes amplitudes are often higher and could affect the machining process quality. Then the spindle is more flexible in this position.

In order to show the difference between the spindle nose FRF and the spindle nose equipped by a rigid tool, one complimentary hammer test was then made in the tool tip following the same experimental previous description. Results depicted in FIG. 3 shows that tool tip FRF magnitude is approximately three times higher than spindle FRF magnitude. As expected, the spindle with tool is more flexible than the spindle nose itself [START_REF] Koppka | Thesis: A contribution to maximization of productivity and workpiece quality of the reaming process by analyzing its static and dynamic behavior[END_REF]. In addition, it shows that mode frequencies are the same for spindle noise itself and spindle nose 3 Machining tests:

Preliminary calculations:

This part of study will focus on spindle positions on Z axis. In order to compare two spindles or two machining center spindle positions regarding the tool, the stability axial depth limit was the chosen criterion. This limit "ap lim " is defined as the maximal axial depth of cut where chatter does not occur for any chosen spindle rotation speed [FIG. 4].

In order to compare the two spindle configurations regarding the tool, an SDOF approximation was made, considering that the system dynamic behavior is led by the highest mode on FRF. A simplified chatter limit prediction model invented by Tlusty [START_REF] Tlusty | Self-excited vibrations on machine tools[END_REF] was used to compare both tested spindle position ability to remove materials without chattering "ap lim ". This is one dimensional model which considers that tool is excited by one radial periodic force regarding teeth in contact with the work piece and gives an approximation of the limit of axial depth of cut:

ap lim ( ) ( )

Equation 1

Where Kt (N/mm²) is the specific cutting force, Re Stability lobe is a chart which makes distinction between two domains steady machining area and instable machining one. Around high frequencies when machining soft material like aluminum, stability lobes contain areas that allow increasing the depth of cut. However, when working around low frequencies, high productivity areas disappear. The stability limit becomes somehow horizontal line [FIG. 4] which chatter could not occur behind for any rotation speed [START_REF] Tobias | Dynamic acceptance tests for machine tools[END_REF], this line is often set by the mode with maximum amplitude. Otherwise, the spindle/tool system dynamic behavior depends largely on the tool flexibility. The spindle dynamic behavior could be completely hidden and dominated by the tool dynamic behavior when the later is flexible [START_REF] Gagnol | Thesis: Modélisation du comportement dynamique des électro-broches UGV[END_REF]. The experimental method described below is used to find a correlation with tapping test and validate spindle stiffness qualification criterion set above. In this experimental test, the used tool has been chosen very rigid "milling tool with 9 carbide inserts" in order to maintain spindle dynamic behavior during machining. Cast Iron workpieces especially designed to be very rigid when clamped were used for this machining trial. Moreover machining was made for the two defined spindle positions on z axis. Experimental setup is depicted in FIG. 5-a) the work piece was mounted in a low position work holding device. The spindle goes up to make machining operation on it. In FIG. 5-b) a work holding device raise the workpiece up 500 mm above the table reference plan, then the spindle remain in high position for face milling. Acceleration and sound acquisition are made during machining operation as described below:

-A low mass, wide bandwidth accelerometer was attached on spindle headstock perpendicular to the feed direction -The acceleration signal was recorded using Daysilab software at 5 kHz sampling frequency. All trials are carried on following the algorithm depicted FIG. 6, using cutting conditions displayed in Table 2 The first cutting conditions are Vc 0 = 100 m/min and ap 0 = 2 mm, the depth of cut is incremented gradually until unstable machining state is reached. Then cutting speed could be incremented, and the same cycle is reproduced many times.

Results and interpretation:

From the transfer function of both the spindle in high and low position on Z axis, the results depicted on the stability chart [FIG. 7] were expected. Indeed machining stability and quality is better when the spindle is less flexible i.e. spindle in high position. Moreover the stability limits as defined in section 3.1, is a straight horizontal line passing by 2.5 mm. Besides of acoustic emissions which help on instable machining detection, the surface quality was also an important indicator which was always good after steady milling. However, when chatter occurs it was most of time inacceptable [FIG. 8]; tool vibration print was left on it and was detected by seeing and touching. Otherwise, in the second machining position: spindle situated in the higher position, the system was very rigid; machining was steady all the time and surface quality was very good. For all tested machining, chatter did not occur. The maximum tested and reached axial depth of cut without chattering was 5 mm. Moreover, any strong acoustic emission or acceleration was detected respectively by the sound or acceleration FFT analysis. Otherwise roughness and surface visual state show a big difference between steady machining and instable ones. When milling is instable, produced surface is deteriorated and roughness level is inacceptable.

Given the fact that chatter did not occur for high spindle position configuration, chatter frequencies analysis concerns only the low machining spindle position. Besides to acceleration analysis, acoustic emmissions was recorded during machining in order to analyze machining frequencies.

FIG.8 -Surface state

The two cases depicted in FIG. 9 show the machining transition from steady milling to chatter. In the FIG. 9-a), the inserts passing frequency (75 Hz) is detected in the FFT signal: before chattering, peaks amplitude levels are very low, then when machining become instable a very high peak is detected around 139 Hz. In the FIG. 9-b), the second analyzed case is about an insert passing frequency of 187 Hz, then just when chatter occurs a very high peak around 139 Hz also appears. The 100 Hz peak appearing in all depicted FFT charts corresponds to the electrical network frequency detected by the accelerometer. Otherwise, for almost, machining chatter occurred around 139 Hz.

Discussion:

FRF chart [FIG. 2-b)] have shown that spindle when located in high position have a better dynamic behavior than when located in low position. Indeed, FRF low spindle position curve is often above on the analyzed frequency domain, which means that the spindle is more flexible. Firstly, this result was confirmed by calculation which has shown that stability limit was higher for the higher spindle position [FIG. 1]. Secondly, it was confirmed by milling tests, since stability chattering limits shown in FIG. 7 have been reached easily after 3mm of axial depth of cut in low spindle position configuration and have not been reached for tests when the spindle was in high position. Otherwise, acceleration values, sound intensity and surface visual and measured states was used like indicators to detect chatter occurrence. Moreover, machining frequencies analyses have shown that high level vibrations amplitude frequencies were around 139 Hz. In FIG. 3 it is obvious that this frequency is very close to the spindle natural frequency at 133 Hz which is detected by a high amplitude peak and a phase change. This frequency corresponds to the minimum FRF real part Re[FRF] min [START_REF] Wang | Suppression of chatter vibration on a CNC machine center-an example[END_REF]. Otherwise, there was no chatter when machining was made in high spindle position configuration. The 139 Hz mode frequency which has the most influence on vibration has disappeared. This fact is demonstrated through FIG. 2-b) which shows that there is no more mode around 139 Hz in this position.

Stability limit Horizontal line is set by the more critical mode (i.e.higher amplitude) in the FRF. This limit was found experimentally equal to 2,5 mm; it was also found by calculation equal to 2 mm, there is a small difference between the two results depending on model reliability and experimental parameters accuracy; however both methods could be used for setting stability limits to compare machining center spindles.

Conclusion:

The study presented in this paper was focused on comparing the dynamic behavior of machining center spindle, and the spindle regarding the tool from avoiding chatter point of view. The proposed criterion "ap lim " was obtained by calculation and by experiment. It was tested on two spindle positions presenting two different behaviors depending on spindle position in the machine workspace. Moreover, tests were carried on in low machining frequencies, which excited the structure modes. This fact was revealed during machining by chatter frequencies detection which was very close to the spindle natural frequency.

Finally, according to the mathematical model calculation presented before and milling results, dynamic machining center behavior could change on workspace. However, stability limit criterion gives a good representation of machining center behavior regarding the tool and could be generated and used in order to compare various spindles in different machining centers. This method is useful to estimate machining center structural limitation by only one hammer test and calculation.

FIG. 2 -FIG. 1 -

 21 FIG.2 -Spindle nose FRF, a) Through Y axis, b) Through Z axis

High

  [FRF] min (mm/N) is the minimum Real part of the tool point FRF and m is the number of teeth cutting simultaneously. Then, results are shown in the table below: Kt= 1200 MPa, m = 4 and Kr = 90°

FIG. 4 -

 4 FIG.4 -Stability Chart

D

  -A microphone is used for sound frequency analyses and helps on vibration frequency detection

FIG. 6 -FIG. 7 -

 67 FIG.6 -Trials algorithm

FIG. 9 -

 9 FIG.9 -Machining frequencies a) Teeth passing frequency = 75 Hz b)Teeth passing frequency = 187 Hz

Table 1 -Stability limits 1.1. Experimental Method:

 1

Table 2 -Cutting conditions

 2


E Masi
email: masi@imft.fr

O Simonin
email: simonin@imft.fr

An algebraic-closure-based moment method for unsteady Eulerian modeling of non-isothermal particle-laden turbulent flows in very dilute regime and high Stokes number

An algebraic-closure-based moment method (ACBMM) is developed for unsteady Eulerian particle simulations coupled with direct numerical simulations (DNS) of non-isothermal fluid turbulent flows, in very dilute regime and for large Stokes numbers. It is based on a conditional statistical approach [1] which provides a local instantaneous characterization of the dynamic of the dispersed phase accounting for the effect of crossing between particle trajectories which occurs for large Stokes numbers.

Introduction

Particle-laden turbulent flows are encountered in several industrial applications and need to be modeled. In order to model the turbulence at industrial scales, Reynolds-averaged Navier-Stokes (RANS) methods or large-eddy simulations (LES) are usually used depending on the nature of the flow. The LES are particularly adequate for modeling unsteady turbulences of internal flows confined in complex geometries as, for instance, in combustion chambers of aero engines. For these two-phase flow applications, the entire mixture needs to be modeled and a question about the choice of the most appropriate method for modeling the particulate phase arises. In the literature, a method which has been extensively used over the years is the Eulerian-Lagrangian LES approach [START_REF] Wang | Large eddy simulation of particle-laden turbulent channel flow[END_REF][START_REF] Boivin | On the prediction of gas -solid flows with two-way coupling using large eddy simulation[END_REF]. According to this approach, the turbulence is resolved by means of the LES and the particle properties are obtained by using a Lagrangian method. So far, the Eulerian-Lagrangian LES approach has presented two major limits: the prohibitive computational cost of the numerical simulation when this method is used for real industrial applications; the difficulty to model the coupling between particles and turbulent fluid flow because of the small non-resolved scales of the turbulence [START_REF] Armenio | FIOROTTO Effect of the subgrid scales on particle motion[END_REF][START_REF] Fede | Numerical study of the subgrid fluid turbulence effects on the statistics of heavy colliding particles[END_REF]. An alternative to the Eulerian-Lagrangian LES approach is the Eulerian-Eulerian LES approach. The latter implies that the particulate phase is described as a continuum medium at the same level as the turbulence so that a straightforward coupling between phases is possible in an Eulerian framework. This approach is considerably less expensive than a Lagrangian technique and thus suitable provided that a satisfactory level of accuracy is ensured. In the literature, several successfully models are available in the framework of the Eulerian-Eulerian RANS modeling [START_REF] Elghobashi | A two-equation turbulence model for twophase flows[END_REF][START_REF] Chen | A turbulence closure model for dilute gas-particle flows[END_REF][START_REF] Reeks | On a kinetic equation for the transport of particles in turbulent flows[END_REF][START_REF] Simonin | Prediction of the dispersed phase turbulence in particle-laden jets[END_REF][START_REF] Zaichik | Modelling of particle dynamics and heat transfer in turbulent flows using equations for first and second moments of velocity and temperature fluctuations[END_REF]. In contrast, the Eulerian-Eulerian LES (or DNS, as direct numerical simulations) of dilute or very dilute highly inertial particle-laden turbulent flows are still a timely topic of research. Indeed, the existing DNS/LES Eulerian models [START_REF] Druzhinin | Direct numerical simulations bubble-laden turbulent flows using the two-fluid formulation[END_REF][START_REF] Ferry | A fast Eulerian method for dispersed two-phase flow[END_REF][START_REF] Rani | Evaluation of the equilibrium Eulerian approach for the evolution of particle concentration in isotropic turbulence[END_REF][START_REF] Ferry | Equilibrium Eulerian approach for predicting the thermal field of a dispersion of small particles[END_REF][START_REF] Shotorban | A Eulerian model for large-eddy simulation of concentration of particles with small stokes numbers[END_REF][START_REF] Ferrante | On the accuracy of the two-fluid formulation in direct numerical simulation of bubble-laden turbulent boundary layers[END_REF] have proven to be able to correctly predict the dispersed phase behaviour only for Stokes numbers smaller than the critical one (St K << 1 based on the Kolmogorov timescale). In this work an appropriate model for predicting unsteady Eulerian particle simulations coupled with DNS of non-isothermal fluid turbulent flows in very dilute regime and for large Stokes numbers is proposed. It should be considered as the starting point for developing the Eulerian-Eulerian LES approach [START_REF] Moreau | Development of gas-particle Euler-Euler LES approach: A priori analysis of particle sub-grid models in homogeneous isotropic turbulence[END_REF].

Unsteady Eulerian Modeling

The algebraic-closure-based moment method (ACBMM) is based on a conditional statistical approach [START_REF] Évrier | Partitioning of particle velocities in gas-solid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] which provides a local instantaneous characterization of the dynamic of the dispersed phase accounting for the effect of crossing between particle trajectories which occurs for large Stokes numbers. This phenomenon has a great impact on the thermal characterization of the dispersed phase as well. Indeed, the occurrence of crossing trajectories entails that particles convey information of their interactions with very distant, and independent, turbulent eddies, that is with different dynamic and thermal turbulent scales. This effect involves many different velocities and temperatures in the same volume of control, violating the assumption of the uniqueness of the particle velocity and temperature distributions. An extension of the statistical approach accounting for non-isothermal conditions may be found in Ref. [START_REF] Masi | The mesoscopic Eulerian approach for evaporating droplets interacting with turbulent flows[END_REF]. Particle velocity and temperature are partitioned in two contributions: i) an Eulerian particle (velocity and temperature) field, referred to as mesoscopic field, which is spatially correlated and shared by all the particles and which accounts for correlations between particles and between particles and fluid; ii) a spatially-uncorrelated particle (velocity and temperature) contribution, referred to as random uncorrelated motion (RUM) contribution, associated with each particle and resulting from the chaotic behaviour of the particles. The RUM contribution is characterized in terms of Eulerian fields of particle velocity and temperature moments; the larger is the particle inertia the more important is the RUM. The proposed ACBMM computes low-order moments of the conditional probability density function (PDF) using a set of governing equations which are derived from the PDF kinetic equation in the general frame of the kinetic theory of dilute gases. It consists of a set of three equations describing the evolution of the mesoscopic particle number density n p , the mesoscopic velocity u p,i and the mesoscopic temperature T p , which needs to be closed by providing constitutive relations as algebraic closures. Neglecting gravity and radiant sources and assuming non-interacting and non-colliding spherical particles in translation, with a same diameter equal or smaller than the Kolmogorov lengthscale, a density much larger than that of the fluid, stationaries inter-phase heat transfers and constant heat capacity, the set of Eulerian equations is [START_REF] Simonin | On the spatial distribution of heavy particle velocities in turbulent flow: From continuous field to particulate chaos[END_REF][START_REF] Masi | The mesoscopic Eulerian approach for evaporating droplets interacting with turbulent flows[END_REF]

∂ n p ∂t + ∂ n p u p,i ∂x i = 0, (1)
∂ n p u p,i ∂t + ∂ n p u p,i u p,j ∂x j = - n p τ p ( u p,i -u f,i ) - ∂ n p δR p,ij ∂x j , (2)
∂ñ p T p ∂t + ∂ñ p ũp,j T p ∂x j = - ñp τ θ T p -T f - ∂ñ p δΘ p,j ∂x j , (3)
where u f,i and T f are the fluid velocity and temperature and

τ p = (1/τ p )|x p (t) = x; H f -1
and τ θ = (1/τ θ )|x p (t) = x; H f -1 are the mesoscopic dynamic and thermal response times obtained by conditional ensemble averaging the particle response times

τ p = 4ρ p d p 3ρ f C D ||u p -u f @p || , τ θ = P r ρ p d 2 p C pp 6Nu µ f C pf ; (4)
u f @p is the fluid velocity at the particle location, P r is the Prandtl number, µ f and ν f are the dynamic and kinematic fluid viscosity respectively, ρ f is the density of the fluid, ρ p and d p are the particle density and the particle diameter, and C pp /C pf is the particle-to-fluid heat capacity ratio. C D and Nu are the corrected drag coefficient [START_REF] Schiller | A drag coefficient correlation[END_REF] and Nusselt number [START_REF] Marshall | Evaporation from drops: II[END_REF] C D = 24 Re p (1 + 0.15Re 0.687 p ), Nu = 2 + 0.552Re 1/2 p P r 1/3 , (

formulated in terms of the particle Reynolds number, Re p = ||u pu f @p ||d p /ν f . The highestorder moment appearing in Eq. ( 2) is the second-order velocity moment referred to as RUM particle kinetic stress tensor (RUM-ST), δR p,ij (x, t) = δu p,i (t)δu p,j (t)|x p (t) = x; H f . The highest-order unclosed moment in Eq. ( 3) is the velocity-temperature correlation δΘ p,i (x, t) = δu p,i (t)δT p (t)|x p (t) = x; H f referred to as RUM particle heat flux (RUM-HF) for the sake of simplicity (although this term has no longer units of a heat flux since the particle mass and heat capacity were removed because constant quantities). The first term on the right hand side (r.h.s.) of Eq. ( 2) accounts for the effects of the drag force and the second term is the transport due to the RUM-ST. The first term on the r.h.s. of Eq. ( 3) represents the interphase heat exchanges and the second contribution is the transport of the mesoscopic temperature due to the RUM-HF. At the second order, the modeling leads to an equation for the RUM-ST [START_REF] Moreau | Modélisation numérique directe et des grandes échelles des écoulements turbulents gaz-particules dans le formalisme eulérien mésoscopique[END_REF] ∂ n

p δR p,ij ∂t + ∂ n p δR p,ij ũp,k ∂x k = -2 n p τ p δR p,ij -n p δR p,jk ∂ u p,i ∂x k -n p δR p,ik ∂ u p,j ∂x k - ∂ ∂x k n p δQ p,ijk , (6)
and an equation for the RUM-HF [START_REF] Masi | Theoretical and numerical study of the modeling of unsteady non-isothermal particle-laden turbulent flows by an Eulerian-Eulerian approach[END_REF]

∂ ∂t ñp δΘ p,i + ∂ ∂x j ñp ũp,j δΘ p,i = -ñp 1 τ p + 1 τ θ δΘ p,i -ñp δΘ p,j ∂ ũp,i ∂x j -ñp δR p,ij ∂ T p ∂x j - ∂ ∂x j ñp δ∆ p,ij . (7)
The highest-order moments appearing in Eqs. ( 6)-( 7) are third-order RUM correlations defined by the conditional ensemble average as δQ p,ijk (x, t) = δu p,i (t)δu p,j (t)δu p,k (t)|x p (t) = x; H f and δ∆ p,ij (x, t) = δu p,i (t)δu p,j (t)δT p (t)|x p (t) = x; H f . The first term on the r.h.s. of Eq. ( 6) represents the dissipation of the RUM-ST while second and third terms are productions by mesoscopic velocity gradients. The first term on the r.h.s. of Eq. ( 7) accounts for the dissipation of the RUM-HF due to heat and momentum transfers and second and third terms are productions by both mesoscopic velocity and temperature gradients. In the framework of the ACBMM approach proposed in this study, Eqs. ( 6) and ( 7) are not resolved but rather used for the development of the second-order moment closures. Preliminary two-dimensional Eulerian-Eulerian DNS using modeled RUM-ST in conjunction with resolved RUM-HF may be found in Ref. [START_REF] Dombard | Modeling heat transfer in dilute two-phase flows using the mesoscopic Eulerian formalism[END_REF]. With the purpose of deriving algebraic closures, the equation of the evolution of the RUM kinetic energy (δθ p (x, t) = 0.5 δu p,i (t)δu p,i (t)|x p (t) = x; H f ) and the RUM temperature variance (δθ θ (x, t) = 0.5 δT p (t)δT p (t)|x p (t) = x; H f ) are also written [START_REF] Simonin | On the spatial distribution of heavy particle velocities in turbulent flow: From continuous field to particulate chaos[END_REF][START_REF] Masi | Theoretical and numerical study of the modeling of unsteady non-isothermal particle-laden turbulent flows by an Eulerian-Eulerian approach[END_REF]:

∂ n p δθ p ∂t + ∂ n p δθ p ũp,m ∂x m = -n p δR p,nm ∂ u p,n ∂x m -2 n p τp δθ p - 1 2 ∂ ∂x m n p δQ p,nnm , (8)
∂ n p δθ θ ∂t + ∂ n p δθ θ ũp,m ∂x m = -n p δΘ p,m ∂ T p ∂x m -2 n p τθ δθ θ - 1 2 ∂ ∂x m n p δΩ p,m . (9)
The last term in Eq. ( 9) is defined as δΩ p,m (x, t) = δu p,m (t)δT p (t)δT p (t)|x p (t) = x; H f . The ACBMM bases its efficiency and accuracy on the algebraic closures provided for modeling the RUM-ST and the RUM-HF; in this study, this concern is addressed.

Modeling The RUM Particle Kinetic Stress Tensor (RUM-ST)

The RUM stress tensor is composed of spherical and deviatoric contributions

δR p,ij = 1 3 δR p,kk δ ij + δR * p,ij = 2 3 δθ p δ ij + δR * p,ij . (10)
The spherical part may be obtained by means of an additional transport equation [START_REF] Simonin | On the spatial distribution of heavy particle velocities in turbulent flow: From continuous field to particulate chaos[END_REF][START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF] while the deviatoric contribution needs to an algebraic closure. In this work the deviatoric RUM-ST is modeled by using a polynomial representation for tensor functions [START_REF] Masi | Theoretical and numerical study of the modeling of unsteady non-isothermal particle-laden turbulent flows by an Eulerian-Eulerian approach[END_REF]. The polynomial representation is used in the framework of an assumption of equilibrium of the local anisotropy. Indeed, one of the most important finding of this study was to observe that the RUM-ST is a self-similar tensor which means that its temporal and spatial evolutions are related to that of its trace, involving equilibrium of anisotropy. Such an equilibrium leads to an algebraic stress model (ASM) which is an implicit and nonlinear system of equations referred to as 2ΦASM (by analogy with turbulence models, see, e.g. [START_REF] So | An explicit algebraic Reynolds stress and heat flux model for incompressible turbulence: Part II Buoyant flow[END_REF] and references cited in). It is obtained defining a local RUM anisotropy tensor as

b * p,ij = δR p,ij 2δθ p - 1 3 δ ij . (11)
Then the equilibrium involves

D Dt b * p,ij = 0 with D Dt = ∂ ∂t + u p,j ∂ ∂x j , (12)
which using the definition [START_REF] Druzhinin | Direct numerical simulations bubble-laden turbulent flows using the two-fluid formulation[END_REF] gives the relation

D Dt δR p,ij = δR p,ij δθ p D Dt δθ p . (13)
Injecting Eqs. ( 6) and (8) into Eq. ( 13), and assuming equality between left hand side (l.h.s.) and r.h.s. third-order correlations, the equation takes the form

δR p,ij - δR p,nm 2δθ p ∂ u p,n ∂x m = - 1 2 δR p,jk ∂ u p,i ∂x k - 1 2 δR p,ik ∂ u p,j ∂x k . ( 14
)
Defining the mesoscopic particle rate-of-strain and vorticity tensors as, respectively,

S p,ij = 1 2 ∂ u p,i ∂x j + ∂ u p,j ∂x i , Ω p,ij = 1 2 ∂ u p,i ∂x j - ∂ u p,j ∂x i , (15)
and using the decomposition of S p,ij in deviatoric and spherical contributions S p,ij = S * p,ij + 1 3 S p,kk δ ij and the RUM anisotropy definition, Eq. ( 14) may be rearranged as

b * p,ij -2b * p,nm S * p,nm = - 2 3 S * p,ij -b * p,ik S * p,kj + S * p,ik b * p,kj - 2 3 b * p,nm S * p,nm δ ij + b * p,ik Ω p,kj -Ω p,ik b * p,kj , (16)
which represents the implicit nonlinear 2ΦASM model proposed for closing the RUM-ST. Explicit solutions are then provided using techniques well known in turbulence [START_REF] Gatski | On explicit algebraic stress models for complex turbulent flows[END_REF][START_REF] Girimaji | Fully explicit and self-consistent algebraic Reynolds stress model[END_REF][START_REF] Wallin | An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows[END_REF]. Hereinafter, bold notation denotes three-dimensional second-rank tensors, curly brackets {.} represent the tensor trace and the asterisk means traceless tensor. The matrix multiplication is then defined in a matrix notation as C = AB = A ik B kj = C ij and B 2 = BB. According to Pope [START_REF]POPE A more general effective-viscosity hypothesis[END_REF], the anisotropy tensor may be expressed using a polynomial representation for tensor functions b * = ς G ς T (ς) which represents the linear combination of a set of non-dimensional independent, symmetric and deviatoric second-order tensors T (ς) , using scalar coefficients G ς which are functions of the invariants of the dimensionless S + and Ω + . Using the Cayley-Hamilton theorem, Pope [START_REF]POPE A more general effective-viscosity hypothesis[END_REF] showed that a set of ten (ς = 10) tensors T (ς) is needed to form an integrity basis [START_REF] Spencer | Theory of Invariants[END_REF] in order to express every symmetric deviatoric three-dimensional secondorder tensor formed by S + and Ω + ; the problem then reduces to model the ten coefficients G ς (see, e.g., Refs. [START_REF] Gatski | On explicit algebraic stress models for complex turbulent flows[END_REF][START_REF] Wallin | An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows[END_REF]). However, the integrity basis and coefficients associated with may be reduced if some approximations are introduced. Using a two-dimensional flow approximation (three tensor basis) Girimaji [START_REF] Girimaji | Fully explicit and self-consistent algebraic Reynolds stress model[END_REF] developed a self-consistent solution technique in order to model the Reynolds stresses in turbulent flows. In this work, we will use the same technique applied to our 2ΦASM model. The algebraic, implicit and nonlinear system for the dispersed phase is rewritten using the Girimaji's notation as

b * L 0 1 -L 1 1 {b * S + } = L 2 S + + L 3 b * S + + S + b * - 2 3 {b * S + }I -L 4 b * Ω + -Ω + b * ,( 17
)
where

L 0 1 = 0, L 1 1 = 2, L 2 = -2 3 , L 3 = -1 and L 4 = -1, S +
and Ω + are dimensionless tensors by the quantity II 1/2 S = {S * 2 } 1/2 and I is the identity matrix. According to Girimaji [START_REF] Girimaji | Fully explicit and self-consistent algebraic Reynolds stress model[END_REF], the general representation of the anisotropy tensor under the two-dimensional approximation is

b * = 3 ς=1 G ς T (ς) = G 1 S + + G 2 S + Ω + -Ω + S + + G 3 S +2 - 1 3 {S +2 }I , (18)
where the three coefficients are functions of the two invariants η 1 = {S +2 } and η 2 = {Ω +2 }. Using Eq. ( 18) and the two-dimensional hypothesis, the contracted product b * p,ij S + p,ij may be written as

b * p,ij S + p,ij = G 1 S + p,ij S + p,ij = G 1 η 1 . (19)
Inserting equations ( 18) and (19) into the system (17) leads to

[G 1 S + + G 2 (S + Ω + -Ω + S + ) + G 3 (S +2 - 1 3 {S +2 }I)](L 0 1 -η 1 L 1 1 G 1 ) = (20)
(L 2 + 1 3 η 1 L 3 G 3 -2η 2 L 4 G 2 )S + + 2L 3 G 1 (S +2 - 1 3 {S +2 }I) -L 4 G 1 (S + Ω + -Ω + S + ).
Then, comparison between homogeneous l.h.s. and r.h.s. terms in Eq. ( 20) leads to explicit solutions for the unknown coefficients. In our particular case in which L 0 1 = 0, the coefficients G 2 and G 3 are directly obtained as follows

G 2 = L 4 η 1 L 1 1 , G 3 = - 2L 3 η 1 L 1 1 , (21)
and G 1 reduces to the solution of a pure quadratic

G 2 1 = - 1 η 1 L 1 1 [L 2 + 1 3 η 1 L 3 G 3 -2η 2 L 4 G 2 ]. (22)
Eq. ( 22) with the "L" coefficients leads to G 1 = ± √ 2η 1 + 2η 2 /2η 1 , which admits real solutions only for η 1 + η 2 ≥ 0. Unfortunately, as η 1 and η 2 have positive and negative sign respectively, real solutions are not ensured. In order to use this model, local negative values of the discriminant are set to zero. The legitimacy of such an approximation was investigated in Ref. [START_REF] Masi | Theoretical and numerical study of the modeling of unsteady non-isothermal particle-laden turbulent flows by an Eulerian-Eulerian approach[END_REF]; it was observed that negative values correspond to small magnitudes of both the particle rateof-strain and vorticity tensors justifying the approximation G 1 = 0. Concerning the sign of the coefficient G 1 , according to Eq. ( 19) and the definition of the normalized production of the RUM particle kinetic energy by shear (-{b * S * }, from Eq. ( 8)), in the simplest case it should be taken as negative. Negative sign for G 1 is obtained in turbulence, in the domain of applicability of an equilibrium assumption of anisotropy [START_REF] Girimaji | Fully explicit and self-consistent algebraic Reynolds stress model[END_REF]. However, for the dispersed phase interacting with turbulent flows, it is usual to have reverse energy exchanges from the RUM to the mesoscopic contribution [START_REF] Moreau | Development of gas-particle Euler-Euler LES approach: A priori analysis of particle sub-grid models in homogeneous isotropic turbulence[END_REF][START_REF] Masi | The mesoscopic Eulerian approach for evaporating droplets interacting with turbulent flows[END_REF] corresponding to a reverse sign of the first-order coefficient G 1 . In the special case in which the deviatoric RUM-ST and the deviatoric particle rate-of-strain are axisymmetric tensors (as reported in Refs. [START_REF] Masi | Theoretical and numerical study of the modeling of unsteady non-isothermal particle-laden turbulent flows by an Eulerian-Eulerian approach[END_REF][START_REF] Masi | Modeling of the randomuncorrelated velocity stress tensor for unsteady particle Eulerian simulationin turbulent flows[END_REF]), assuming that the principal direction of the RUM-ST is aligned with the most extensive or compressive direction of the particle rateof-strain, then the sign of the contracted product b * p,ij S + p,ij , which is the same of G 1 , is given by the sign of the third invariant of the particle rate-of-strain tensor defined as η 3 = {S +3 }. In this case the solution becomes

G 1 = sign(η 3 ) √ 2η 1 + 2η 2 /2η 1 . Injecting the coefficients G 1 ,
G 2 , G 3 into Eq. ( 18) leads to an explicit solution for the RUM anisotropy tensor and thus for the RUM-ST. Hereinafter, this model will be referred to as 2ΦEASM. Results from preliminary isothermal Eulerian-Eulerian DNS using 2ΦEASM may be found in Ref. [START_REF] Sierra S Ánchez | Modeling the dispersion and evaporation of sprays in aeronautical combustion chambers[END_REF].

Modeling The RUM Particle Heat Flux (RUM-HF)

In a similar way, the RUM-HF is modeled assuming similarity between the evolution of the RUM heat flux and that of the square root correlation between the RUM kinetic energy (δθ p ) and the RUM temperature variance (δθ θ ), which involves the equilibrium of the normalized heat flux. This assumption may be written as which developed leads to

D Dt δΘ p,i δθ p δθ θ = 0 (23)
1 δθ p δθ θ D Dt δΘ p,i = 1 2 δΘ p,i δθ p δθ θ 1 δθ p D Dt δθ p + 1 2 δΘ p,i δθ p δθ θ 1 δθ θ D Dt δθ θ . (24)
Then injecting Eqs. ( 7), ( 8) and ( 9) into Eq. ( 24) and assuming equality between l.h.s. and r.h.s. diffusion terms, the equation takes the form

δΘ p,i δθ p δθ θ - 1 2 δΘ p,l δθ θ ∂ T p ∂x l - 1 2 δR p,kl δθ p ∂ u p,k ∂x l = - δR p,il δθ p δθ θ ∂ T p ∂x l - δΘ p,l δθ p δθ θ ∂ u p,i ∂x l , (25)
which may be rearranged as follows

δH p,i - 1 2 δH p,l Kp,l -b * p,kl S * p,kl = -2 b * p,il + 1 3 δ il Kp,l -(S * p,il + Ω p,il )δH p,l , (26)
by defining the quantities

δH p,i = δΘ p,i δθ p δθ θ , Kp,i = δθ p δθ θ ∂ T p ∂x i . (27)
A fully dimensionless system may be obtained using non-dimensionalized quantities Kp,i , S * p,ij

and Ω p,ij by the same inverse timescale, then referred to as K + p,i , S + p,ij and Ω + p,ij . Eq. ( 26) is an implicit and nonlinear particle algebraic heat flux model (AHFM), referred to as 2ΦAHFM (by analogy with turbulence models). Explicit self-consistent solutions, 2ΦEAHFM, are then provided using a technique suggested for turbulent heat flux by Wikström et al. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF] and by introducing a regularization procedure for ensuring non-singular solutions during the process of matrix inversion. Eq. ( 26) was indeed rearranged in a similar way that the analogous equation of Wikström et al. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF] in order to straightforwardly apply their technique to our 2ΦAHFM. Using the same notation as in Ref. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF], we define the scalar quantity into the parentheses on the l.h.s. of Eq. [START_REF] So | An explicit algebraic Reynolds stress and heat flux model for incompressible turbulence: Part II Buoyant flow[END_REF] as

N θ = - 1 2 δH p,l K+ p,l -b * p,kl S + p,kl . (28)
The system (26) may then be solved by inverting the matrix

A p,ij = N θ δ ij + S + p,ij + Ω + p,ij , (29)
obtaining the explicit solutions

δH p,i = -2A -1 p,ij b * p,jl + 1 3 δ jl K+ p,l . (30)
According to the Cayley-Hamilton theorem, an analytic expression for A -1 may be written (see, e.g., Ref. [START_REF] Pope | Turbulent Flows[END_REF] ) and the operation of matrix inversion is ensured provided that A admits real non-zero eigenvalues so to avoid the occurrence of singularities when its determinant (denominator of the analytical solution) vanishes. The technique of Wikström et al. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF] involves injecting the relation of the flux (Eq. 30 in our case) into the definition of N θ (Eq. 28 in our case), replacing A -1 by its analytical expression. For fully three-dimensional flows, this leads to a fourth-order polynomial in N θ difficult to handle. An alternative is to resolve a thirdorder polynomial corresponding to a two-dimensional flow approximation and then using the solution for three-dimensional flows as well [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF]. Assuming two-dimensional flows (but using three-dimensional tensors) the analytic expression for the inverse of the matrix is

A -1 = N θ I -(S + + Ω + ) N 2 θ -1 2 (η 1 + η 2 ) . (31)
Injecting Eq. ( 31) into Eq. ( 30) and the latter into the definition (28) the following third-order polynomial is obtained

2N 3 θ + 2{b * S + }N 2 θ -(Q 1 + R 1 )N θ -{b * S + }Q 1 + R 2 = 0, (32)
where Q 1 = η 1 + η 2 and the terms R 1 and R 2 are, respectively,

R 1 = 2 b * p,il + 1 3 δ il K+ p,l K+ p,i R 2 = 2 S + p,ij + Ω + p,ij b * p,jl + 1 3 δ jl K+ p,l K+ p,i . (33)
A unique real solution is then ensured using transcendental functions as in Ref. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF]. At this point, we need to address the concern of stable solutions by ensuring the determinant of the matrix does not vanish and singularities do not occur. From the a priori analysis, we observed that the determinant is mainly positive and that negative values correspond to small magnitudes of it. In order to avoid singular solutions the determinant should not change its sign and the following condition should be satisfied

det(A) = N 2 θ - 1 2 (η 1 + η 2 ) > 0 (34)
while, in general, it is not. Solutions proposed in Ref. [START_REF] Öm | Derivation and investigation of a new explicit algebraic model for the passive scalar flux[END_REF] for turbulent flows do not apply to the dispersed phase and a regularization procedure is therefore necessary. The reciprocal of the determinant is re-written as follows where η = √ η 1 , ξ = √ -η 2 and x 2 = η 2 /(2N 2 θ ). Then, using the first order MacLaurin series of the function 1/(1 + x 2 ) we write x 2 ≈ 1 -1/(1 + x 2 ) which inserted into Eq. ( 35) leads to

det(A) -1 = 1 N 2 θ -1 2 (η 1 + η 2 ) = 2/(2N 2 θ ) 1 -η 2 /(2N 2 θ ) + ξ 2 /(2N 2 θ ) = 2/(2N 2 θ ) 1 + ξ 2 /(2N 2 θ ) -x 2 (35)
det(A) -1 = 2(2N 2 θ + η 2 ) (2N 2 θ ) 2 + ξ 2 (2N 2 θ + η 2 ) . ( 36
)
The stability is thus ensured as the denominator of Eq. ( 36) never vanishes for N θ and ξ 2 or η 2 = 0. The same regularization procedure was used by Gatski and Speziale [START_REF] Gatski | On explicit algebraic stress models for complex turbulent flows[END_REF] for modeling the Reynolds stress tensor in the framework of the EASM approach. Hereinafter, explicit regularized solutions for the RUM-HF will be referred to as 2ΦEAHFM.

A Priori Analysis

An evaluation of the ACBMM by assessing the proposed algebraic closures is given by an a priori analysis using particle Eulerian fields which are extracted from a Eulerian-Lagrangian DNS of a temporal particle-laden non-isothermal turbulent planar jet. The numerical simulation corresponds to the dispersion of a cold particle-laden jet into a hot homogeneous isotropic decaying turbulence [START_REF] Simonin | Numerical study of the inter-phase heat transfers in a nonisothermal particle laden slab flow[END_REF]. The simulation domain is a cube of length size L box = 2πL ref and mesh 128 3 cells. The initial gas velocity has hyperbolic-tangent profile supplemented with statistically homogeneous and isotropic velocity fluctuations [START_REF] Vermorel | Numerical study and modelling of turbulence modulation in a particle laden slab flow[END_REF]. The initial gas temperature has hyperbolic-tangent mean profile and zero fluctuations. As we assume very dilute regime, collisions and turbulence modulation are not accounted for. Within the slab, of width d = 0.25L box , 13 millions of solid particles are randomly embedded at the same mean velocity and temperature as the carrier flow and zero fluctuations. The simulation corresponds to a Stokes number St ∼ 1 computed over a characteristic timescale of the turbulence seen by the particles representing an estimate of the integral Lagrangian timescale T L in mean-sheared flows [START_REF] Deutsch | Large eddy simulation applied to the motion of particles in stationary homogeneous fluid turbulence[END_REF]. The Reynolds number based on the energetic lengthscale is Re Le ∼ 73. The initial set of dimensionless parameters is: kinematics viscosity detailed in Ref. [START_REF] Kaufmann | Comparison between Lagrangian and mesoscopic Eulerian modelling approaches for inertial particles suspended in decaying isotropic turbulence[END_REF]. The evaluation of the accuracy of the models is performed at tensor/vector level by assessed each component of the RUM-ST and the RUM-HF against the actual quantity provided by the particle Eulerian database. Moreover, an assessment at scalar level is also given. Figure 1 shows the mean profiles of the RUM-ST components as predicted by 2ΦEASM compared with the actual stresses. Modeled, by 2ΦEAHFM, and actual mean profiles of the RUM-HF components are depicted in Figure 2. For both the second-order moments, results at tensor/vector level are very satisfactory. At scalar level, the scalar quantities investigated are the productions of the RUM kinetic energy and the RUM temperature variance (first terms on the r.h.s. of Eqs. ( 8) and (9), respectively). Figure 3 shows the mean profiles of the two productions. In Figure 4, their PDFs evaluated at the periphery of the jet are depicted. They are computed by multiplying the local value of each modeled production by the actual-to-modeled magnitude ratio; in this way, only the shape of the PDFs is investigated [START_REF] Moreau | Modélisation numérique directe et des grandes échelles des écoulements turbulents gaz-particules dans le formalisme eulérien mésoscopique[END_REF]. Figure 3 shows that the two productions are successfully predicted, even if the RUM kinetic energy production is somewhat overestimated. Figure 4 shows that both the models are able to match the shape of the actual PDF and to reproduce reverse energy exchanges. Globally, results show that models compare very well with the Eulerian-Lagrangian DNS and properly reproduce all crucial trends of the computation. The proposed algebraic closures enable the ACBMM to correctly predict unsteady non-isothermal particle-laden turbulent flows, in dilute regime, for large Stokes number as well.

ν f = 1.82e-4u ref L ref , turbulent kinetic energy q 2 f = 3.37e-4u 2

Figure 1 :

 1 Figure 1: Profiles of actual (solid line) and modeled (line with symbols) mean deviatoric RUM stresses, at the time t = 6.2T f @p .

Figure 2 :

 2 Figure 2: Profiles of actual (solid line) and modeled (line with symbols) mean streamwise and normalwise RUM-HF components, at the time t = 6.2T f @p .

2 - 2 -Figure 3 :

 223 Figure 3: Mean profiles of actual (solid line) and modeled (line with symbols) RUM temperature variance production (left) and RUM kinetic energy production (right), at the time t = 6.2T f @p .

Figure 4 :

 4 Figure4: PDF of actual (solid line) and modeled (symbols) productions of RUM temperature variance (left) and RUM kinetic energy (right), at the periphery of the jet, at t = 6.2T f @p .

Acknowledgements

This work received funding from the European Community through the TIMECOP-AE project (Project # AST5-CT-2006-030828). It reflects only the author's views and the Community is not liable for any use that may be made of the information contained therein. Numerical simulations were performed by the IBM Power6 machine; support from Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) is gratefully acknowledged.


M.-C Gauffre

H Neau

O Simonin

R Ansart

N Meyers

S Petitot

Numerical simulation of a 3D unsteady two-phase flow in the filling cavity in oxygen of a cryogenic rocket-engine

The feeding of the LOX dome of a cryogenic rocket-engine is a decisive stage of the transient engine ignition. However flight conditions are difficult to reproduce by experimental ground tests. The work reported here is part of an ongoing research effort to develop a robust method for prediction and understanding the LOX dome feeding. In the framework of this project, experiments with substition fluids (air and water) are conducted, without mass and energy transfer. This work presented here intends to reproduce these experiments through incompressible two-phase flow CFD simulations, in an industrial geometry equivalent to the experimental mock-up, made up of a feeding piper, a dome and 122 injectors. More precisely, the aim is to compare the numerical results obtained with NEPTUNE CFD code with the experimental results, through the dome pressure and the mass flow rate of water at the outlet. An important work was made to obtain the same inlet conditions in NEPTUNE CFD code as the experimenters, in order to compare the numerical results with the experimental results for the best. The influence of the interfacial momentum transfer modeling and turbulence modeling are also studied here. The turbulence modeling plays no macroscopic or local role on the mass flow rate of water, on the mass of water in dome and on the dome pressure. The drag model has a major impact on our results as well globally as locally, unlike the turbulence modeling. The Simmer-like model is prefered in comparison to the Large Interface called LIM, because it is in better agreement with experimental data. Moreover, it has to be highlighted that the Simmer-like model is very sensitive to its parameter d, the inclusion diameter.

I. Introduction

The feeding of the LOX dome of a cryogenic rocket-engine is a decisive stage. Indeed, the flow in the cavity governs the way the combustion chamber is supplied from the opening to the closure of the main valve. In particular, in the first moments of the filling of the cavity, the distribution of the flow in the injection plate region controls the ignition phase of the engine. Thus, it is necessary to know the transient flow within the injection cavity in order to provide the correct boundary conditions to simulate correctly the ignition of a cryogenic engine. As soon as the main valve is opened, the liquid oxygen enters the cavity, already swept by a helium venting at ambient temperature. In contact with the cavity walls and the helium, the liquid oxygen vaporizes strongly at first and then in a more moderate way afterwards. Therefore, the topology of the flow within the cavity evolves significantly in the first instants of its filling. Moreover, in the case of an upper-stage rocket-engine, the filling is affected by both the flashing phenomenon and the microgravity which can modify the heat transfers. Flight conditions are difficult to reproduce by experimental ground tests: CNES (Centre National d' Études Spatiales) and SAFRAN Snecma set up a research program based on both experimental and numerical studies.

It involves comparing numerical results obtained with CFD codes with the experimental work performed at LEGI (Geophysical and Industrial Flows Laboratory at Grenoble, France). Several CFD codes were foreseen: on the one hand, industrial codes like Fluent, Star-CD or CFX and on the other hand, R&D codes like NEPTUNE CFD, CEDRE, AVBP or LEONARD. Finally two R&D codes were retained: LEONARD 1,2 code developed at Polytech Marseille (France) and NEPTUNE CFD code developed at IMFT (Fluid Mechanics Institute at Toulouse, France). Previous works [START_REF] Ansart | Dossier d'analyse physique du code NEPTUNE CFD : Remplissage de la cavit d'injection[END_REF] showed good potential of NEPTUNE CFD code for the simulation of the feeding of the LOX dome. In this paper the numerical results obtained with NEP-TUNE CFD code are presented. The framework is the following: substition fluids (air and water) are used and the cases studied are adiabatic, without mass and energy transfer. The aim of this research work is to simulate the filling of the cavity.

II. Experimental setup

An experimental program has been set up at LEGI in order to study the flow in the cavity, without flashing phenomenon. The matter at issue is both to understand and to identify the phenomena at stake in the transient feeding of the cavity and to build a database, so as to validate the developments in NEPTUNE CFD models. Thus, a series of test campaigns was carried out on an experimental test bed set up at LEGI.

A. Description of the experimental setup

The LOX cavity of the rocket-engine represented in a simplified manner at LEGI by the mock-up with a toric volume which keeps the volume of the real cavity can be seen in figure 1. In addition, this experimental model is made up of a cavity (also called dome), a feeding pipe upstream from the cavity, 122 injectors downstrean from the cavity and an igniter pipe in its centre.

B. Experimental conditions

The experiment consists of a representative model of the injection cavity filled with substitution fluids respecting Weber number defined by W e = , where ρ 1 is the density of water, V max is the supposed maximum inlet velocity of water when the valve is completely open, D is the diameter of the valve and σ is the surface tension. Water is the filling fluid, whereas air is the fluid initially in the cavity and also used to sweep this former. Moreover, these experiments were carried out without heat and mass transfer: dome wall are not heated and the two fluids are at the same temperature. Water flow at feeding pipe inlet level is controled by a spherical bushel valve, which ideally opens between 0°and 90°linearly during the 100 ms of the opening stage. Then during the plateau stage of 500 ms, the valve remains open at its maximum angle 90°, before closing in 500 ms until the complete closing. It has to be mentioned that during the motion of the bushel valve, the liquid entry only begins after a valve rotation of 19.6°from the closed position of the valve: this is the real beginning of the valve.

C. Measurements carried out at LEGI

During the experiments at LEGI, 4 several measurements were carried out with a complete opening of the bushel valve. Pressure was measured in the cavity. Moreover, imaging techniques in white light were used to visualize the flow in the transparent, PMMA-made cavity. Imaging techniques by laser sheet also permitted to visualize the injectors outlet: the laser sheets detect if water comes out of the injector outlet and therefore it is possible to determine if the injectors are "wet" or not. It is to be noted that this process can not detect more than 70 "wet" injectors, for the moment.

III. Description of CFD model

The three dimensional numerical simulations were performed with NEPTUNE CFD v1.08.

A. Solver and models

Overall presentation of NEPTUNE CFD

NEPTUNE CFD is developed in the framework of the NEPTUNE project, financially supported by CEA (Commissariat à l' Énergie Atomique), EDF ( Électricité de France), IRSN (Institut de Radioprotection et de Sûreté Nucléaire), and AREVA-NP. This code is a Finite Volume Eulerian multi-phase solver parallelized 5 designed for nuclear engineering.

The behavior of a fluid continuum made of several physical phases or components can be modeled using the general Eulerian multi-field balance equations [START_REF] Drew | Mathematical modeling of two-phase flow[END_REF][START_REF] Simonin | Modélisation numérique des écoulements turbulents diphasiques à inclusions dispersée, École de printemps de Mécanique des Fluides Numérique[END_REF][START_REF] Simonin | Eulerian prediction of a turbulent bubbly flow in a sudden pipe expansion[END_REF] . Considering the adiabatic case without mass and energy transfer, the two-fluid balance equations (mass conservation and momentum conservation) are obtained from the fundamental conservation laws of physics. In our simulations, k = 1 is for liquid phase and k = 2 is for gas phase.

In NEPTUNE CFD code [START_REF] Méchitoua | Numerical simulation of multiphase flow with an elliptic[END_REF]13 , the conservation laws are written in a classical differential form that is valid at any time and location within the continuum, except across the interfaces between two physical phases. At the interfaces, jump conditions are derived from the continuous equations.

Balance equations

Mass balance

The multi-fluid mass balance equation for the field k writes:

∂ ∂t (α k ρ k ) + ∂ ∂x i (α k ρ k U k,i ) = 0, (1)
where α k is the volume fraction of phase k, ρ k is the density of phase k and U k,i is the i th component of the mean velocity of phase k. Total volume conservation leads to k α k = 1.

Momentum balance equation

The multi-fluid momentum balance equation for the field k is written in its semi-conservative form (all the contributions are conservative, except the pressure gradient one) as:

∂ ∂t (α k ρ k U k,i ) + ∂ ∂x j (α k ρ k U k,i U k,j ) = ∂ ∂x j (α k τ k,ij + Σ Re k,ij ) -α k ∂P ∂x i + α k ρ k g i + p =k I (p→k) + α k S k , (2)
where:

• P is the mean pressure, µ k is the dynamic viscosity and g i is the acceleration due to gravity.

• τ k,ij = µ k ∂Ui ∂xj + ∂Uj ∂xi -2 3 div(U )δ ij is the viscous stress tensor. • Σ Re k,ij = -α k ρ k U ′ k,i U ′ k,j k
is the turbulent stress tensor.

• S k is the i th component of the external source term of head losses.

• I (p→k) represents the i th component of the average interfacial momentum transfer rate from phase p to phase k, that accounts for the drag force. It verifies I (p→k) + I (k→p) = 0.

Closure laws in NEPTUNE CFD

After having averaged the instantaneous local balance equations, the two-phase flow system has four main unknowns in our case: the volume fraction α k of phase k, the mean pressure P and the mean velocity U k of phase k. Two terms resulting from the transition to the time-average of the instantaneous local balance equations have to be modeled in NEPTUNE CFD in order to resolve the balance equations: the turbulent stress tensor Σ Re k,ij and the interfacial momentum transfer rate I (p→k) .

Turbulence modeling In our case, the k-ε turbulence model is activated for both phases. The Reynolds stresses tensor is closed for each phase, using a Boussinesq-like hypothesis:

ρ k U ′′ k,i U ′′ k,j k = -µ t k ∂U k,i ∂x j + ∂U k,j ∂x i + 2 3 δ ij ρ k q 2 k + µ t k ∂U k,m ∂x m (3)
where µ t k is the turbulent viscosity and

q 2 k = 1 2 U ′′ k,i U ′′ k,i k
is the turbulent kinetic energy of phase k.

As the flow is complicated, this kind of model should be sufficient, in particular for the stratified flow.

Interfacial momentum transfer and drag closure law In our case, only the drag force is activated for the modeling of the average interfacial momentum transfer rate I (p→k) .

The standard choice provided by NEPTUNE CFD is the Simmer-like model. [START_REF] Méchitoua | Assessment of NEPTUNE CFD code for some free Surface Flows interesting Fluvial Hydraulic[END_REF] The Simmer-like law is fitted to have a physical behaviour in the limits: it considers either dispersed gas bubbles in a continuous liquid flow or dispersed liquid droplets in a continuous gas flow with regard to the volume fraction. Thus, the Simmer-like law corresponds to the liquid droplets drag law in a continuous gas flow when the volume fraction is lower than 0.3 and to the gas bubbles drag law in a continuous liquid flow when the volume fraction of water is greater than 0.7. For intermediate volume fractions, the drag law is calculated by cubic interpolation between these two limits. The drag force F by mass unit is given by the following equations:

F = F l = - 3 4
ρ g ρ l C Dl d l |U rel |U rel if α l ≤ 0.3, (4)
F = F g = - 3 4 ρ l ρ g C Dg d g |U rel |U rel if α l ≥ 0.7, (5)
F = f (F l (α l = 0.3), F g (α l = 0.7)) if 0.3 ≤ α l ≤ 0.7, (6)
where α l is the volume fraction of liquid, ρ l and ρ g are respectively the density of liquid and gas, d g and d l represent respectively the characteristic diameter of droplets and bubbles, U rel is the relative velocity between the two phases, f defines the cubic interpolation, and C D is the drag coefficient. This coefficient is defined by:

C D = 24ν c |U rel |d p 1 + 0.15 |U rel |d p ν c 0.687 , (7)
where the index p corresponds to the dispersed phase and the index c to the continuous phase, and where ν is the kinematic viscosity. Nonetheless, this model is not adapted to pure stratified flows, because the Simmer-like law may lead to a significant over-evaluation or under-estimation of the friction between the gas and the liquid. Indeed, the drag model based on the friction between a dispersed and a continuous phase is sensitive to the diameter value. So an alternative to model the friction is the Large Interface Model or LIM 7 implemented in NEPTUNE CFD code and validated by EDF R&D from experimental measurements. [START_REF] Fabre | Stratified flows, part i: local structure[END_REF] This method allows to locate the free surface and takes into account momentum and turbulence exchanges between phases. First the free surface is located and built from the local gradient of volume fraction, thanks to the refined-gradient method. Then in this three-cell layer built around the interface, liquid and gas characteristic tangential velocities are determined, in order to calculate interfacial velocity and momentum exchange, considering wall-function on both sides.

B. Numerical setup

For the simulations, incompressible two-phase flow (with water and air) are considered for each phase, without mass and energy transfer. It should be remembered that NEPTUNE CFD code solves the mass balance and the momentum balance equations for each phase. The k-ε turbulence model is activated for both phases. The dome is initially filled with air. Water is injected at pipe inlet level.

Computational domain

The simulations are conducted on an unstructured three-dimensional mesh depicted in figure 2.

This mesh is the result of the merging of non coincident meshes, composed of 1,124,000 cells (mainly hexahedra). Because of the complexity of the geometry, the use of an unstructured mesh is needed. The geometry is equivalent to the experimental mock-up. This industrial geometry is composed of a feeding pipe, a cavity with 122 injectors and an igniter pipe in its centre.

Initial conditions

The mock-up is initially filled with air (in phase k = 2). The density of air is ρ 2 = 1.2 kg/m 3 and the dynamic viscosity of water is µ 2 = 1.85.10 -5 Pa.s.

Boundary conditions

Inlet boundary conditions

Water is injected at pipe inlet, with a density ρ 1 = 1, 000 kg/m 3 and a dynamic viscosity µ 1 = 10 -3 Pa.s. As stated before in this paper, a Weber number is imposed at the inlet:

W e = ρ1V M ax D σ
, where V max is the supposed maximum inlet velocity of water when the valve is completely open, D is the diameter of the valve and σ is the surface tension. If the k-ε turbulence model is activated, this Weber number imposes the turbulent kinetic energy q 2 1 = 3 2 (V M ax I) 2 and the turbulent dissipation

ε 1 = C µ (q 2 1 ) 3/2 l -1
, where I is the turbulent rate imposed at inlet, C µ = 0.09 and l = 0.03D H with D H the hydraulic diameter. This injection takes into account the complete opening between 0°and 90°of the bushel valve to be in the same conditions as the experiments at LEGI. It implies that the evolution of the direction of the velocity-vectors, of the inlet opening surface area and of the mass flow rate have to be considered.

The inlet experimental data are the inlet mass flow rate of water and the opening angle θ of the bushel valve. Thus, there is still to determine the evolution of the opening surface area of this valve. As a first approximation, the 3D geometric problem of the valve opening is turned into a 2D problem in the same plane as the inlet mesh of the feeding pipe (figure 3).

The intersection of a fixed disc and a mobile disc in two different planes becomes the intersection of a fixed disc and a mobile ellipse in the same plan. This ellipse is the projection of the mobile disc onto the plane of the fixed disc and evolves during the time (in particular the position of its centre 0 ′ and the dimension of its half-small axis called r). After some calculations, the area S Exp of the intersection surface is:

S Exp = R 2 cos -1 d R + r - R 2 d R + r sin cos -1 d R + r + 1 2 cos -1 d R + r (r 2 + R 2 ) - R 2 -r 2 4 sin 2cos -1 d R + r - rdR R + r sin cos -1 d R + r , (8)
where:

• θ is the opening angle of the bushel valve,

• D is the diameter of the valve, • R = D 2 is the radius of the valve, • L is the path length of the valve,

• d = L
2 cosθ is the distance between the centre of the disc and the ellipse one,

• r = D 2 sinθ is the half-small axis of the projected ellipse. Actually the evolution of the opening surface area S Exp is expressed in function of θ, the opening angle of the bushel valve. As the opening law of the valve is provided by the experimental data at LEGI, the time evolution of the experimental opening surface area S Exp (calculated "theoretically") can be deduced.

The second issue is that the inlet mesh of the pipe is discrete. Thus, it is needed to determine which cells are well in the open inlet area of the valve, called S N um . For that, the open inlet area cells are detected thanks to a second reduction: these cells are in the open area if they are a part of an ellipse area whose centre is I, whose half-small axis is [START_REF] Abgrall | How to present pressure oscillations in multi-component flow calculations: A quasi-consevative approach[END_REF] 2 JK and whose half-big axis is IM . Therefore, it amounts to say that the coordinates (Y,Z) of these open inlet area cells verify:

Y -Y I 1 2 JK 2 + Z -Z I IM 2 ≤ 1, (9)
which is then used directly by the code to assess the possible numerical opening surface area S N um . By comparison of the values of S N um and S Exp , it can be deduced the opening times of the valve t i , then the inlet mass flow rate of water corresponding Q Exp (t i ) (which is interpolated in the code) and finally the numerical velocity V N um (t i ) given by:

V N um (t i ) = Q Exp (t i ) S N um (t i )ρ 1 cos(α N um ) , (10)
where

α N um = π 2 -θ ′ is the inlet injection angle, with θ ′ = θ + γ = θ + arctan -1 cos(θ) 1+sin(θ)
(θ and γ are represented in figure 3 (a)). A study (not presented here) was led and reveals that the injection angle α N um of the inlet velocity has few effects on our results. Therefore, α N um was calculated to be the most close to the real opening of the bushel valve.

Outlet boundary condition

At each injector outlet, are imposed a free outlet pressure and a singular head loss ∆P = 1 2 ζρ k U 2 k distributed on the four bottom cells of each injector, with ζ = 1 in order to model the experimental sudden expansion at the injectors outlet.

Outlet boundary condition

Along the walls, standard turbulent friction functions are used for both phases.

IV. Results

The aim of this work is to compare the numerical results obtained with NEPTUNE CFD code with the experimental data of the LEGI, in order to study the filling rate of the cavity, the mass flow rate distribution at injection and the pressure in the cavity. In this section, a reference case is first presented with comparisons to experimental results, before focusing on the influence of the turbulence modeling and on the influence of the drag models, globally and locally.

A. Reference case and compararisons with experimental data

Validation of the valve opening modelised

This first part presents the preliminary validation of the valve opening thanks to comparisons between the inlet data implemented in NEPTUNE CFD code at IMFT and the experimental data given by the LEGI.

The figure 4 compares the experimental opening valve angle θ and the corresponding opening surface area deduced S Exp to the numerical data implemented. As written before, the inlet mesh of the feeding pipe is discrete, that is why the numerical surface area S N um is crenellated. The red dashdotted lines on these graphs represent the different instants of the position of the bushel valve: the first one for the effective beginning of the valve opening (θ = 19.6°at t = 0.123 s), the second one when the valve is completely open (θ = 90°), the third one when the valve begins to close up to θ = 19.6°, which is marked by the last red line. When θ < 19.6°, the valve is effectively closed and there is no inlet mass flow rate of water. Therefore, the numerical simulations really start when θ = 19.6°(at t = 0.123 s) to spare time, whereas the experiments begin before at t = 0 s. So in the following figures, the time origin is t = 0.123 s, so that the figures are more clear. the opening valve and in perfect agreement when the valve is completely open. Nonetheless it can be noticed deviation from the experimental inlet velocity during the closing phase. Indeed, the inlet opening surface area being discrete, gaps of inlet velocity can be observed during this phase.Future improvements will allow S N um to be even closer to S Exp , like this these velocity jumps will be reduced. • : Experimental inlet velocity deduced, : numerical inlet velocity, : valve position.

Points of comparison

To compare the numerical results with the experimental results, it is interesting to set up common points of comparison.

The first one is the pressure measured in the dome. A probe was placed at the level of the top of the dome behind the igniter in NEPTUNE CFD code, like the experimenters at LEGI.

A second point of comparison is the evolution of the fraction of wet injectors number named χ. Experi- mentally, an injector is said to be "wet" if a drop of water is detected at injector outlet by the laser sheet. Numerically this experimental criterion is not valid: it is impossible to detect a drop of water at injector outlet because the volume fraction of water is never totally null during the eulerian simulations. Therefore a numerical "wet" level is equally defined: an injector is defined as "wet" if the mass flow rate of water is bigger than a certain mass flow rate level Q Level , defined as a part of the maximum outlet mass flow rate by injector. Knowing that this mock-up is made up of 122 injectors, the maximum mass flow rate by injector is Q M ax /122, where Q M ax is the maximum mass flow rate of water when the bushel valve is completely open. Thanks to the time evolution of the fraction of wet injectors number χ, it is possible to determine three criteria for χ :

• the date of the first wet injector, named a,

• the date of the beginning of the slope, named b,

• the slope of injection, named c.

As an example, the evolution of the fraction of wet injector number is represented in figure 7, with the three wet criteria a, b and c.

Besides, the LEGI has determined a wet law which permits to have an approximation of the evolution of the fraction of wet injectors number. [START_REF] Cartellier | Activité LOX dôme : investigation expérimentale du remplissage de la cavité aux temps courts aux temps longs[END_REF] This law is given by :

χ(t) = 0 si t ≤ τ 0 inj , (11)
χ(t) = t -τ 0 inj τ inj si τ 0 inj ≤ t ≤ τ 0 inj + τ inj , (12)
χ(t) = 1 si t ≥ τ 0 inj + τ inj , (13)
where: τ 0 inj = 3.4 τ 0 ( T0 τ0 ) 0.75 , τ inj + τ 0 inj = 3.87 τ 0 ( T0 τ0 ) 0.68 , T 0 = L0 V M ax with L 0 the length between the feeding pipe and the first wet injector, and V M ax the supposed maximum inlet velocity of water when the valve is completely open. It has to be noted that in the theorical wet law, the criterion a is the same as the criterion b. That is why, the comparisons to the numerical results will be led only on the criteria b and c.

Results of the reference case

The reference case presented here was performed for a Weber number W e = 4.72.10 4 , that is to say V M ax = 11.5 m/s and Q M ax = 5.65 kg/s. The k-ε turbulence model is activated. The maximum inlet velocity V M ax imposes the turbulent kinetic energy and the turbulent dissipation with a turbulent rate chosen at inlet I = 5%. The drag law chosen is the Simmer-like model with the inclusion diameter d = 10 -3 m.

In figure 8 are represented the inlet and outlet mass flow rates and the mass of water in dome. Until t = 0.823 s, the outlet mass flow rate of water is lower than the inlet one: the mass of water increases in dome and approximately occupies 90% of the dome volume to the maximum. This is the dome filling phase. Then from t = 0.823 s until t = 1 s, the outlet mass flow rate is equal to the inlet one, so the mass of water in the dome is constant. Finally, from t = 1 s, the outlet mass flow rate becomes higher than the inlet one: the mass of water slowly decreases in the dome. This is the emptying dome phase. The experimental and the numerical dome pressures are depicted in figure 9. First it can be observed that the both pressures increase and decrease at the same time. Nonetheless, the numerical plateau value is different from the experimental one: on the one hand, the numerical pressure is higher than the experimantal one and on the other hand, the numerical pressure fluctuates much more.An hypothesis would be that the numerical surface area of water at injectors outlet would differ from the experimental one. A work is in perspective to study the injectors outlet and to try to understand this difference. Moreover simulations will be carried out with a compressible gas. After having studied globally the flow, it would be interesting to locally focus on the evolution of wet injectors to compare the numerical thresholds with the theoretical wet law. The figure 10 shows the fraction of wet injectors number χ for different level mass flow rate of water (from 0.0001% to 100%) and the theorical experimental wet law. Excepted for the threshold at 100%, there is rather good agreement between the experimental wet law and the numerical threshold chosen. Besides, the visualization of all thresholds shows that some injectors have a mass flow rate of water very big: at t = 0.475 s, 20% of injectors have a mass flow rate of water at 100%, whereas the 80% of one remaining have a mass flow rate at 10% or less. There is also an imbalance of the mass flow rate of water at injectors outlet. In order to have more information on the wet injectors, it can be refered to the table 1 which lists the experimental and the numerical values for criteria b and c according to the threshold chosen, and to the associated figure 11. Except the threshold at 1%, the numerical results are in very good agreement for the criterion b, which represents the beginning of the slope at the first order. Nonetheless, the relative deviation for the criteria c, which represents the slope of injection, is very great. It has to be mentioned that the experimental law is a first-order model and that the numerical slope remains subjective. Thus this criterion is not necessarily very appropriate for the comparison. It is more interesting to base oneself on the shape of the curves. For the following comparisons, all the thresholds will be considered.

B. Influence of the turbulence modeling

A step of this work focuses on the influence of the turbulence modeling. More precisely the impact of the k-ε turbulence model is globally and locally analysed, throught three cases:

• one without turbulence modeling,

• another with k-ε turbulence model with a turbulent rate imposed at inlet I = 5%,

• a last one with k-ε turbulence model with a bigger turbulent rate imposed at inlet I = 10%.

The case presented here is performed for a Weber number W e = 4.72.10 4 , that is to say V M ax = 11.5 m/s and Q M ax = 5.65 kg/s. The drag law chosen is the Simmer-like model with the inclusion diameter d = 10 -3 m.

The turbulence modeling and the turbulent rate I have no decisive impact globally on the evolution of the outlet mass flow rate of water (in figure 12(a)). Nonetheless, observing more precisely the mass of water in dome (presented in figure 12(b)) shows a difference of behavior when the valve begins to close. Emptying the dome is a little faster without turbulence modeling. There is also little influence of the turbulence modeling and of the inlet turbulent rate I on the evolution of the dome pressure, depicted in figure 13. Moreover, the turbulence modeling or the turbulent rate I have also few impact on the evolution of the fraction of wet injectors number χ depicted in figure 14 and on the criteria b and c listed in table 2.

In conclusion, activating the turbulence modeling or modifying the turbulent rate I play no macroscopic or local role.

C. Influence of the interfacial momentum transfer modeling

In this part the influence of the interfacial momentum transfer modeling is studied. First two models are tested: the Simmer-like law and the Large Interface Model called LIM. Then the influence of the particule diameter d, parameter of the Simmer-like law, is more precisely studied, the diameter at d = 10 -3 m. The case presented here is performed for a Weber number W e = 4.72.10 4 , that is to say V M ax = 11.5 m/s and Q M ax = 5.65 kg/s. The k-ε turbulence model is activated with the inlet turbulent rate I = 5%.

Influence of the drag law

In this section, the influence of two drag laws are studied: the Simmer-like law (with an inclusion diameter d = 10 -3 m) and the Large Interface Model called LIM.

The evolutions of the outlet mass flow rates of water and the mass of water in dome (in figure 15) have the same behavior until 0.400 s. Then their evolutions differ: more fluctuations of the outlet mass flow rate can be observed with LIM. Moreover, the outlet mass flow of water is bigger with LIM, that is why the mass of water in dome is smaller with LIM. Until t = 0.823 s, the dome is filling before reaching a maximum filling to 80% for LIM and 90% for the Simmer-like model. After t = 0.823 s, the difference of behavior is increased between the two drag models: the dome drains off water much faster with LIM. This different behavior can be explained by a mechanical friction between liquid and gas phases lower with LIM. Indeed the exchange surface area is lower with LIM. Thus, the friction is lower and that is why the outlet mass flow rate is higher with this model. A similar assessment can be made for the pressure in dome, depicted in figure 16. After a similar phase, the two pressures differ: the pressure with LIM is more fluctuating and bigger, because the dome is emptying faster with LIM.

The global effects highlighted before also find oneselves locally in the evolution of the fraction of wet injectors number χ, depicted in figure 17. All the injectors are not wet at the same time with LIM, contrary to the Simmer-like model. Even if the values of the criterion b which represents the date of the beginning of the slope are rather close, the values of the criterion c, the slope of injection, are always lower with LIM, as it is confirmed in the table 3.

To conclude, notable differences of behavior between the two drag models are noticed globally or locally. Thus the drag law model has a major impact on our results. Even if the Simmer-like model seems to be less physical than LIM, it would be prefered. Indeed LIM, based on a free surface detection, is in course of validation at EDF R&D: for the moment, LIM is more adapted to stationnary flows with structured mesh, contrary to our cases whose water / air interface moves fast in an unstructured mesh. After having tested the two drag models (the Simmer-like law or the Large Inferface model), the influence of the particule diameter d, which is a parameter of this model, is studied in this part. Three diameters are tested: 10 -5 , 10 -4 and 10 -3 m. It has to be noted that reducing the inclusion diameter amounts to increase the drag force (as shown in Eq. ( 5), ( 6) and ( 6)).

In figure 18(a) are represented the evolution of the mass flow rates of water. It has to be pointed out that the mean mass flow rate of water is the same for the three cases, let the value of the inclusion diameter. However the behavior differs as regards the mass of water in dome (figure 18(b)): there is no monotonous effect of the particule diameters on the evolution of the mass. Besides it can be noticed that the evolution of the mass of water is smoother with the lower inclusion diameter d = 10 -5 m, because the mechanical friction is bigger. Similarly, the evolution of the dome pressure depicted in figure 19 is not monotonous with the inclusion diameter. When d decreases, the interfacial friction increases and the pressure decreases, except for d = 10 -4 m. Actually, it would be expected that the pressure decreases with the inclusion diameter d. Indeed reducing d shoud homogenize the water / air mixture, because the inclusions are smaller, so the exchange surface area and the friction are higher and the pressure should decrease.A case at d = 10 -6 m would confirm the trend to increase the pressure when the inclusion diameter decreases. The evolution of the fraction of wet injectors number χ depicted in figure 20 and the criteria b and c listed in table 4 show a weak trend of the injectors wetting to arrive earlier when the inclusion diameter deacreases. Besides, more the diameter is low, more the evolution of χ is smooth and more the maximum of χ is reached early. To conclude, the diameter inclusion at d = 10 -3 m seems to be in good agreement with experimental data in our case. : Experimental wet law, : Threshold 0.0001%, : Threshold 0.01%, △ : Threshold 0.1%, : Threshold 1%, • : Threshold 10%, : Threshold 100%, : valve position.

V. Conclusion

The feeding of the LOX dome of a cryogenic rocket-engine is a decisive stage of the transient engine ignition. However flight conditions are difficult to reproduce by experimental ground tests. The work reported here is part of an ongoing research effort to develop a robust method for prediction and understanding the LOX dome feeding. In the framework of this project, experiments with substition fluids (air and water) are conducted, without mass and energy transfer. This work presented here intends to reproduce these experiments through incompressible two-phase flow CFD simulations, in an industrial geometry equivalent to the experimental mock-up, made up of a feeding pipe, a dome with 122 injectors and an igniter.

An important work was made to obtain the same inlet conditions in NEPTUNE CFD code as the experimenters (in particular in the modeling of the valve opening at pipe inlet level), in order to compare for the best the numerical results obtained at IMFT with the experimental results at LEGI. This work focuses also on the influence of the drag models (the Simmer-like law and the Large Interface Model called LIM) and the turbulence modeling. The turbulence modeling or the inlet turbulent rate play no macroscopic or local role on the mass flow rate of water, on the mass of water in dome and in the dome pressure. The drag model has a major impact on our results as well globally as locally, unlike the turbulence modeling. The Simmer-like model is prefered in comparison to the Large Interface called LIM, because it is in better agreement with experimental data. Moreover, it has to be highlighted that the Simmer-like model is very sensitive to its parameter d, the inclusion diameter. In our case, the diameter inclusion at d = 10 -3 m seems to be in good agreement with experimental data.

A work on the inlet boundary condition is in sight, so that the numerical inlet opening surface area S N um is even closer to the experimental one S Exp theoretically calculated. This will improve the numerical inlet velocity imposed in NEPTUNE CFD code. Besides, further work will include a more specific study of the impact of the outlet boundary condition, with the addition of a tank at injector outlet, which would be more representative of the experiments at LEGI. In addition, simulations will be carried out with a compressible gas.

oriented fractional step method, Third International Symposium on Finite Volumes for Complex Applications, Porquerolles, France, 2002.

13 N. Méchitoua, M. Boucker, J. Laviéville, J.-M. Hérard, S. Pigny and G. Serre. An unstructured finite volume solver for two-phase water / vapour flows modeling based on ellipstic oriented fractional step method, NURETH 10, Seoul, South Corea, 2003

Figure 1 :

 1 Figure 1: Experimental mock-up at LEGI.

  (a) The complete mesh of the mock-up.(b) Zoom on the feeding pipe -dome joint.

Figure 2 :

 2 Figure 2: Three-dimensional mesh for the simulations.

  The simplified 2D approach.

Figure 3 :

 3 Figure 3: The 3D geometric problem of the valve turned into a 2D problem in the same plane.

2 ]

 2 Opening surface area of the valve [m (b) Experimental and numerical surface areas of the valve.

Figure 4 :

 4 Figure 4: Comparisons between the experimental inlet data and the numerical input data for the opening valve angle and surface area. • : Experimental inlet data, : numerical inlet data, : valve position.

Figure 5 :

 5 Figure 5: Comparison between the experimental inlet data and the numerical input data for the mass flow rate of water. • : Experimental inlet data, : numerical inlet data, : valve position.

  Inlet velocity of water [m/s]

Figure 6 :

 6 Figure 6: Comparison between the experimental inlet velocity deduced and the numerical inlet velocity.• : Experimental inlet velocity deduced, : numerical inlet velocity, : valve position.

Figure 7 :

 7 Figure 7: Representation of the three wet criteria (a, b and c) on the evolution curve of the fraction of wet injectors number.

  Mass of water in dome [kg] and mass flow rate of water [kg/s]

Figure 8 :

 8 Figure 8: Numerical mass of water and mass flow rates of water. • : Inlet mass flow rate, : outlet mass flow rate, △ : mass of water in dome, : valve position.

Figure 9 :

 9 Figure 9: Comparison of the dome pressure.: Experimental pressure, : numerical pressure, : valve position.

Figure 10 :

 10 Figure 10: Comparison of the fraction of wet injectors number χ.: Experimental wet law, : Threshold 0.0001%, : Threshold 0.01%, △ : Threshold 0.1%, : Threshold 1%, • : Threshold 10%, : Threshold 100%, : valve position.

Figure 11 :

 11 Figure 11: Influence of the threshold on the criteria b and c of the fraction of wet injectors number χ. : Experimental values, ⋄ : numerical criterion b, * : numerical criterion c.

Figure 12 :

 12 Figure 12: Influence of the turbulence modeling on the evolutions of the mass flow rate of water and the mass of water in dome. : No turbulence modeling, : k-ε turbulence model with I = 5%, : k-ε turbulence model with I = 10%, : valve position.

Figure 13 :

 13 Figure 13: Influence of the turbulence modeling on the evolution of the dome pressure. : Experimental pressure, : No turbulence modeling, : k-ε turbulence model with I = 5%, : k-ε turbulence model with I = 10%, : valve position.

  injectors number χ (b) k-ε turbulence model with I = 5%. injectors number χ (c) k-ε turbulence model with I = 10%.

Figure 14 :

 14 Figure 14: Influence of the turbulence modeling on the evolution of the fraction of wet injectors number χ.: Experimental wet law, : Threshold 0.0001%, : Threshold 0.01%, △ : Threshold 0.1%, : Threshold 1%, • : Threshold 10%, : Threshold 100%, : valve position.

  Mass flow rate of water [kg/s] (a) Evolution of the mass flow rate of water. • : Inlet mass flow rate, or : Outlet mass flow. Mass of water in dome [kg](b) Evolution of the mass of water in dome.

Figure 15 :

 15 Figure 15: Influence of the drag law on the evolution of the mass flow rate of water and the mass of water in dome. : Simmer-like model (d = 10 -3 m), : Large Interface Model, : valve position.

Figure 16 :

 16 Figure 16: Influence of the drag law on the evolution of the dome pressure. : Experimental pressure, : Simmer-like model (d = 10 -3 m), : Large Interface Model, : valve position.

Figure 17 :

 17 Figure 17: Influence of the drag law on the evolution of the fraction of wet injectors number χ. : Experimental wet law, : Threshold 0.0001%, : Threshold 0.01%, △ : Threshold 0.1%, : Threshold 1%, • : Threshold 10%, : Threshold 100%, : valve position.

  Mass flow rate of water [kg/s] (a) Evolution of the mass flow rate of water. Mass of water in dome [kg] (b) Evolution of the mass of water in dome.

Figure 18 :

 18 Figure 18: Influence of the inclusion diameter in the Simmer-like law on the evolution of the mass flow rate of water and the mass of water in dome. : d = 10 -5 m, : d = 10 -4 m, : d = 10 -3 m, : valve position.

Figure 19 :

 19 Figure 19: Influence of the inclusion diameter in the Simmer-like law on the evolution of the dome pressure. • : Experimental pressure, : d = 10 -5 m, : d = 10 -4 m, : d = 10 -3 m, : valve position.

Table 4 :

 4 Influence of the inclusion diameter in the Simmer-like law on the criteria b and c, according to the threshold chosen. Threshold Criteria d = 10 -5 m d = 10 -4 m d = 10

  injectors number χ (c) d = 10 -5 m.

Figure 20 :

 20 Figure 20: Influence of the inclusion diameter in the Simmer-like law on the evolution of the fraction of wet injectors number χ.: Experimental wet law, : Threshold 0.0001%, : Threshold 0.01%, △ : Threshold 0.1%, : Threshold 1%, • : Threshold 10%, : Threshold 100%, : valve position.

Table 1 :

 1 Influence of the threshold chosen on the criteria b and c.

	Threshold Criteria Numerical results Experimental data Relative deviation [%]
	0.0001%	b [s]	0.253	0.266	4.9
		c [χ/s]	40.8	28.6	-42.7
	0.01%	b [s]	0.268	0.266	-0.8
		c [χ/s]	108.5	28.6	-279
	0.1%	b [s]	0.271	0.266	-1.9
		c [χ/s]	127.7	28.6	-347
	1%	b [s]	0.150	0.266	43.6
		c [χ/s]	126.8	28.6	-343
	10%	b [s]	0.273	0.266	-2.6
		c [χ/s]	71.7	28.6	-151

Table 2 :

 2 Influence of the turbulence modeling on the criteria b and c, according to the threshold chosen.

	Threshold Criteria No turbulence	k-ε turbulence	k-ε turbulence	Experimental
			modeling	model with I = 5% model with I = 10%	data
	0.0001%	b [s]	0.254	0.253	0.256	0.266
		c [χ/s]	42.5	40.8	46.3	28.6
	0.01%	b [s]	0.262	0.268	0.264	0.266
		c [χ/s]	62.5	108.5	60.1	28.6
	0.1%	b [s]	0.266	0.271	0.271	0.266
		c [χ/s]	80.8	127.7	104.2	28.6
	1%	b [s]	0.148	0.150	0.151	0.266
		c [χ/s]	10.7	126.8	115.3	28.6
	10%	b [s]	0.273	0.273	0.276	0.266
		c [χ/s]	85.8	71.7	93.4	28.6

Table 3 :

 3 Influence of the drag law on the criteria b and c, according to the threshold chosen.

	Threshold Criteria Simmer-like law Large Interface Model Experimental data
	0.0001%	b [s]	0.253	0.248	0.266
		c [χ/s]	40.8	35.7	28.6
	0.01%	b [s]	0.268	0.254	0.266
		c [χ/s]	108.5	38.5	28.6
	0.1%	b [s]	0.271	0.261	0.266
		c [χ/s]	127.7	51.3	28.6
	1%	b [s]	0.150	0.141	0.266
		c [χ/s]	126.8	51.8	28.6
	10%	b [s]	0.273	0.271	0.266
		c [χ/s]	71.7	66.4	28.6
	2. Influence of the inclusion diameter in the Simmer-like law

Acknowledgments

The authors would like to acknowledge Snecma ans CNES for funding this study through a PhD CNES program and the scientific support of P. Boivin (Snecma), J. Laviéville and N. Méchitoua (EDF R&D). The simulations were performed at the Hight Performance Computing centre CALMIP (CAlcul en MIdi Pyrénées) under the project P0111.