L Germes Martinez

B Duret
email: duret@coria.fr

J Reveillon

F X Demoulin

Vapor mixing in turbulent vaporizing ows

Keywords: DNS, Compressible, Turbulent mixing, Two-phase ows, Evaporation, CLSVOF

Scalar mixing and evaporation processes are analyzed in a compressible twophase homogeneous and isotropic turbulence conguration. To this end, the fully compressible Navier-Stokes equations are solved using a pressure-based method and a mass-conservative interface capturing method (CLSVOF). Temperature and vapor mass fraction transport equations are coupled with the momentum equation via the evaporation rate: i.e. heat and mass transfer has an eect on the ow dynamic through the generation of a Stefan ow at the interface. In this formalism, a pressure equation is developed to include the acoustic, thermal dilatation and compressibility eects. Moreover, the method can consider several gas structures in the same domain, each with independent temperature, density, and thermodynamic pressure. The latter is important to simulate atomization, where interaction among the liquid structures (such as breakup or coalescence) can generate multiple gas inclusions. This work is one of the rst to describe with high delity the vaporization and turbulent mixing occurring in dense two-phase ows.

A statistical analysis of the vapor mass fraction eld is performed. The obtained results reveal the inuence of the non-homogeneous scalar source at the interface, which depend directly on the local interface temperature, on the PDF shape. Under these conditions, the vapor mass fraction PDF deviates from the Gaussian and beta PDF shapes, frequently used for turbulent combustion modeling. Furthermore, access to the local surface averaged vapor mass fraction and

Introduction

Mixing and evaporation have a signicant role in many natural and engineering processes. For instance, the performance and pollutant emissions of internal combustion engines depend highly on these two processes. Generally, mixing is enhanced by a turbulent ow, when compared with a process driven mainly by diusion. The ability of these ows to increase the surface between the uids intensies the molecular diusion [START_REF] Dimotakis | Turbulent mixing[END_REF]. An example of the importance of turbulent mixing in nature is the release of pollutants (such as chemicals, solid particles, or heat) into the atmosphere or aqueous environments.

Furthermore, in industry, turbulent mixing is the driving mechanism for heat and mass transfer in a chemical reactor, especially in cases where the turbulent and reaction time scales are comparable. During the mixing process, three distinct stages can be identied [START_REF] Eckart | An analysis of the string and mixing processes in incompressible uids[END_REF]. First is the initial state, in which the entrainment of one uid into another occurs. This is characterized by large concentration gradients at the interface and a near-zero gradient elsewhere. In the intermediate state, the motion created by the larger eddies distorts the interface and increases molecular transport between the uids. As the uids become more homogeneous at the nal state, the concentration gradients decrease until there is no more heat or mass transfer.

Normally, scalar mixing in turbulent ows is studied with dierent strategies depending on how the scalar is released into the carrier phase [START_REF] Hunt | Turbulent diusion from sources in complex ows[END_REF]. For instance, in experiments the scalar is released as a single point source in the ow by adding a relatively small amount of dye to a liquid. Such scenarios have been studied in an experimental framework using the evolution of single-point and two-point statistics of the scalar elds.

Experimental measurements of the probability density function (PDF) of a passive scalar in a homogeneous turbulent shear ow performed by [START_REF] Tavoularis | Experiments in nearly homogenous turbulent shear ow with a uniform mean temperature gradient. part 1[END_REF] indicated a quasi-normal distribution for the scalar signal.

Subsequently, several experimental and theoretical works show that, under certain conditions, the scalar PDF, with a mean scalar gradient, exhibits an exponential tail due to anomalous mixing [START_REF] Shraiman | Lagrangian path integrals and uctuations in random ow[END_REF]; [START_REF] Chertkov | Intermittent dissipation of a passive scalar in turbulence[END_REF]; [START_REF] Pumir | Exponential tails and random advection[END_REF]). Moreover, [START_REF] Jayesh | Probability distribution, conditional dissipation, and transport of passive temperature uctuations in grid-generated turbulence[END_REF] studied the evolution of the scalar PDF and the dissipation rate for passive temperature uctuations in decaying grid-generated turbulence; concluding that the scalar PDF has a pronounced exponential tail for turbulence Reynolds number (Re λ ) above 70 but below this value the PDFs can be modeled by a normal distribution.

Several authors have also resorted to direct numerical simulations (DNS) to investigate passive scalar mixing in turbulent ows. For example, [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF] used DNS to study the mixing of a passive scalar in an isotropic turbulent ow. Here, the authors observed the eect of the initial conditions on the temporal evolution of the scalar rms and the scalar dissipation rate. Additionally, they analyzed the evolution of the scalar PDF shape from a double-Dirac function distribution to a Gaussian distribution. Similar simulations were carry out by [START_REF] Overholt | Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence[END_REF]; [START_REF] Juneja | A dns study of turbulent mixing of two passive scalars[END_REF]; [START_REF] Yeung | Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations[END_REF]. Furthermore, [START_REF] Yeung | Schmidt number eects on turbulent transport with uniform mean scalar gradient[END_REF] used DNS of isotropic and homogeneous turbulence to investigate the eects of the Schmidt number on the inertia and dissipation length scales of the scalar energy spectrum. More recently, [START_REF] Yeung | High-reynolds-number simulation of turbulent mixing[END_REF] studied a 2048 3 DNS of passive scalar mixing in a turbulent ow, observing a deviation from small-scale isotropy. [START_REF] Yeung | Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations[END_REF]Sreenivasan (2013, 2014) used DNS to study the mixing of passive scalars with a very low Schmidt number; the goal was to conrm the existence of the -17/3 power law of the scalar spectral density in the inertial-diusive range.

In turbulent combustion modeling, one way for approximating the local averaged chemical production rate is the use of a presumed scalar PDF. In this method, the scalar pdf shape is assumed a priori, and the parameters of such PDF are determined by solving the transport equation of the rst couple of moments [START_REF] Girimaji | Assumed β-pdf model for turbulent mixing: Validation and extension to multiple scalar mixing[END_REF]). The PDFs commonly used for this purpose are the Gaussian [START_REF] Lockwood | The prediction of the uctuations in the properties of free, round-jet, turbulent, diusion ames[END_REF]) and the beta-pdf [START_REF] Gerlinger | Investigation of an assumed pdf approach for nite-rate chemistry[END_REF]; [START_REF] Lempke | Inuence of spatial discretization and unsteadiness on the simulation of rocket combustors[END_REF]). An example is the work of [START_REF] Gerlinger | Investigation of an assumed pdf approach for nite-rate chemistry[END_REF], which performs simulations of a supersonic hydrogen-air diusion ame assuming the multivariable beta-pdf proposed by [START_REF] Girimaji | Assumed β-pdf model for turbulent mixing: Validation and extension to multiple scalar mixing[END_REF]. A major drawback of this approach is that the PDF shape is chosen from a mathematical perspective instead of a physical one. To improve this method, [START_REF] Bray | Finite rate chemistry and presumed pdf models for premixed turbulent combustion[END_REF] investigated the accuracy of predicting the mean reaction rate in turbulent premixed ames using three presumed pdf: beta-pdf, twin delta function, and a pdf based on unstrained ame properties. They found a large deviation with the rst two assumed PDFs compared with DNS calculations. Similar results were found by [START_REF] Salehi | Presumed pdf modeling for rans simulation of turbulent premixed ames[END_REF] when simulating a premixed ame with an assumed beta-pdf approximation and a modied laminar amelet PDF.

More recently, a new approach was proposed by de [START_REF] De Frahan | Deep learning for presumed probability density function models[END_REF].

Machine learning techniques were used to develop assumed PDF models for reacting ows. When comparing the proposed method with a standard beta pdf model, signicant improvements were found. Furthermore, [START_REF] Lipatnikov | Assessment of a amelet approach to evaluating mean species mass fractions in moderately and highly turbulent premixed ames[END_REF] used the PDF directly extracted from DNS and presumed beta-function using the rst two moments yielded by DNS data for modeling the combustion process variable and the mean source term in a premixed turbulent ame.

The studies mentioned above on passive scalar mixing in a turbulent ow use a single-phase conguration. In most cases, the scalar source is dened during the simulation initialization as one or several vapor patches arbitrarily distributed in the domain. However, scalar mixing has been studied using a Lagrangian approach to account for the two-phase ows. An example is the work of [START_REF] Reveillon | Eects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF], where the eect of the droplet's spatial distribution on evaporation and scalar mixing is investigated. They use a spectral DNS solver for the evolution of the turbulent carrier phase and a Lagrangian model for the dispersed phase, which acts as scalar point sources. Several publications focus on the study of turbulent mixing in similar scenarios [START_REF] Reveillon | Dns study of spray vaporization and turbulent micro-mixing[END_REF]; [START_REF] Demoulin | Modeling of turbulent spray combustion with application to diesel like experiment[END_REF]; [START_REF] Colin | A new scalar uctuation model to predict mixing in evaporating two-phase ows[END_REF]). The main interest of this type of work is to understand the eect of the vaporization source term in the scalar transport equation on the evolution of the mean vapor mass fraction, uctuations and dissipation rate.

A limitation when using the point-particle Lagrangian approximation is that the simulations are restricted to particles smaller than the Kolmogorov scale.

A fully resolved DNS coupled with an interface capturing method must be performed to simulate larger particles. This type of conguration was adopted by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]. They studied the mixing process in a two-phase HIT conguration with a vaporizing liquid. Statistical and spectral analysis were presented to demonstrate the eects of liquid volume fraction on the velocity and scalar elds. Nevertheless, this work does not consider the coupling between vaporization rate and ow dynamics through the velocity jump condition at the interface. Consequently, the study is limited to low vaporization rates, where the phase change process hardly aects the interface and ow dynamics. [START_REF] Bouali | Dns analysis of smallscale turbulence-scalar interactions in evaporating two-phase ows[END_REF] used a similar setup to study the inuence of the liquid/gas interfaces on the scalar dissipation rate statistics.

In the literature, several numerical strategies include the relation between temperature, vapor mass fraction and uid dynamics in their formalism (Es-maeeli and Tryggvason, 2004;[START_REF] Tanguy | A level set method for vaporizing two-phase ows[END_REF][START_REF] Malan | Direct numerical simulation of free-surface and interfacial ow using the VOF method: cavitating bubble clouds and phase change[END_REF][START_REF] Palmore | A volume of uid framework for interfaceresolved simulations of vaporizing liquid-gas ows[END_REF][START_REF] Scapin | A volume-of-uid method for interfaceresolved simulations of phase-changing two-uid ows[END_REF]. Since the conservation equations are coupled through the interface jump conditions, the above methods can be subdivided according to the interface representation used. To implement one of these methods to study the scalar mixing and vaporization processes in a turbulent and conned environment, the eect of compressibility must be considered in order to conserve the initial total mass. A step toward this kind of simulation was taken by [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] and [START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF], who used a low-mach formulation of the Navier-stokes equations. Consequently, the gas density is a function of pressure and temperature in their simulations. However, this formalism considers a spatially constant thermodynamic pressure, which could aect the result if more than one gas structure is present in the domain. To overcome this limitation, a fully compressible formulation of the Navier-Stokes equations, where the acoustic and compressibility eect are resolved, could be used. For instance, the method proposed by [START_REF] Urbano | A semi implicit compressible solver for two-phase ows of real uids[END_REF] is a good candidate for this type of simulation. They use a semi-implicit projection method to solve the compressible Navier-Stokes equations coupled with cubic equations of state, which allows them to consider real uid interactions. However, mass transfer is not yet implemented in their formalism.

In this work, the all-Mach formalism proposed by [START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF] is used to study the scalar mixing and evaporation process. To this end, a twophase HIT conguration with phase change is considered. In this formalism, the interface velocity jump is handled explicitly during the resolution of the compressible Navier-Stokes equations, and the interface regression induced by the liquid evaporation is taken into account using a sink term in the VOF transport equation. Consequently, the Stefan ow and its eect on the system dynamics are captured. Moreover, the momentum, energy, and species conservation equations are coupled by the evaporation rate present in the jump conditions, resulting in a non-homogeneous temperature and vapor mass fraction at the interface. Another consequence of the coupling among the conservation equations is that the scalar is not passive since its eventual variation will directly aect the turbulent velocity eld. For the reasons mentioned above, the description of the physical phenomena in our work could be considered a major step forward compared with the available publications in the literature. Moreover, this study is considered as an important evolution of our previous study which used a passive scalar to model vaporization [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]. Such improvement of the description of phase change has been achieved with the development of numerical methods pursued in the last decade in our team.

In the following, the governing equations, jump conditions, and numerical methods are rst described. The numerical conguration is then explained alongside the initial simulation parameters. Then, a statistical analysis of the vapor mass fraction and evaporation rate is performed based on the temporal evolution of the mean, variance and PDFs of both variables. Ultimately, the role of the interface curvature on the vaporization rate magnitude is discussed by analyzing their joint PDF.

Governing equations

This work considers a two-phase ow composed of a monocomponent liquid and a mixture of ideal gases. Both phases are fully compressible and are separated by an innitely thin interface, where heat and mass transfer occurs. The gas phase is initially described as dry air. When evaporation starts, the inert air mixes with the vapor generated around the interface. The system of equations solved for both phases is the following:

                         ∂ρ ∂t + ∇ • (ρu) = 0 ∂ρu ∂t + ∇ • (ρu ⊗ u) = ∇ • Ω + σHδ Γ n + ρf vol ρc p ∂T ∂t + u • ∇T = ∇ • (λ∇T ) + α T T DP Dt ∂Y v ∂t + u • ∇Y v = ∇ • (ρ g D m ∇Y v ) ρ g (1a) (1b) (1c) (1d)
where

D[•] Dt = ∂[•] ∂t + u • ∇[•]
is the material derivative operator. ρ, u, and P are the density, velocity, and pressure, respectively. Ω = -P Ī + τ is the total stress tensor and Ī the identity matrix. The viscous stress tensor is dened as:

τ = 2µ Ξ - 2 3 µ (∇ • u) Ī
where Ξ = 1 2 ∇u + ∇u T and µ the dynamic viscosity. σ is the surface tension, H is the mean curvature, and the surface tension forces are considered based on the Dirac function (δ Γ ). The mean turbulent kinetic energy is maintained near a prescribed value by a forcing method similar to [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]. It consist in employing the volume forces (f vol ) as a source term dened as

f vol = Au
where A is a forcing coecient and u represent the velocity uctuations.

T is the temperature, and Y v is the vapor mass fraction, only dened in the gas phase because the liquid is considered a monocomponent uid. Thermophysical properties, such as the thermal conductivity (λ), thermal dilatation coecient (α T ), liquid heat capacities (c pl ), and mass diusivity (D m ) are considered constants. However, the gas heat capacity (c pg ) is an averaged weighted by the vapor and air mass fractions:

c pg = k=a,v Y k c pk
Since both phases are considered compressible, equations of state for each phase are needed to close the system 1. For the liquid, a modied version of Tait's equation is chosen:

ρ l = ρ 0 P -P 0 B + 1 1 γ l (2)
where ρ 0 is a density at a reference state, P 0 is a reference pressure, B is a constant, which depends on the isothermal compressibility of the uid, and γ l is a parameter inherent to the uid. In the gas phase, the ideal gas equation is used:

ρ g = P M g RT ( 3
)
where R is the ideal gas constant and the gas molar mass (M g ) dened as:

M g =   k=a,v Y k M k   -1
.

Similar to [START_REF] Caltagirone | A multiphase compressible model for the simulation of multiphase ows[END_REF]; [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF]; [START_REF] Duret | A pressure based method for vaporizing compressible two-phase ows with interface capturing approach[END_REF], the following equation is implemented for the spatial and temporal evolution of the pressure:

DP Dt = -ρc 2 ∇ • u + α T c 2 c p (∇ • (λ∇T )) ( 4
)
where c is the speed of sound. More detail about de development of eq. 4 can be found in Germes [START_REF] Martinez | Direct Numerical Simulation of turbulent liquid/gas compressible ows with phase change[END_REF]. To consider the phase change in our system, proper jump conditions at the interface are included, ensuring the conservation of mass, momentum, and species in the whole domain. The additional jump conditions due to phase change are [START_REF] Ishii | Thermo-uid dynamics of two-phase ow[END_REF][START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF]:

[u] = ω 1 ρ Γ P -n • 2µ Ξ • n Γ = σH -ω2 1 ρ Γ [λ∇T ] Γ = -ω h lv ρ g D m ∇Y v • n| Γ g = ω(1 -Y vs ) (5) (6) (7) (8)
where [A] Γ = A l -A g for a given variable A, ω is the vaporization rate and h lv the latent heat.

Numerical methods

In this work, an interface capturing method is coupled with a projection method to solve the system 1. A detailed explanation of the numerical methodology can be found in [START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF]. This section presents a general overview of the main numerical methods used for our simulations.

Interface representation and projection method

The coupled level-set/volume of the uid (CLSVOF) method for the interface representation proposed by [START_REF] Ménard | Coupling level set/VOF/ghost uid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet[END_REF] was modied to handle compressible uids and phase change in the ARCHER code. To this end, two additional terms are added to the RHS of the liquid volume fraction (α l ) transport equation (see [START_REF] Duret | A pressure based method for vaporizing compressible two-phase ows with interface capturing approach[END_REF] or [START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF] for more details about the derivation) :

∂α l ∂t + ∇ • (α l u l ) = - α l ρ l Dρ l Dt - ṁ ρ l (9)
where u l is the liquid velocity and ṁ = ρ l Σ s d is the volumetric evaporation rate. ṁ is estimated using the local surface density

Σ = SΓ V cell
where S Γ is the local surface area computed after the PLIC reconstruction and V cell is the volume of the grid cell. Additionally, the level set function (φ) is used to compute the geometrical properties of the interface, such as, the normals n = ∇φ |∇φ| and the curvature

(κ = -∇ • n).
This work uses the projection method based on two velocities presented in [START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF] to consider compressibility eects and phase change. An Eulerian staggered grid is considered, where the velocities are dened on the faces and the scalars (e.g., P, α l , T, Y v ) are located in the center of the cell. The two velocity elds obtained with this method are the liquid velocity eld and its extension in the gas phase and the gas velocity eld and its extension in the liquid phase, which are continuous at the interface. Then, a Ghost Fluid method [START_REF] Fedkiw | The ghost uid method for deagration and detonation discontinuities[END_REF] is applied to consider the velocity jump at the interface and couple the two velocity elds.

Temperature and vapor mass fraction transport

The temporal evolution of the temperature eld is obtained by solving eq. 1c. Here, a 5th-order WENO scheme is implemented to discretize the convective term, and a 2nd-order central dierence scheme is used for the diusion term.

Moreover, the energy jump condition at the interface is handled using a Ghost Fluid method. To avoid numerical heating/cooling during the solution of the convective term of eq. 1c or the interface transport, the temperature at each phase is extrapolated linearly using Aslam's extension method [START_REF] Aslam | A partial dierential equation approach to multidimensional extrapolation[END_REF]).

Concerning the vapor mass fraction, eq. 1d is solved in the gas phase. The diusive and convective term of eq. 1d are discretized using the same numerical schemes used for the temperature. At the interface, a Dirichlet boundary condition is implemented using the Aslam-Chiu extension method described in [START_REF] Bouali | Dns analysis of smallscale turbulence-scalar interactions in evaporating two-phase ows[END_REF]. To compute the vapor mass fraction at the interface (Y vs ), thermodynamic equilibrium is assumed. The saturation pressure (P sat ) is computed using the Clausius-Clapeyron relation:

P sat = P ref exp - h lv M v R 1 T Γ - 1 T ref ( 10
)
where P ref is the saturation pressure at T ref , M v is the molar mass of the vapor and T Γ is the interface temperature. Then, Y vs is estimated with the following relation:

Y vs = XM v XM v + (1 -X) M a (11)
where X = Psat P is the vapor molar fraction and P is the mean pressure of the system.

Finally, the evaporation rate is computed from eq. 8, using the normal gradient of the vapor mass fraction at the interface:

ω = ρ g D m ∇Y v • n| Γ g (1 -Y vs ) (12)

Surface averaged quantities

The following sections investigate surface averaged variables to describe the vaporization process. The numerical methodology used to obtain such quantities is detailed in this section.

First, the interface is reconstructed by discretizing the zero level-set surface with a triangulated mesh using a Marching Cube algorithm. Then, its geometrical properties, such as surface area and mean curvature, are extracted. To this end, we use the "surf ace_operators" routines of the in-house PyArcher code for post-processing. These routines are based on earlier versions of the routines described by [START_REF] Mohamed | Statistical modeling of the gasliquid interface using geometrical variables: toward a unied description of the disperse and separated phase ows[END_REF] and Di Battista et al. ( 2019) available through the project Mercur(v)e1 .

Then, interpolation schemes of order 0 and 1 are implemented to estimate the evaporation rate and vapor mass fraction at the zero level-set surface, respectively. Finally, the surface averaged quantities are calculated as follows:

ψ Γ = 1 S r ψdS r ( 13
)
where ψ is a scalar (here ω or Y v ), dS r is the local surface in a given cell and S r is the total surface area reconstructed from the zero level-set surface.

Numerical conguration

The well-known two-phase homogeneous isotropic turbulence [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]; [START_REF] Duret | Simulation numérique directe des écoulements liquide-gaz avec évaporation: application à l'atomisation[END_REF]; [START_REF] Duret | A pressure based method for vaporizing compressible two-phase ows with interface capturing approach[END_REF]) is retained to study mixing and vaporization. However, heat and mass transfer was not considered in previous works. Here, the uid thermophysical properties, initial temperature of each phase, and vapor concentration in the box must be specied carefully: they will dictate the vaporization regime observed until a saturation state is reached.

Consequently, the following section aims to estimate the equilibrium state that should be reached with a given set of physical parameters and initial temperatures. Based on this equilibrium, two dierent regimes of vaporization and mixing are chosen and investigated. The equilibrium is deduced from an iterative algorithm based on the rst law of thermodynamics.

Equilibrium Conditions

Phase equilibrium is achieved when there is no mass transfer between the phases [START_REF] Cengel | Thermodynamics: an engineering approach[END_REF]. In other terms, the system has reached phase equilibrium when there is no vapor mass fraction gradient, which is the driving force for mass transfer. Since the vapor mass fraction at the interface is a function of the temperature, thermal equilibrium is implied when the system reaches phase equilibrium. Therefore, in this work, the word equilibrium refers to phase equilibrium.

An energy balance is performed to estimate the system's conditions at the equilibrium state. By considering a closed system where three components (air, vapor, and liquid) coexist at the same time, we can write:

∆U l + ∆U a + ∆U v = 0 (14)
where

∆U l = m l,eq c vl T eq -m 0 l c vl T 0 l -m l,eq -m 0 l h lv ∆U a = m a,eq c va T eq -m 0 a c va T 0 a ∆U v = m v,eq c vv T eq -m 0 v c vv T 0 v
represent the internal energy change between the initial and the nal state of the liquid, air, and vapor phase, respectively. Here, m 0 i , m i,eq and T 0 i are the initial mass, the mass at the equilibrium and the initial temperature of the component i = l, a, v. An expression for the equilibrium temperature (T eq ) is determined from eq. 14, considering the following assumptions: Heat capacities of each component are constant.

Air is not soluble in the liquid.

Initial temperature of the air and vapor are the same.

There is no chemical reaction.

T eq = m 0 l c pl T 0 l + m 0 a c pa T 0 a + m 0 v c p,v T 0 g -m evap h lv m l,eq c pl + m 0 a c pa + m v,eq c pv (15)
The amount of evaporated mass can be determined by m evap = m v,eqm 0

v . An iterative algorithm is used to solve eq. 15 with the condition that the vapor mass fraction calculated with the equilibrium temperature is the same as the vapor mass fraction calculated with the equilibrium vapor mass.

Parameters choice

Similar to previous works [START_REF] Duret | A pressure based method for vaporizing compressible two-phase ows with interface capturing approach[END_REF]; [START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase ows with phase change[END_REF]), a compressible two-phase HIT conguration with phase change is considered to study the scalar turbulent mixing and evaporation. Simulations are carried out in a 3D periodic box with a length size l x = 1.5 × 10 -4 m, a liquid volume fraction α l = 10% and a mesh of 256 3 . The mean turbulent kinetic energy is k = 3.6 m 2 • s -2 , the liquid has the physical properties of n-decane, and the gas is initially considered dry air. Additionally, once evaporation starts, the gas phase is considered an ideal mixture of air and n-decane vapor. The physical properties at the initial state of the simulation are summarized in the table 1. The parameters for the Tait's equation are B = 10 9 P a, P 0 = 10 5 P a, ρ 0 = 750 kg•m -3 and γ l = 1.215. Moreover, the initial values of the simulation's dimensionless numbers are listed in table 2. Two simulations are performed to investigate the inuence of the evaporation regime on scalar turbulent mixing, where only the initial liquid temperature is varied. The initial liquid temperatures of the simulations are T 0 l = 540 K, 460 K and the initial gas temperature in both cases is T 0 g = 573 K. The rst simulation is expected to present a larger evaporation rate than the second since the liquid temperature is closer to the boiling point. Hence, the equilibrium state and the temporal evolution of the temperature and vapor mass fraction are expected to dier in both simulations. Using the method presented in section 4.1, the variables at the equilibrium are estimated for each case, and the results are summarized in the table 3.

Cases T 0 l (K) T 0 g (K) Teq(K) Yv,
In the hypothetical case where evaporation is not considered in the system (only heating), the equilibrium temperature is between the initial liquid and gas temperatures. The table 3 shows that the equilibrium temperature for case 1 is lower than the initial liquid temperature. This is explained by the large evaporation rate, where the liquid's latent heat consumes an important part of the initial energy of the system. In contrast, dierent results are observed for case 2, where the equilibrium temperature is higher than the initial liquid temperature. Since case 2 has a lower evaporation rate, the system's energy in the initial state is sucient to evaporate the liquid and increase its temperature.

However, the increase in liquid temperature is lower than in the hypothetical case without evaporation because some energy is still consumed. In both cases, the equilibrium temperature is smaller than the gas temperature.

Moreover, the impact of the evaporation regime can be seen in the amount of liquid evaporated. In case 1, the equilibrium vapor mass fraction is 2.3 times larger compared to case 2. This indicates a more signicant amount of evaporated liquid in case 1, caused by the increased evaporation rate. Another parameter is the liquid mass at equilibrium normalized by its initial value. Again, a lower liquid mass is observed for case 1, which agrees with the aforementioned discussion.

The variables at the equilibrium are used to normalize the temperature (T * = T Teq ) and the vapor mass fraction (Y * v = Yv Yv,eq ). In addition, to analyze dierent stages of the turbulent mixing, four times are dened (t 1 , t 2 , t 3 and t 4 ), each corresponding to a time where a specic value of the normalized mean vapor mass fraction is reached (Y v * = 0.18, 0.32, 0.54 and 0.72). It is important to note that the physical time represented by each dened time varies depending on the simulation. For example, since the vapor production rate is higher in case 1, t 1 represents a shorter time compared to case 2. Furthermore, the values of Y v * were chosen considering the time required to reach the equilibrium state and the computational cost involved. Nevertheless, at least 20 eddy turnover times are considered for each simulation to have a well-developed mixture of the vapor mass fraction.

Initialization procedure

The simulation initialization is similar to [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]; Bouali et al.

(2016) for the velocity elds. Regarding heat and mass transfer, the temperature is constant in each phase during the initialization procedure, meaning that no diusion of species or heat occurs during this stage. Then, sub-iterations of pure diusion are performed, where the following equations are solved:

T k+1 = T k + ∆t dif f ρ k c k p ∇ • λ∇T k Y k+1 v = Y k v + ∆t dif f ∇ • ρ k g D m ∇Y k v ρ k g (16a) (16b)
where k is the current sub-iteration. The number of diusive sub-iterations (k end ) is computed by:

k end = t BL ∆t dif f ( 17
)
where t BL = δ Y v s dif f is the time to generate a vapor boundary layer close to the interface. More precisely, this procedure move a vapor mass fraction iso-contour to an imposed distance δ Y v = 6dx with an estimated speed of

s dif f = Dm dx . ∆t dif f = 0.2
V cf l is the time step computed with the thermal and mass diusion restrictions:

D cf l = 6Dif f dx 2
where Dif f = max λg ρgcpg , λ l ρ l c pl , D m . This procedure avoids numerical instabilities in the vaporization rate computation created by high temperature and vapor mass fraction gradients in the initial elds. Finally, the energy and vapor mass fraction transport equations are activated, and the simulation statistics are taken from this point.

Results

Heat and mass transfer

The instantaneous temperature elds at the domain boundaries with interface visualization (zero level set isocontour) at t 3 for each case are presented in Figure 1 (top). At this point in the simulations, the eects of the turbulent velocity eld are easily appreciated in the temperature transport. The interface visualization opacity is reduced to show the bubbles enclosed within the larger liquid structures. Due to the reduced gas volume, these bubbles are characterized by a higher saturation level than the gas far from the interface. The inuence of the gas volume on time required to reach equilibrium is investigated in more detail by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF].

Additionally, a wide range of vapor mass fraction gradients is observed near the interface. These variations depend mainly on the interface's deformation and the surrounding gas's temperature. For example, a liquid structure moving toward a high-temperature gas will exhibit a smaller thermal boundary layer thickness at the front due to the larger temperature gradients. Similar phenomena are observed in Figure 1 (bottom) for the instantaneous vapor mass fraction elds at the domain boundaries at t 3 . Compared to the temperature elds, larger gradients are present in the gas phase: convection eects are more dominant for this scalar because the mass diusivity (D m = 3.82 × 10 -7 m 2 • s -1 ) is smaller than the initial thermal diusivity (D th = 1.70 × 10 -6 m 2 • s -1 ).

Figure 2 shows the volume rendering of the vapor mass fraction at the same time as Figure 1. Here, we can observe the interactions of the vapor mass fraction and the turbulent velocity eld within the domain. There is vapor accumulation near the interface region, and the vapor begins to saturate the entire gas volume until the system reaches equilibrium. Another interesting phenomenon is the vapor wake behind the smaller liquid structures as they are convected through the domain. There are several possible outcomes for the lifetime of these smaller liquid structures, e.g., they may remain convected by the velocity eld until their complete evaporation, or they could merge with another liquid structure. These possible scenarios are considering that the system is not at equilibrium. Furthermore, it is noted that the size and shape of vapor wakes depend on the size of the structure, its average velocity, and the distance from another interface.

Recently, similar observations have been made for a single evaporating droplet in a homogeneous and isotropic turbulent environment [START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF]) and for dispersed evaporating droplets in a homogeneous turbulent shear environment [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF]). In [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF], the authors show variations in the evaporation rate due to small deformations caused by the turbulent velocity eld. In addition, their results also reveal an increase in the evaporation rate due to the interaction of the vapor mass fraction with the turbulent veloc-ity eld. These results highlight the limitation of classical vaporization models based on spherical droplets with a constant evaporation rate.

The temporal evolution of the liquid and interface temperature normalized by the equilibrium temperature for both cases are presented in Figures 3. t * = t τ is the time normalized by the eddy turnover time dened as:

τ = k
where (•) is the volume averaging operator, and is the mean turbulent dissipation rate. As predicted, the temperature of the liquid for case 1 decreases as it approaches the equilibrium state (see Figure 3a). In addition, the liquid temperature follows the temperature at the interface, decreasing despite the higher gas temperature because of the heat jump caused by the latent heat. This result shows the accuracy of the Ghost Fluid Method implemented to account for the energy jump condition and the coupling between the liquid and gas temperature.

At the rst instants of the simulation, a sudden drop in temperature is observed before it starts to stabilize after two eddy turnover times. It can be explained by the maximum value of the evaporation rate during the early moments of the simulation. Subsequently, the vapor generated accumulates at the interface, decreasing the evaporation rate and slowing the temperature decrease. At the end of the simulation, there is a negative slope in the evolution of the liquid temperature, which means that heat is still exchanged with the gas phase, i. e., equilibrium has not yet been reached.

An opposite behavior is observed for the liquid temperature in case 2. Figure 3b shows an increase in the liquid temperature caused by the initial temperature gradient between the liquid and gas. Then, the temperature rise slows down as the heat ux from the gas to the liquid and the evaporation rate decrease.

Consequently, the behavior of both cases is the same as the one predicted by the iterative algorithm described in section 4.1.

The temporal evolution of liquid and interface temperature is strongly linked to the evaporation rate due to the energy jump condition at the interface. Thus, the temporal evolution of the surface-averaged evaporation rate ωΓ for both cases is shown in Figure 4. By comparing each case, a signicant dierence between the initial magnitudes of ωΓ is observed, which is explained by the initial liquid temperatures chosen for each case. Nevertheless, both curves have the same general behavior; a maximum value at the beginning of the simulation followed by a decrease as the system approaches equilibrium.

The temporal evolution of the normalized volume-averaged gas temperature

(T * g = T g
Teq ) is shown in Figure 5. Here, a more pronounced temperature decrease is observed in case 2 compared to case 1. The eect of the energy jump and the Stefan ow on the heat transfer across the interface could explain this behavior.

In case 1, the temperature decrease near the interface caused by latent heat results in an important local decrease in the temperature gradient. In other words, the increased evaporation rate slows down the heat transfer. The decrease in the heat transfer due to the Stefan ow has also been observed in the literature by [START_REF] Jayawickrama | The eect of stefan ow on the drag coecient of spherical particles in a gas ow[END_REF][START_REF] Jayawickrama | The eect of stefan ow on nusselt number and drag coecient of spherical particles in non-isothermal gas ow[END_REF] in vaporizing droplet congurations.

Another way to investigate the inuence of the evaporation rate on the evolution of a system from the initial condition to the equilibrium state is to study the temporal evolution of Y * v shown in Figure 6 for the two cases studied. Here, we can observe a more rapid increase of Y * v in case 1. These results are consistent with our previous analysis of the temporal evolution of the mean liquid and gas temperatures and evaporation rate: higher evaporation rates indicate faster vapor generation at the interface, contributing to the increase of Y * v in case 1. In addition, the increased Stefan ow pushes the vapor away from the interface, contributing to vapor mixing. However, the time evolution of Y * v depends not only on the conditions near the interface (T Γ , ωΓ , etc.), but it also depends on the turbulent velocity eld and the interaction between the liquid and gas structures. This last part is investigated in the following sections by analyzing the evolution of the second moment and the PDFs of Y * v . From Figures 5 and6 gaseous thermal diusion coecient D th g = λg ρgcg induced by the increase of the gas density and heat capacity. Consequently, the evaporation regime inuences the dierence between the thermal and mass characteristic times. For example, these times will tend to approach each other more rapidly in case 1 due to the faster increase of the mean vapor mass fraction. These results show that large evaporation rates induce a relatively high Lewis number variation, limiting the assumption that the Lewis number is constant and equal to one, normally used in turbulent combustion mixing models.

The temporal evolution of the normalized mean pressure P * = P Peq for each case is shown in Figure 8. The pressure variations in both cases can be explained by analyzing the time evolution of temperature and vapor mass fraction. In our simulations, the vapor molecular mass is greater than the air molecular mass, increasing the mean molecular mass of the gas phase, which will tend to increase the system pressure as the liquid evaporates. In contrast, the decrease in gas temperature in both cases will tend to decrease the system pressure. In case 1, a slight increase in initial pressure is observed, resulting from the amount of vapor produced, which compensates for the pressure drop due to the decrease in temperature. However, as less liquid evaporates in case 2, the pressure experiences a signicant decrease. These results agree well with the behavior predicted by the algorithm presented in the section 4.1, summarized in the table 3.

Scalar turbulent mixing analysis

Evolution of the scalar variance

This section studies scalar mixing using the normalized vapor mass fraction second moment. This has been previously done for single-phase [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF]; [START_REF] Juneja | A dns study of turbulent mixing of two passive scalars[END_REF]), dispersed [START_REF] Reveillon | Dns study of spray vaporization and turbulent micro-mixing[END_REF][START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF]; [START_REF] Reveillon | Eects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF]) and dense [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]) two-phase ows. In the rst scenario, several power laws for the scalar variance decay can be found in the literature, depending on whether an imposed mean scalar gradient feeds the scalar. This behavior is explained by the dissipation rate term in the variance transport equation. In the DNS of dispersed two-phase ows, Lagrangian techniques are implemented to track the temporal evolution of droplet diameter, velocity, and position. These droplets create vapor pockets as the turbulent velocity eld convects them. This kind of simulation provides good insight into the eects of a non-homogeneous scalar source on the temporal evolution of the scalar variance. However, some limitations arise, for instance, the interactions between the scalar and the liquid/gas interface are not solved, and an idealized model for liquid vaporization is used.

In this work, the interface and internal ows of liquid structures are solved, similar to [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]. Moreover, our formalism coupled temperature, vapor mass fraction, and uid dynamics through the evaporation rate. Consequently, the vapor mass fraction is not a passive scalar, and the mixing and the uid dynamics of the system depend on the mass and thermal boundary layers simultaneously.

Figure 9 shows the temporal evolution of the RMS (Root Mean Square) of the normalized vapor mass fraction, computed as follows:

Y * v = (Y * v -Y * v ) 2
for the two cases studied. The curves in Figure 9 have similar behavior. An increase in the RMS is observed at the beginning of the simulation. Then, it reaches a maximum value after approximately an eddy turnover time, before its decay towards zero. This behavior diers from the results obtained with the single-phase ow but is similar to those obtained with the two-phase dispersed ow simulations. As discussed by [START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF], the initial increase in scalar RMS is caused by the source terms appearing in its transport equation when an evaporating liquid generates vapor. The maximum value of the RMS is reached when the scalar dissipation rate balances the dispersion generation. After reaching the maximum, turbulent mixing becomes more inuential in the system than vaporization, resulting in a decay of the scalar variance.

Moreover, [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF] observed similar behavior in dense two-phase ow simulations. In this case, the initial increase in dispersion is inuenced by vapor generation at the interface and the diversity of boundary layer thickness formed around the liquid structures. Furthermore, the initial liquid volume fraction on the time evolution of the scalar variance is investigated. It is observed that a decrease in the available gas volume causes an increase in the maximum value of the scalar RMS and a sharper subsequent decay, resulting in a shorter saturation characteristic time.

In a more detailed analysis of the curves presented in Figure 9, we observe the evaporation regime's inuence on the scalar variance's temporal evolution.

For a higher evaporation rate, the maximum value of the variance increases and is reached in a slightly shorter time than in case 2. Moreover, the subsequent decay of the scalar dispersion is also more pronounced for the case with a large evaporation rate. A similar observation was made by [START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF] for dispersed two-phase ows. This is explained by the faster generation of vapor in the dry air environment and the subsequent mixing of this vapor. In this work, the scalar variance depends not only on the magnitude of the evaporation rate but also on the temporal evolution of the surface density, the inuence of the Stefan ow on the system dynamics, and the curvature of the interface, resulting in a more realistic model of the atomization, mixing, and evaporation processes. For example, the Stefan ow pushes the vapor away from the interface, contributing to the decay of the scalar RMS and consequently enhancing the mixing process. In addition, the non-homogeneity of the evaporation rate considered in the present work dramatically aects the temporal and spatial evolution of the scalar dissipation rate and, consequently, the scalar variance.

The following sections perform a more detailed analysis of scalar generation and mixing by computing scalar PDFs and joint PDFs.

Evolution of the scalar PDFs

The asymptotic behavior of scalar mixing has been extensively studied and is well known in the literature. For a conned two-phase environment with a vaporizing liquid surrounded by dry air, the scalar PDF starts as a shape similar to a Dirac function where Y * v = 0 if evaporation has not started yet. Then, vaporization occurs and when the equilibrium state is nally reached, the scalar PDF has a similar shape at Y * v = 1, where vapor generation is no longer possible due to the absence of a vapor gradient.

In this work, we are interested in studying the temporal evolution of scalar PDF shapes. For this purpose, the PDF of the normalized vapor mass fraction for four successive instants for case 1 is plotted in Figure 10. The general tendency described above is easily seen. An interesting behavior that is observed for the curves presented in Figure 10a is the presence of values of Y * v > 1. Hence, the vapor mass fraction is not upper-bounded to 1 when normalized by its saturation state value, even though it is its nal value. This behavior has at least two explanations; rst, from the table 3, we can observe that the temperature at equilibrium is lower than the liquid temperature. Consequently, the equilibrium interface vapor mass fraction is smaller than the initial vapor mass fraction.

Secondly, there is a wide range of liquid structures with dierent characteristic lengths, resulting in dierent temporal evolution of temperature proles in each liquid structure. This means that some smaller liquid structures do not have the same behavior as the overall system. For example, a small convected droplet in contact with high-temperature gas will evaporate and increase in temperature simultaneously, even though the system's average temperature tends to decrease. This is also a consequence of solving the temperature and internal ow of the liquid structure. Furthermore, these results demonstrate the strong coupling between the energy and species conservation equations and the variety of phenomena captured by our compressible formalism.

In the analysis of passive scalar mixing in single-phase ows, values of the normalized vapor mass fraction greater than unity are generally not observed.

In the literature, the total amount of vapor is usually imposed at the initial condition of the simulation. In the simulation performed by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF], the vapor mass fraction at the interface is constant and imposed during simulation initialization. Consequently, the temporal evolution of the temperature and its relation with the vapor mass fraction is not solved. As a result, the equilibrium vapor mass fraction equals the imposed vapor mass fraction at the interface. The following PDF corresponds to t * 2 = 2.88 and an RMS of Y * v = 0.22, which means that the decay has already started; thus, the contribution of the scalar dissipation rate is more signicant than the contribution of the vaporization source term in the variance transport equation. In this state, the eects of the turbulent velocity eld on the scalar mixture are more critical. As a result, the vapor is transported into the gas eld away from the interface, resulting in a rightward shift and a decrease in the PDF peak. In addition, the tail of the PDF has decreased, indicating a reduction in the vapor mass fraction at the interface.

P (Y * v ) Y v * t * 1 -→ t * 4 (a) 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 Y * v 0.0 0.5 1.0 1.5 2.0 2.5 P (Y * v ) (b)
Subsequently, at the third PDF, which corresponds to a time of t * 3 = 6.21, regrowth of the peak and a shrinking of the PDF are observed, showing a homogenization of the vapor mass fraction in the gas phase away from the interface, which is also characterized by a decay of the RMS (Y * v = 0.14). In addition, the probability of values smaller than Y * v = 0.25 is close to zero, which means that the vapor is spread over the entire gas volume. A tail is also present on the right side of the PDF, corresponding to the larger values close to the interface.

Finally, the shrinking continues in the fourth PDF: Y * v = 0.077 at t * 4 = 12.0, which is consistent with the asymptotic behavior described above. The PDF peak is larger, and its shape has fewer irregularities, indicating a higher level of homogenization of the vapor mass fraction. An interesting behavior is observed at the end of the right tails of the PDFs.

A zoom of the highest values of the Y * v is presented in Figure 10b. This way, the temporal evolution of the PDF tails is easily seen. Here, we see that the PDFs tails move to the left, towards unity. In addition, there is an abrupt behavior change in the last values of the PDFs. This could be interpreted as the result of the temporal evolution of the vapor mass fraction gradient near the interface. It is worth mentioning that these results describe a general behavior of the system and that the mass boundary layer thickness varies depending on the local thermodynamic conditions in each structure and each cell close to the interface.

The PDFs of the vapor mass fraction for the case 2 at the same level of mixing are presented in Figure 11. In this case, values of Y * v > 1 can also be found but with a smaller probability than in case 1. Here, the equilibrium temperature is higher than the initial liquid temperature (see table 3). Consequently, a probable reason for having values of Y * v > 1 is that the smallest liquid structures are heated at temperatures higher than the equilibrium temperature, which generates a high vapor mass fraction at the interface. This occurs especially in the rst moments of the simulation. The probabilities for the values of Y * v > 1 decrease with time until they become nonexistent in the last PDF.

Comparing the two studied cases, the inuence of the evaporation rate in the PDF shapes can be observed. In case 2, the physical time needed to arrive at the same value of Y * v (t * 1 = 1.79) is almost 13 times higher than in case 1 due to the lower evaporation rate. As a consequence, the turbulent velocity eld has more time to advect the vapor mass fraction eld, resulting in a reduction of the PDF peak. Similar behavior is observed in the other PDFs. The second PDF correspond to a time of t * 2 = 4.73 and a Y * v = 0.137. At this time, the peak of the PDF continues to decrease. In the third and fourth PDFs (t * 3 = 8.99 and t * 4 = 15.26), the width has been reduced considerably. This indicates that a more ecient turbulent mixing is achieved compared to case 1, but their behavior over time remains similar.

Self-similarity

Figure 12 shows the normalized PDF of

Y * v (P (Y * v )Y * v with respect to Y * v -Y * v Y * v
). Similar to the results obtained by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF], a signicant asymmetry is observed in all the PDFs. The peaks are on the negative side of the x-axis, indicating that the smaller values far from the interface are dominant. In Figure 12a, the shape of the PDF is aected by the initial conditions in the early stages of the simulation. However, the PDF adopts similar shapes once the scalar mixing is more developed.

Nonetheless, around the peak, we can notice a dierence in the shapes of the PDFs, especially in the low evaporation case. Compared with the results of the single-phase simulations, it is evident that the interface acting as a vapor source strongly inuences the PDF shape. For example, the tail present on the right of all PDFs is the primary indicator of these dierences.

The shape observed in Figure 12 is similar to the results obtained by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF]. However, there is a noticeable dierence on the right-hand side of the PDFs. In the work of [START_REF] Duret | Dns analysis of turbulent mixing in two-phase ows[END_REF], a secondary peak is formed at the maximum value of the scalar due to the interface boundary condition of the vapor mass fraction. In this work, the tail shape of the PDFs varies as a function of the vapor mass fraction at the interface, which is non-homogeneous and dependent of the interface temperature. Consequently, the vapor mass fraction at the interface possesses its own PDF, which will be investigated in the next section. In both cases, the shapes of the PDFs are similar to a log-normal distribution, conrming that the PDF shape is not adequately represented by a beta distribution when the interface acts as a scalar source.

Surface conditional PDFs

This section analyzes the temporal evolution of the evaporation rate and the vapor mass fraction at the interface. For this purpose, the PDFs for the same instants presented in the previous section are used.

Figure 13 shows the PDFs of the surface-averaged normalized vapor mass fraction Y * vs = Yvs Yv,eq at t * 1 , t * 2 , t * 3 and t * 4 for both cases. For the interface vari--4 ables, the asymptotic behavior of the PDF is similar to the gas phase variables.

-2 0 2 4 6 (Y * v -Y * v ) Y * v 0.0 0.2 0.4 0.6 0.8 1.0 1.2 P (Y * v ) Y * v (a) -4 -2 0 2 4 6 (Y * v -Y * v ) Y * v 0.0 0.2 0.4 0.6 0.8 1.0 1.2 P (Y * v ) Y * v (b)
At rst, the PDF is spread around the value corresponding to the initial conditions. Then, the PDF tends to a Dirac function when the system is close to the equilibrium state.

Large dispersion of the vapor mass fraction is observed in the rst PDFs due to a wide range of interface temperatures. Moreover, the mean vapor mass fraction moves toward its value at equilibrium (Y * vs = 1) in the successive instants. In contrast, the dispersion decreases due to the homogenization of the interface temperature, resulting in the increase of the PDF peak with time.

Comparing Figures 13a and13b, while the general trend is similar as the system approaches equilibrium, a dierence in PDF behavior is observed early in the simulation. In case 2 there is an initial increase in the mean value of Y * vs followed by a successive decrease. On the contrary, in case 1, the mean value of Y * vs decreases from the beginning of the simulation to the end. This initial increase of Y * vs in the second case could be attributed to the pressure drop of the system, which has the tendency to increase the vapor molar fraction (see equation 11).

The PDFs of the evaporation rate for the four studied times are presented in Figure 14 for both cases. In this context, the PDF is expected to evolve to a Dirac function shape at ω = 0 when the system reaches equilibrium. An important dispersion of the evaporation rate is observed at the earlier moments due to the high temperature and vapor mass fraction gradients. In addition, a signicant asymmetry is observed in all PDFs with a peak on the left side and a long tail to the right, showing the low probability of the highest values of the evaporation rate. For successive instants, the decrease of the mean evaporation rate and the tail width is observed, causing an increase of the PDF peak. The main dierence between both cases is that the magnitude of ω is more signicant in case 1.

Compared to the interface vapor mass fraction, the evaporation rate also depends on the interface temperature, but it also depends on the vapor concentration around the interface, and the interface geometry, such as the curvature

P ( ω) ωΓ t * 4 ←-t * 1 (a) 0.5 1.0 1.5 2.0 2.5 3.0 ω kg m 2 s 0 1 2 3 4 P ( ω) ωΓ t * 4 ←-t * 1 (b)
Figure 14: Instantaneous PDFs of ω at t * 1 (solid t * 2 (dashed line), t * 3 (dashed-dotted line) and t * 4 (dotted line) for case 1 (a) and case 2 (b).

and characteristic lengths of the liquid and gas structure.

Inuence of mean curvature on the evaporation rate

To better understand the simulation state at the studied times, the interface visualization colored by the evaporation rate magnitude for the same instants as the PDFs are presented in Figure 15. As for Figure 1, the eect of the turbulent velocity eld on the interface deformations is observed. Moreover, the overall decrease of the evaporation rate as time passes is shown. Furthermore, the larger values of the evaporation rate are found in the negative values of the mean curvatures. Hence, the probability of having a large evaporation rate is higher in small liquid structures than in the concave parts of the interface, usually present in gas structures. This could be explained by the formation of local saturation zones in concave interfaces. Indeed, the vapor generation is limited by the magnitude of the scalar gradient at the interface: in small gas structure the scalar generation is rapidly stopped when the gas structure is completely at saturation, which happens quickly if the gas structure is closed or quasi-closed. Consequently, the scalar gradient tends to zero, reducing the magnitude of the vaporization rate. Even if the gas structure is large and not saturated, local saturation can occurs in strongly curved concave interface region : the vapor is "trapped" and scalar gradients are greatly reduced as well as the vaporization rate. In the second joint PDF (top right of Figure 17), there is an increase in the probability of ω in the negative values of the mean curvature (convex interface) whereas the probability is reduced in the positive curvature (concave interfaces), which indicates an increased saturation state close to concave interface as time evolve. Then, in the third joint PDF (bottom left of Figure 17), the decrease of ω is more noticeable: the simulation is closer to the equilibrium state. Moreover, the values ω are less dispersed because of the homogenization of the vapor mass fraction and temperature elds. Finally, in the last joint PDF (bottom right of Figure 17), the probability of both convex and concave interfaces start to become comparable, and the tendency of ω to tend to zero is more evident.

In Figure 18, the temporal evolution of the joint PDF of the evaporation rate for case 2 is presented. Similar behavior to that shown in the Figure 17 is observed. However, the eects of the evaporation regime are easily seen.

For example, the magnitude of the evaporation rate is considerably smaller, as noted in the previous section (see Figure 4). In the earlier times, the shape of the joint PDFs for both cases is similar despite the dierence in the magnitude of the evaporation rate. Once the turbulent mixing begins homogenizing the scalars, the dierences in the joint PDFs become more apparent. For the same value of the mean vapor mass fraction, we can observe the attened shape of the joint PDFs, which indicates a higher saturation level for the same value of the mean vapor mass fraction.

Conclusion

In this work, a fully compressible DNS formalism is applied to analyze the evaporation process and scalar turbulent mixing in an atomization regime. The method is able to capture the coupling between temperature, vapor mass fraction, and momentum conservation equations; resulting in a non-passive scalar mixing and a realistic representation of the vaporization process. Hence, our work is one of the few available in the literature that considers all these phenomena simultaneously.

The scalar turbulent mixing and the evaporation process are analyzed. To this end, the temporal evolution of the vapor mass fraction RMS for both cases is studied. An initial increase in the vapor mass fraction RMS results from the vapor production at the interface in a primarily dry air environment. Afterward, the scalar RMS decays as the equilibrium state is reached. Moreover, it was observed that the evaporation regime aects this behavior by increasing the maximum value when large evaporation rates are considered, and it is accompanied by a faster decay.

Furthermore, the temporal evolution of the PDFs of the vapor mass fraction is studied. It is shown that the PDFs shape strongly depends on the mixing stage, but remains close to a log normal distribution. An interesting behavior is observed in the right tail of the PDFs. Here, the shape of the PDFs tails depends on the vapor accumulation and generation near the interface. In addition, it is shown that the evaporation rate aects the temporal evolution of scalar PDFs;

resulting in an increased physical time to reach the equilibrium value of the vapor mass fraction when the evaporation rate decrease.

Finally, a statistical analysis of the surface averaged vapor mass fraction and evaporation rate is performed. The temporal evolution of their PDFs is analyzed. Both variables have shown a more signicant dispersion in the rst instants of the simulation. In the second case, the important pressure variation aects the vapor mass fraction generated at the interface: initially the PDFs means increase due to the pressure drop that enhances the scalar source, then a decrease is observed due to the lower scalar generation induced by lower temperatures in the system.

In addition, the behavior of the evaporation rate is studied from the curvature point of view, observing that the probability of the highest values of the evaporation rate is found in the convex interface curvatures. These results illustrate the impact of the assumptions typically made in the literature, where constant variables such as vapor mass fraction and evaporation rate are assumed at the interface. These assumptions appears to be not adequate in turbulent two-phase environments with high evaporation rates. Further work will be focus on improving vaporization modelling for dense two-phase ows, based on our DNS data.

Figure 1 :

 1 Figure 1: Instantaneous temperature (top) and vapor mass fraction (bottom) elds at the domain boundaries with interface visualization (zero level set isocontour) at t 3 for case 1 (left) and case 2 (right).

Figure 2 :

 2 Figure 2: Volume rendering of the vapor mass fraction eld with interface visualization (Level Set 0 isocontour) at t 3 for the two cases studied.

Figure 3 :

 3 Figure 3: Normalized liquid (solid line) and interface (dashed line) temperatures temporal evolution for case 1 (a) and case 2 (b).

Figure 4 :

 4 Figure4: Evaporation rate temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 5 :

 5 Figure5: Normalized gas temperature temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 6 :

 6 Figure 6: Normalized vapor mass fraction temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 7 :

 7 Figure7: Gas Lewis number temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 8 :

 8 Figure8: Normalized pressure temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 9 :

 9 Figure 9: Temporal evolution of Y * v for case 1 (dash line) and case 2 (solid line).

Figure 10 :

 10 Figure 10: (a) Instantaneous PDFs ofY * v at t * 1 (solid line), t * 2 (dashed line), t * 3 (dashed-dotted line) and t * 4 (dotted line), and (b) Zoom of the higher values of Y * v at the same times for case 1.

  Figure 11: (a) Instantaneous PDFs ofY * v at t * 1 (solid line), t * 2 (dashed line), t * 3 (dashed-dotted line) and t * 4 (dotted line), and (b) Zoom of the higher values of Y * v at the same times for case 2.

Figure 12 :Figure 13 :

 1213 Figure 12: Instantaneous normalized PDFs of Y * v at t * 1 (solid line), t * 2 (dashed line), t * 3 (dasheddotted line) and t * 4 (dotted line) for case 1 (a) and case 2 (b).

Figure 15

 15 Figure15aims to understand the inuence of the curvature on the evaporation rate. In this work, the convention represented in Figure16is used: red represents negative values of the mean curvature corresponding to a convex interface, and blue represents positive values of the mean curvature corresponding to a concave interface. In Figures15, adecrease in the evaporation rate is observed when the interface is concave. Also, we can see some gas structures enclosed in the liquid and the low evaporation rate in these areas.

Figure 17

 17 Figure17shows the joint PDF evolution of the evaporation rate ω as a function of the normalized mean curvature H * . Here, the mean curvature is normalized with the radius of the droplets used in the initialization of the HIT conguration for φ = 0.1 (see[START_REF] Luret | Modeling collision outcome in moderately dense sprays[END_REF];[START_REF] Canu | Curvature-based interface resolution quality (irq) indicator to assess simulation accuracy[END_REF] for more details). The joint PDFs are computed at the time instants corresponding to Y * v = 0.18, 0.32, 0.54 and 0.72, and are presented in Figure17from left to right and from top to bottom. The rst joint PDF (top left of Figure17) represents the initial stages of the mixing process, and a wide dispersion of ω over a large range of mean curvatures is observed. Additionally, ω dispersion is related to the variations of local vapor concentration close to the interface, the interface deformations, and the diversity of liquid and gas structure characteristic lengths.Similar observations were made by[START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] for dispersed twophase ow simulations in homogeneous shear turbulence. In their simulations, the liquid/gas interface and internal ow of the liquid droplets are resolved, and the evaporation rate depends on the vapor mass fraction and the temperature elds.[START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] observed dispersion of ω due to the interface deforma-

Figure 15 :

 15 Figure 15: Instantaneous interface visualization (zero level set iso-contour) colored by the magnitude of the evaporation rate at t * 1 , t * 2 , t * 3 and t * 4 (from left to right, and top to bottom) for case 1.

Figure 16 :

 16 Figure 16: Graphical representation of the curvature sign of a deformed liquid structure.

Figure 18 :

 18 Figure 18: Instantaneous joint PDFs for ω in function of the curvature for case 2 at t * 1 , t * 2 , t * 3

Table 1 :

 1 Physical properties for the air, vapor, and n-decane at the initial state of the simulation.

	Fluid	Sc	P r	Le
	Air	1.96 0.44 4.44
	N-decane	-	8.90

-

Table 2 :

 2 Dimensionless numbers at the initial state of the simulations.

  eq

						m l,eq	Peq
						m 0 l	P 0
	1	540	573	522	0.44	0.76	1.03
	2	460	573	464	0.19	0.93	0.83

Table 3 :

 3 Dierent variables at the equilibrium state for the two cases studied.

http://docs.mercurve.rdb.is/

Acknowledgements

This work was granted access to the HPC resources of IDRIS, TGCC, and CINES under the allocation (A0112B10101) made by GENCI (Grand Equipement National de Calcul Intensif). The CRIANN (Centre Regional Informatique et d'Applications Numeriques de Normandie) is also gratefully acknowledged (project No. 2017004) for its CPU resources. Also, the authors thank the MESCYT (Ministerio de Educación Superior, Ciencia y Tecnología, Dom. Rep.) and professor Bruno Renou who have supported and managed this work in the framework of the Caliope program.


Helen Mair Rawsthorne
email: helen.rawsthorne@ign.fr

Nathalie Abadie

Eric Kergosien

Cécile Duchêne

Éric Saux

ATONTE: towards a new methodology for seed ontology development from texts and experts

Keywords: application ontology, low-level ontology development, ontology development methodology, seed ontology, technical domain

ATONTE (ATlantis methodology for ONtology development from Texts and Experts) is a methodology for the manual development of low-level seed ontologies. The modelling process is based on a combination of knowledge from non-fiction text corpora such as manuals, information guides or sets of instructions, and the knowledge of domain experts. This article presents the five key steps of the ATONTE process. Seed ontologies created with ATONTE can be used to develop and populate knowledge graphs for use in specific applications within given technical domains.

Introduction

ATONTE (ATlantis methodology for ONtology development from Texts and Experts) is a methodology for the manual development of low-level seed ontologies based on non-fiction text corpora such as manuals, information guides or sets of instructions, as well as the knowledge of domain experts. Seed ontologies created with ATONTE can be used as to develop and populate knowledge graphs for use in specific applications within given technical domains.

Current methodologies for ontology construction from text primarily offer automatic and semi-automatic approaches [START_REF] Al-Arfaj | Ontology Construction from Text: Challenges and Trends[END_REF]. They employ techniques in machine learning or natural language processing to perform statistical or linguistic analyses on the text. However, automatised approaches pose several problems. Firstly, they assume that the corpus contains precisely the necessary knowledge to satisfy the application without surrounding noise. Secondly, the type of knowledge found in technical corpora can be complex, which is difficult to model automatically. Thirdly, using an automatic approach leaves the door open to unexpected errors that could put the integrity of the content at risk, and errors in the modelling of texts that include information relating to safety and security could put human lives at risk. On top of these three issues, automatised approaches do not allow for the easy integration of domain experts' knowledge. Such knowledge can be obtained by conducting targeted interviews with domain experts, and the knowledge they provide can supplement and clarify knowledge extracted from textual sources. A stable and dependable manually-created seed ontology could nevertheless still be enriched semi-automatically or automatically.

We thus present ATONTE as a reverse methodology that integrates domain experts' knowledge and manual modelling techniques during the early stages, and that depends on automation only towards the end in order to complete, instantiate and verify the ontology. To develop ATONTE, we took inspiration from SAMOD (Simplified Agile Methodology for Ontology Development) [START_REF] Peroni | A Simplified Agile Methodology for Ontology Development[END_REF], MOMo (Modular Ontology Modeling) [START_REF] Shimizu | Modular Ontology Modeling[END_REF] and NeOn (Networked Ontologies) [START_REF] Suárez-Figueroa | The NeOn Methodology for Ontology Engineering[END_REF].

ATONTE is still under development and is being refined empirically. We are currently testing the methodology via the creation of the ATLANTIS (coAsTaL mAritime NavigaTion InstructionS) ontology [START_REF] Rawsthorne | ATLANTIS : Une ontologie pour représenter les Instructions nautiques[END_REF]. ATLANTIS models the knowledge contained in the Instructions nautiques, a series of books published by the Shom (the French Naval Hydrographic and Oceanographic Service) that contain information on navigating in coastal waters.

Proposed methodology

In this section we present the five key steps of ATONTE, which are illustrated in Figure 1. Some elements from SAMOD are reused in steps 3 to 5, concepts from MOMo are incorporated into step 4, and elements from NeOn are integrated into step 5.

Choosing ATONTE

Each of the three following requirements should be satisfied for it to be a good choice of ontology development methodology for a given task.

1. The aim of the ontology is to model some or all of the knowledge contained in a nonfiction text corpus, and combine it with the knowledge of domain experts, with a given application in mind. 2. The end users of the ontology application are known and can be consulted, along with domain experts (can be the same people), on an ad-hoc basis during the development process. 3. No other ontologies that could be used for the final application, or part of it, already exist and have been published on the Web.

Groundwork

Define the ontology application. Interview end users of the application to find out what knowledge from the corpus needs to be modelled in the ontology, what knowledge needs to be added (if any) and how it all needs to be structured in order to be useful to them. Unless the knowledge spans only a very small domain, divide it into coherent subdomains. This will facilitate the documentation and modelling phases.

Documentation

For each subdomain, write a motivating scenario, a list of competency questions and a glossary. Use the results of the interviews to guide the writing of this documentation.

In the motivating scenario, include:

• a name for the subdomain, • an explanation of the motivation behind modelling the subdomain for the application in question and with the end users' needs in mind, • a list of the principal characteristics of the concepts within the subdomain found in the corpus, and • a set of excerpts from the corpus that demonstrate all the ways in which the subdomain knowledge is represented in the text.

Next, draw up a list of natural language competency questions that encapsulate all the uses for the knowledge mentioned by the end users during the interviews. The answers to the questions should figure in the excerpts in the motivating scenario. Finally, compile a glossary defining all the technical terms used in the subdomain. Enrich and validate the documentation with the help of domain experts.

Modelling and testing

For each subdomain, take the excerpts from the motivating scenario one at a time and semiformalise the knowledge contained within them that is relevant to the subdomain. Do this by breaking down the relevant knowledge in the excerpt into finer-grained chunks, favouring a subject-predicate-object structure where possible. Once all the excerpts have been semiformalised in this way, group together the subjects/objects and the predicates that serve the same purpose to create the first set of classes and properties for the subdomain: this is the first iteration of the subdomain model. Rewrite the semi-formal expressions formally as triples. This is the first iteration of the set of triples. Figure 2 shows an example of this process of semi-formalisation and then formalisation of two excerpts of text. Now search for another set of subdomain-relevant excerpts in the corpus and try to formalise their content by creating triples using the newly-created classes and properties. If this task cannot be performed satisfactorily, modify the model accordingly whilst ensuring that it still fits all the other excerpts. Continue in an iterative fashion with unseen sets of excerpts until the model has stabilised and you have a solid set of triples specific to the subdomain. Submit each subdomain model to a series of three tests: a model test, a data test and a query test, in that order. For the model test, use a reasoner to verify the consistency of the model and then manually check that the model corresponds to the motivating scenario written for it. For the data test, verify the validity of the model by populating it with the triples created from the corpus. For the query test, translate the natural language competency questions into SPARQL queries and run them on the manually-created set of triples for that subdomain to check that the results match the answers specified in the documentation. Move on to the next test only once the previous test has been passed. If the subdomain model fails a test, return to the modelling phase to fix the issue before testing again.

Merging and refactoring

Now the subdomain models can be merged to create the final full model. Start with the largest subdomain model, merge the second largest into it and then perform the series of tests on this intermediate model. Repeat this process of merging the next-largest subdomain model into the intermediate model and then testing until all subdomain models have been integrated and the Subdomain Geographical and Meteorological Features Excerpt 1 "The current in the bay flows eastward, but the wind is northerly. " Excerpt 2 "The port features a lighthouse. " Semi-formalisation current 1 -is a -current / bay 1 -is a -bay / current 1 -is in -bay 1 / current 1 -flowseastward / wind 1 -is a -wind / bay 1 -has -wind 1 / wind 1 -is -northerly / port 1 -is a -port / lighthouse 1 -is a -lighthouse / port 1 -features -lighthouse 1 Formalisation id:current_1 :hasType :Current / id:bay_1 :hasType :Bay / id:current_1 :isLocatedIn id:bay_1 / id:current_1 :hasDirection :east / id:wind_1 :hasType :Wind / id:bay_1 :contains id:wind_1 / id:wind_1 :hasDirection :north / id:port_1 :hasType :Port / id:lighthouse_1 :hasType :Lighthosue / id:port_1 :contains id:lighthouse_1 final full model has been created and tested. The merging process is done by exporting the .owl files of the models and manually combining them, removing all duplicate classes, properties and individuals in the process.

The refactoring process involves reusing existing knowledge in semantic resources, annotating the model and enriching it using the capabilities of the OWL language, for example to create inferences. The elements of the refactoring process are given in SAMOD [6, p. 11]. A detailed description of how to reuse existing semantic knowledge resources is given in NeOn [START_REF]Modularizing Ontologies[END_REF]. Figure 3 shows some examples of this refactoring process. After the refactoring process has been carried out, the model should undergo a final testing cycle.

Conclusion and perspectives

The ontology development methodology ATONTE allows the creation of application seed ontologies from non-fiction text corpora and the knowledge of domain experts. One characterising feature of this methodology is the early creation of an exemplar dataset, upon which the modelling of the ontology is based. Another is the use and combination of knowledge from a text corpus with the knowledge of domain experts.

ATONTE is currently being used to create the ATLANTIS ontology, which models the contents of the Shom's Instructions nautiques. ATLANTIS will offer new possibilities for the access and the verification of the knowledge contained within the series of Instructions nautiques. For example, an online platform giving access to ATLANTIS could allow future navigators to search directly for the specific information they are looking for, by category or by geographic area, instead of having to read the full texts. Domain experts' knowledge has been vital in the construction of the ontology for this particular application, in order to understand how best to model specific parts of the knowledge within the corpus. ATLANTIS currently contains 110 classes, 90 object properties, 90 data properties and 2190 axioms in total. A full evaluation of ATLANTIS is underway and will be the subject of a future publication.

The short-term perspectives of this project include improving ATONTE by adding a final ontology evaluation phase. It will entail automatically extracting the relevant information from the corpus being modelled in order to populate a knowledge graph. The seed ontology will thus be enriched with new concepts and properties from other parts of the corpus that were not explored manually. Only once the ontological model is instantiated will it be possible to carry out a rigorous evaluation of the ontology produced with ATONTE.

Figure 1 :

 1 Figure 1: A flowchart showing the steps of the ATONTE process to produce a seed ontology.

Figure 2 :

 2 Figure 2: An example of the semi-formalisation and formalisation of two sample text excerpts belonging to a subdomain called "Geographical and Meteorological Features".

Figure 3 :

 3 Figure 3: An example of refactoring part of a sample model: aligning classes and properties with existing semantic resources and exploiting the capabilities of the OWL language.

Acknowledgments

This work is co-financed by the Shom and the IGN and is being carried out at the LASTIG, a research unit at Université Gustave Eiffel.


Jennifer Mackie
email: jennifer.mackie@mines-albi.fr

Olivier De Almeida
email: olivier.dealmeida@mines-albi.frc

Fabrice Schmidt
email: fabrice.schmidt@mines-albi.frd

Karine Labastie
email: karine.labastie@irt-saintexupery.com

Modeling of IR Lamps with Coated Reflector Used in the Slurry Powder Impregnation Process of Composite Tapes

This study proposes a modelling strategy to simulate the heating stage during the production of thermoplastic composite tapes. Impregnation using a slurry powder technique with carbon fibres and PEKK (PolyEther-Ketone-Ketone) requires a heating step, achieved with an infrared (IR) oven, to evaporate the water and melt the polymer powder. These phenomenon are highly temperature dependant justifying the need to characterise heat transfer within the infrared oven.

. Temperature distribution results were transferred into the commercial software COMSOL Multiphysics® to estimate model parameters: filament temperature and an emissivity distribution function.

INTRODUCTION

In recent years there has been a growing interest for the use of composite materials in the aeronautical industry. Thermoplastic composites are of particular interest as they can be welded or reshaped and recycling is possible. Processes like automated tape laying use thermoplastic composite tape to create larger structure [START_REF] Stokes-Griffin | [END_REF] [3]. Different processes can be used to create these tapes such as fusion impregnation, fiber commingling, film stacking or powder impregnation [4]. In the case of wet powder impregnation (slurry powder) the wet impregnated tape undergoes a series of heating steps to evaporate the water [5] and melt the polymer powder [6]. These phenomenon are highly temperature dependent and to better understand the evolution of the product during these steps it is first necessary to understand the heating conditions applied.

In the slurry powder impregnation process under study, the thermoplastic tape is heated using radiative technology. An oven composed of multiple infrared halogen lamps, reflective panels, glass plates, cooling air injection and air extraction system is used. Here the lamps have an integrated reflector in the form of a fine coating on the underside of the glass tube contrary to clear lamps with Lambertian emission that are more commonly used [7],. This reflector modifies the spatial distribution of the lamps emission. Methods such as ray tracing [8] can be used to model this type of emission but is highly time consuming, especially for a high number of emitters. This motivated the need to find a simpler model capable of predicting the lamps emission with reasonable computation time.

This study then proposes a modeling approach based on the radiosity method for the case of lamps with integrated reflector. The strategy used consisted in applying a weight factor to the emitted power around the circumference of the lamp so as to represent the spatial distribution of the emitted radiation. The values of the different weight factors were identified from the comparison of numerical data with experimental data recorded with a specific experimental set up by using an inverse analysis method.

First, the infrared lamps and experimental procedure are presented. Then, a focus on experimental results is carried out before developing on the numerical model and identification of the model parameters.

EXPERIMENTAL PROCEDURE Infrared Lamps

'Clear' halogen infrared lamps are composed of a tungsten filament, a tubular quartz case and a small quantity of halogen gas. Because of the helical shape of the tungsten filament, the radiative emission of this type of lamp is Lambertian and mostly depends on filament temperature. This temperature in turn is determined by the supply voltage as the lamp is heated using Joule heating. In this study halogen IR lamps (Dr. Fischer-2000W, 400V) with an integrated reflector are studied (Fig. 1). The integrated reflector is a fine coating on the underside of the quartz tube. This coating effects the overall emission of the lamp by redirecting part of the filaments radiation towards the front and away from the back. As the filament temperature varies with supply voltage so does the spectral emission of the lamp. For nominal filament temperature the wavelength of maximal emission is λmax=1.2 µm with 95% of emission in the range [0.6;6] µm. By determining the relationship between supply voltage and filament temperature the evolution of this range will also be determined.

Experimental setup

The experimental setup used to evaluate the infrared lamp's filament temperature and spatial distribution function is presented in Fig. 2. The infrared heating was achieved by controlling the supply voltage U (V) applied to the tungsten filament. The lamp thus heated a flat panel on its front side, and the temperature on the back side of the sheet was measured using an IR camera (FLIR SC325, [7.5-13] µm). The lamp was directed with the 'clear' side facing upwards and the reflector facing downwards as can be seen in Fig. 1, and was positioned centrally to the panel.

FIGURE 2. Experimental set-up

Material

The material used in the experimental set up was a flat sheet of ABS (Acrylonitrile Butadiene Styrene), an amorphous thermoplastic polymer. The main thermal and optical properties of ABS can be found in Table 1. The panel used in this study was 680 x 670 x 1.5 mm.

EXPERIMENTAL RESULTS

The temperature field on the back face of the ABS sheet was measured for supply voltages ranging from 20% to 100% of nominal voltage (400V). For each experiment the sheet was heated until a temperature of 50°C was reached, the lamp was then turned off and the cooling of the plate was recorded. This temperature was chosen to assure the panel maintains a much lower temperature than the glass transition temperature of the material For certain measurements, the power of the lamp did not allow the ABS to reach the maximum temperature desired, in this case the lamp was shut off when the temperature variation was negligible however it is important to note that a steady state was not reached.

Figure 3 shows the temperature field on the back face of the ABS sheet in the last instants of heating. The point of maximum temperature is shifted vertically towards the top of the panel as an effect of the reflector. This point can be seen on the example given in Fig 3(a), it serves as the intersection point for perpendicular and parallel profiles that will be used for the optimization step. The evolution of the temperature of this point versus time can be observed Fig.

NUMERICAL SIMULATION Model

To model the Lambertian emission of a 'clear' lamp only the filament needs to be considered [START_REF] Nakouzi | Modélisation du procédé de cuisson de composites infusés par chauffe infrarouge[END_REF]. The spiraled filament is approximated by an equivalent cylinder, to which a temperature T (K) and integrated emissivity εf (T), is applied in order to define the lamp's radiosity J (W.m -2 ) according to Eq.1 where σ=5.670367e-8 W.m -2 K -4 is the Stephan-Boltzmann constant.

) ( T T J f    (1)
For a lamp with a coated reflector, this model has to be extended. A model capable of taking into account spatial effects in the form of a distributive function on the filaments emissivity was thus proposed. As described with Eq. 2, the filament circumference was divided into n zones (6 zones in Fig. 4(b)), emissivity was weighted by a coefficient ki. The sum of these coefficients should be equal to the number of zones to conserve transmitted power by the filament (Eq. 3).

) (T k

f i i    (2)   6 1 6 1 i i k (3)
Numerical simulation was carried out using commercial software COMSOL Multiphysics®. Boundary conditions were assumed only radiative and convective on the front and back faces of the sheet (Eq.4), thermal insulation was applied to the four edge faces.

) ( ) ( ) ( 4 4 T T T T h T k n air air       (4)
As a first approach, convection in this study was taken to be constant and equal to an average value h=4.5 W.m -2 .K -1 calculated using the data from the cooling phase of the tests.

In this model 400 hexahedral elements were used to mesh the panel. The Biot number was estimated and was found to be inferior to 0.1, surface heat transfer is there for dominant over conduction within the material, simulations were carried out with only one mesh in the thickness of the ABS sheet to reduce computation time. Filament temperature is assumed constant and homogeneous during the simulation. The filament surface was meshed with 360 quadrilateral elements.

Optimization

Parameter optimization was carried out using the SNOPT solver [START_REF] Gill | User's Guide for SNOPT Version 7: Software for Large Scale Nonlinear Programming[END_REF] implemented in COMSOL Multiphysics®. The work presented in this section was carried out for a supply voltage U=160V with a Least-Squares objective function (Eq. 5). Initial and optimized values of the control variables can be found in Table 2. This optimization was carried out considering the coefficients ki positive and inferior or equal to unity, however the condition given by Eq. 3 was not implemented as a first approach. The surface temperature distribution at the end of heating time th is shown in Fig, 5 (a), the overall profile is very similar to that of the experimental results. Figure 5 (b) shows a comparison between the experimental data and the simulated temperature for the perpendicular profile. The agreement is fair between the simulated profile and the experimental data with the exception of the maximal peak. The global error on the temperature taken over the whole face was found to be 2.8%.

  2 exp ) ( min num T T (5)
The reflector redirects a percentage of the filament's emission towards the front of the lamp but the coating is not opaque, as a result the remaining emission is transmitted through the back face of the lamp. The optimized values found in this study do not physically represent this situation as k4 and k5 are null. In the model this is compensated by the remaining coefficients.

The lamp was turned by 90° so that the clear side was fully facing the ABS sheet. Using the optimized coefficients from Table 2, an error of 26.4% was found between experimental and numerical data. In further work optimization of the parameters will be carried out using Matlab®, data from multiple lamp positions will be used for the optimization and the constraint Eq. 3 will be implemented. This should allow to determine more physical coefficients representative of the overall emission of the lamp.

CONCLUSION AND PERSPECTIVES

This work proved promising for modelling infrared halogen lamps with an integrated reflector via an adapted radiosity method. A distributive function was proposed for modelling such a lamps emission and was found to give fair results for a given position. However, work still needs to be carried out taking into account the problem constraints in order to find values for the weighted coefficients that better represent the physical problem.

In this work an average convection coefficient was considered, further work should be carried out to better integrate convection into the problem.

Once the optimised distributive function is found, it will be integrated into a model of the IR oven and simulations will be carried out to compare experimental heating of an instrumented ABS sheet with the heating simulation.

FIGURE 1 ,

 1 FIGURE 1,. (a) Infrared halogen lamp with integrated reflector,(b) schematization

FIGURE 3 ,

 3 FIGURE 3, (a) Temperature field at heating time t=546 for U=160V. Profiles perpendicular and parallel to the lamp intersect at the point of maximum temperature, (b) Temperature versus time for the point of maximum temperature.

4

 4

FIGURE 4 ,

 4 FIGURE 4, (a) Simulation geometry, (b) Filament model

FIGURE 5 ,

 5 FIGURE 5, (a) Optimized surface temperature on the back face of heated ABS sheet at th. (b) Comparison between experimental data and simulated heating with optimized parameters at th.

TABLE 1 ,

 1 Thermal properties of ABS[START_REF] Dauphin | Amélioration des simulations thermiques dans les systèmes d'éclairage automobile[END_REF]

	Properties

TABLE 2 ,

 2 Optimized lamp parameters

	Properties	k1	k2	k3	k4 k5	k6	Tlampe (K)
	Initial value	1/6	1/6	1/6	1/6 1/6	1/6	1900
	Optimized value 0.63023 0.19598 0.10211	0	0	0.36477	1899.9

ACKNOWLEDGMENTS

These results were obtained under the research project "METEOR" at the IRT Saint Exupéry in collaboration with the ICA. We thank the industrial and academic members of the IRT who supported this project through their contributions, both financial and in terms of specific knowledge: Industrial members: ARKEMA, SD TECH, CHOMARAT,AIRBUS, HEXCEL, HUTCHENSON, PORCHER. And Academic members: CNRT Materiaux

We also thank the Commissariat Général aux Investissements and the Agence Nationale de la Recherche for their financial support in the Programme d'Investissement d'Avenir (PIA).