Professeur Amsini Sadiki

Sébastien Tanguy

François-Xavier Demoulin

L Germes Martínez

Réveillon and Thibaut Menard for the interest demonstrated towards my work and the discussions that have led to its improvement. Also, to Jorge Brändle de motta and Fabien Thiesset for their help in the post-processing (by the development of PyArcher) and valorization of the simulations. And the IT department for ensuring the correct operation of the computers.

The team formed by doctoral students of the laboratory, for the countless conversations during the "coffee break", where we were able to share our concerns and, of course, some laughs.

My family, who have always

Résumé

Les écoulements multiphasiques avec changements de phases sont présents de nombreuses applications industrielles telles que la combustion de spray, les échangeurs de chaleur, les réacteurs nucléaires. Par conséquent une meilleure compréhension des transferts de chaleur et de masse est essentielle pour l'optimisation de la consommation d'énergie de ces procédés et donc de la réduction des émissions polluantes.

Pour cela, des méthodes dédiées aux simulations numériques directes des écoulements diphasiques avec changement de phase sont présentées. Le formalisme est basé sur une méthode de représentation de l'interface avec une bonne conservation de la masse appelée CLSVOF, couplée avec une méthode de projection pour la résolution des équations de Navier-Stokes.

Ce formalisme peut être utilisé pour une large gamme de nombre de Mach et permet une formulation implicite des termes liés à l'acoustique. Cela nous permet de capturer simultanément plusieurs phénomènes tels que les changements de masse volumique de la phase gazeuse (due aux effets compressibles et mélange), le taux d'évaporation et la température, tous deux, non homogène au niveau de l'interface, etc. De plus les variations spatiales et temporelles de la pression permettent la considération de différentes structures de liquide et de gaz dans le domaine ; chacune avec différentes conditions thermodynamiques. De plus, nous avons fait une étude des procédés d'évaporation et de mélange turbulent en utilisant la configuration THI diphasique compressible. Dans cette étude, nous pouvons voir les effets du taux d'évaporation et du champ de vitesse turbulent sur les grandeurs statistiques de la fraction massique de vapeur. tours for a circular interface of radius 5 units [START_REF] Vinesh | On derivation and physical interpretation of level set method-based equations for two-phase flow simulations[END_REF]. . . . . . . . . . . . . 9 1.5 Interface representation with the Front tracking method [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. . . . . . . 1.6 Direct numerical simulations of a vaporizing isolated droplet with different methods for the interface capturing and computation of the liquid velocity field (a) [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] (b) [START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF] and (c) [START_REF] Reutzsch | A consistent method for direct numerical simulation of droplet evaporation[END_REF]. . . . . . . . . . . . . . 1.7 Direct numerical simulations of a vaporizing droplets with a lowmach formalism of the Navier-Stokes equations in the gas phase (a) [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] (b) [START_REF] Michael S Dodd | Analysis of droplet evaporation in isotropic turbulence through dropletresolved dns[END_REF]

List of Figures

Introduction

General context and objectives

Heat and mass transfer in multiphase flows has a vital role in many environmental and industrial processes, such as cloud and raindrop formation, spray cooling, nuclear reactors, spray drying, fuel combustion, etc. The latter example is present in a large diversity of internal combustion engines, e.g., the reciprocating engine of a car or an engine for electric power production and a continuous combustion engine of an airplane or a rocket. An example of the fuel spray and combustion on a piston of a compression ignition engine is shown in fig. 1a.

These processes share a common principle; the fuel is injected as a liquid into the combustion chamber. Then, it undergoes several morphological changes from the initial jet formation to the secondary break up to form small liquid structures. These structures interact with a high-temperature gas environment where temperature increases and evaporation starts. The next step is the combustion of the fuel vapor.

Here, evaporation depends significantly on the interaction between the liquid structures and the turbulent mixing of the fuel vapor. Additionally, the combustion efficiency is related to the mixing of the fuel vapor and oxidant. Consequently, understanding the evaporation and turbulent mixing processes could help optimize the design of combustion chambers and thus reduce fuel consumption and pollutant emissions.

Spray drying is another industrial process in which evaporation and turbulent mixing play an essential role and is widely used in the food industry, e.g., milk powder or instant coffee production (see fig. 1). Spray drying is a drying process to transform a liquid feed into a particulate dried product in a single operation [START_REF] Robert H Perry | Perry's chemical engineers' handbook[END_REF]. It increases the liquid's surface area by atomization and is followed by contact with high-temperature air. At this point, the liquid evaporates, and particles are formed while the droplet is still suspended; the final step is the recovery of the dried product. Like spray combustion, the efficiency of the liquid feed evaporation depends on the mixture of the vapor and dry air and the mean This work will focus on Direct Numerical Simulation (DNS) of the evaporation process on turbulent atomization regimes. In these simulations, the phase change occurs at the interface and temperatures below the boiling point. It is mainly controlled by gas density, temperature, and vapor mass fraction gradients at the interface. The vapor generated at the interface creates a velocity divergence that depends on the density ratio between the phases. Consequently, a Stefan flow is produced, which is a mass flow in the normal direction of the interface.

Even though many codes in the literature are dedicated to DNS evaporation, most of them focus on isolated droplets or bubbles and use an incompressible formalism. Consequently, the simulations are limited to open environments, and changes in density and pressure due to vapor production and mixing can not be considered. The latter is important for estimating the influence of Stefan flow on the system dynamics.

A compressible formalism should be implemented for a more realistic representation of the evaporation process in an atomization regime. There are several paths to overcome this task, e.g., Riemann solvers and pressure-based solvers. However, in simulations with a low-Mach number, the accuracy of the Riemann solvers tend to decrease. Hence, this work focus on a pressure-based solver coupled with a mass conservative interface capturing method. This has several advantages, including its compatibility with standard projection methods typically used in incompressible formalism. Also, a implicit formulation of the acoustic terms is allowed, reducing the time step restrictions. Since the method is based on a fully compressible formalism for multiphase flows, the spatial and temporal pressure changes are captured. This allows the simulations to handle many liquid and gas structures, each with its own density and thermodynamic pressure.

Furthermore, the method is suitable for a two-phase HIT configuration with phase change, similar to [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF][START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF]. In this work, we are interested in including the energy and species conservation equations in the formalism. This allows an estimated pressure, density, and vaporization rate dependent on temperature and vapor mass fraction simultaneously. Therefore, a more realistic atomization and evaporation processes model is described.

After the vaporization, the vapor is transported and mixed with the gas phase. In most industrial applications, the gas phase is submitted to a turbulent velocity field that significantly affects the vapor air mixture. This turbulent mixing can be defined as a three-stage process: entrainment, stirring and mixing ( [START_REF] Paul E Dimotakis | Turbulent mixing[END_REF]). The first stage corresponds to the vapor generation, and mixing is controlled by diffusion of the vapor layers close to the interface towards the gas phase. In the second stage, mixing is governed mainly by the turbulent velocity field and is characterized by the distortion of the interface issue of the interaction with the larger eddies in the flow. Finally, as the fluids become more homogeneous in the final state, the concentration gradients decrease until no more heat or mass transfer occurs. Most of the investigation dedicated to this phenomenon occurs in a purely gaseous environment or dispersed two-phase flows, leaving aside the influence of the vapor/interface interactions. This motivates us to use our compressible formalism to contribute to a better comprehension of turbulent mixing.

The work of this Ph.D. thesis is based on the in-house ARCHER code for 2D and 3D incompressible and compressible two-phase flows ( [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]). It has been developed in the Coria laboratory over the last two decades, and it was validated and applied in several configurations related to the atomization process.

Thesis outline

The main objective of this thesis is to develop a compressible formalism for the Direct Numerical Simulation (DNS) of two-phase flows with phase change. The first step in achieving this goal is to describe the challenges related to the interface representation in multiphase flows simulation, accompanied by some of the dedicated numerical methods proposed in the literature in chapter 1. This includes the Coupled Level-set/Volume of fluid method implemented in the inhouse Archer code ( [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]). In addition, the different strategies for studying the evaporation process (experimental, modeling, and DNS) are discussed. Moreover, we present the constraints associated with DNS of two-phase flows with phase change, such as the discontinuities of the velocity field at the liquid/gas interface.

In chapter 2, the modifications implemented to adapt the in-house ARCHER code for the simulation of incompressible two-phase flows with phase change are discussed. Among these changes, we have the additional term considered in the VOF transport equation to evaluate the liquid mass lost due to evaporation and the implementation of the extra interface jump conditions required to respect the mass, energy, momentum, and species conservation. Additionally, several numerical methods proposed in the literature for the treatment of the velocity field discontinuity are implemented and compared. Then, the method chosen as the most appropriate for this thesis objectives is selected. Lastly, the incompressible formalism is validated using the results obtained with the D 2 law and tested in a 2D convected cylinder.

Afterward, the incompressible formalism method is extended to a weaklycompressible formalism using the pressure-based method proposed by [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF]. To this end, additional terms are considered on the governing equations to account for the thermal dilatation and the compressibility effects in each phase. Also, the numerical method implemented is presented in chapter 3. Then, validation cases are performed for the conservation of the total mass in a closed environment configuration. In addition, a compressible HIT two-phase simulation with phase change demonstrates our formalism's ability to handle atomization regimes.

In chapter 4, our compressible formalism is used to analyze the evaporation process and the scalar turbulent mixing. Two simulations are performed, where the evaporation regime varies in a compressible two-phase HIT configuration. Moreover, the scalar turbulent mixing is analyzed using the temporal evolution of the first two moments and the PDFs shapes of the vapor mass fraction. Finally, the general conclusion and some prospects of this work are discussed in chapter 4.6.

Chap. 1 | State-of-the-art

Introduction

This chapter discusses the different approaches found in the literature to study phase change in two-phase flows. Special attention is given to the vaporization process at temperatures below boiling. First, we discuss some of the most popular numerical methods for capturing interfaces in multiphase flows simulations. Then, we turn to experimental works on temperature and evaporation measurements of isolated and dispersed droplets. Finally, DNS analysis of evaporation in incompressible and compressible formalisms is presented.

Three main formalisms simulate two-phase flows: RANS, LES, and DNS. In RANS (Reynolds Averaged Navier-Stokes Simulations), the governing equations are averaged, allowing larger spaces to be simulated with relatively low computational cost. The disadvantage of this kind of simulation is that information about the smallest length scales is lost because only the averaged velocity field is represented. In the Large Eddy Simulation (LES), a spatial filter is applied to the Navier-Stokes equations, which means that only a range of length scales of the flow are resolved. Consequently, the effects of the smallest scales on the flow have to be modeled. Modeled terms become very important in some situations, e.g., near-wall interactions, reactive flows, multi-phase flows, etc. A detailed discussion about this kind of simulation can be found in [START_REF] Wendt | Computational fluid dynamics: an introduction[END_REF].

In the Direct Numerical Simulation (DNS), the governing equations are solved without averaging or filtering. Hence, all the temporal and spatial scales have to be solved, i.e., the mesh size has to be smaller than the Kolmogorov length scale. This restriction limits simulations to academic research because, for the time being, it is too computationally expensive to represent multiscale configurations.

Apart from choosing a model for solving the Navier-Stokes equation, an essential part of the simulation of multiphase flows is the interface's representation. This work will focus on the one-fluid framework, i.e., only one system of governing equations is solved for every phase involved, where the interactions between the phases depend on the method used for capturing the evolution of L. Germes Martínez DNS of phase change June 28, 2022 the interface.

Interface representation

One of the main challenges of simulating the two-phase flows is the interface representation. According to [START_REF] Mirjalili | Interface-capturing methods for two-phase flows: An overview and recent developments[END_REF], among the difficulties in developing a method for interface representation, we have: a) ensuring conservation of mass, momentum, and energy throughout the domain, b) treatment of discontinuous variables at the interface, especially for large density ratios, c) handling complex topologies and scale separation, d) robustness for realistic flow simulation and e) accurate implementation of surface forces. Fig. 1.1 shows the different models and their most popular methods for the interface. In this work, we will focus on the Figure 1.1: Classification of numerical methods for for the interface representation in two-phase flows [START_REF] Mirjalili | Interface-capturing methods for two-phase flows: An overview and recent developments[END_REF].

one-fluid model.

VOF methods

In the Volume-of-Fluid (VOF) methods, a phase indicator function (H) is used to identify if the reference fluid is present in a particular domain location. From (H) we can define the volume fraction as:

C k (t) = 1 V Ω H(x, t)dV
where V is the volume of the kth cell. The volume fraction takes the value of 1 or 0 in the cells filled in one phase or the other. The mixed cells (cells that contain both phases) take fractional values and consequently contain a piece of the interface; for example, in fig. 1.2, the cells containing only the 'material 1' have values of 1; for the 'material 2', the cells have the values of 0, and the cells in between have fractional values. A primary advantage of the VOF method is its capacity to conserve the quantity of volume in the domain. However, computing the geometrical interface properties (curvature and normal vectors) with the volume of fluid as the only information is not straightforward. First, one must reconstruct the interface based on the liquid volume fraction. In the literature, various approaches have been developed for this end. [START_REF] William | Slic (simple line interface calculation)[END_REF] presented the SLIP (Simple Line Interface Calculation) method for the geometric approximation of fluid interface. In this SLIC method, the interface is advected with a time splitting scheme. For example, in a 2D domain, for the advection in the horizontal direction, the cell is divided by a vertical line depending on the volume of fluid of the cell to the left or the right. Then, the same process is repeated with a horizontal line for the advection in the vertical direction. [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] proposed a similar method where the interface is also represented with vertical and horizontal lines. [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] first compute the normal interface vector using the volume of fluid of the neighboring cells and estimate the orientation of the interface. Then, the interface is advected in each direction, keeping the original orientation of the interface. Although, various attends were proposed for further improvement of the SLIC method [START_REF] Joel | Flame advection and propagation algorithms[END_REF][START_REF] Rider | Stretching and tearing interface tracking methods[END_REF]. In addition to distorting the interface, the SLIC methods generate considerable among of "flotsam" and "jetsam," where pieces of the interface get separated in an unphysical way ( [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]).

L. Germes Martínez

DNS of phase change June 28, 2022

A complete approach is the PLIC (Piecewise Linear Interface Calculation) method for the interface representation initially proposed by [START_REF] Rb Debar | Fundamentals of the kraken code[END_REF][START_REF] David | An interface tracking method for a 3d eulerian hydrodynamics code[END_REF]. Here the interface is represented by straight lines (or segments in 3D) in the mixed cells. The main difference with the SLIC methods is that the lines can take any arbitrary orientation concerning the axis. To determine the direction of the line or segment, the normals vectors to the interface has to be computed. Various methods for estimating the interface normals can be found in the literature [START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerian-lagrangian advection[END_REF][START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry[END_REF][START_REF] Edward | Second-order accurate volume-of-fluid algorithms for tracking material interfaces[END_REF]. Generally, the volume of fluid of the current cell and those of the neighboring cells are considered for the computation of the normals. This allows a more realistic representation of the interface, improving mass conservation. In fig. 1 After the reconstruction of the interface, the next step is to choose a method for advecting the volume of fluid. In the literature, two main types of methods can be found for the advection of the interface: split methods and non-split methods. The split methods are based on an operator-splitting numerical scheme with a series of one-dimensional operations, which are combined in the last step [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry[END_REF][START_REF] Gabriel | Conservative volume-of-fluid method for free-surface simulations on cartesian-grids[END_REF]. In the un-split method, the interface is advected in a single step [START_REF] López | A volume of fluid method based on multidimensional advection and spline interface reconstruction[END_REF][START_REF] Owkes | A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (vof) method[END_REF].

Level-set method

The level set (LS) method was first introduced in the context of multiphase flows by [START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. The method consists of the advection of a signed-distance function, representing the shortest distance from the interface to each point in the domain. The level set function is positive in one fluid and negative in the other, making the interface between the fluids the zero level curve. The level set interface representation method and the level set contours are shown in fig. 1.4. Compared to the VOF method, the LS methods have many advantages: an accurate and relatively easy computation of the interface normals and curvature estimated as the normalized gradient of the LS function and with the divergence of the interface norm, respectively. An accurate estimation of the derivatives is possible due to the continuous nature of the LS function. After an advection step, the level set function loses its signed-distance properties. To solve this problem, [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] proposed a reinitialization step based on the resolution of a partial differential equation (PDE) at each time step. Some improvements dedicated to the reinitialization step are available in the literature. For example, [START_REF] Peng | A pde-based fast local level set method[END_REF] proposed a replacement for the sign function used by [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] to decrease the artificial motion of the interface during the reinitialization. The reinitialization can also be performed by solving the eikonal equation, which imposes that the norm of the gradient of the signed function is unity. For solving this equation, several techniques can be found in the literature; for example, [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] developed the FFM (Fast marching method) and the FSM (Fast Sweeping method) proposed by [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF]. An approach for getting rid of the reinitialization has been submitted by [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF]. Their method provides a way to update the interface with an extended velocity field, which allows for keeping the LS function's signed-distance properties. In addition, [START_REF] Sabelnikov | Modified level set equation and its numerical assessment[END_REF] proposed a modification of the level set equation to make the reinitialization steps unnecessary. To this end, a source term is added to the level set equation, the construction of which implies that the eikonal equation is automatically satisfied.

A well-known problem with the level set method is the lack of mass conservation. Over the last two decades, numerous attempts to overcome this issue have been proposed. [START_REF] Russo | A remark on computing distance functions[END_REF][START_REF] Du Chéné | Second-order accurate computation of curvatures in a level set framework using novel highorder reinitialization schemes[END_REF] proposed higher-order accuracy schemes for the reinitialization equation of [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], improving the computation of the interface normals and curvatures. Some authors have used the fact that the mass conservation errors decrease when the mesh is refined. For example, [START_REF] Herrmann | A balanced force refined level set grid method for twophase flows on unstructured flow solver grids[END_REF] followed the interface location in an auxiliary, high-resolution equidistant Cartesian grid. Additionally, for solving the mass conservation problem, [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] proposed the conservative level-set method (CLS). With this approach, [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] relaxes the sharp interface assumption by using a smeared out Heaviside function based on the signed distance function. In the CLS method, the interface is represented as the 0.5 iso-contour of the Heaviside function. It is coupled to the momentum equation via the definition of the discontinuous variables. Despite the interface's diffusive representation, straightforward computation of the interface normals and curvatures is still possible. Then, [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] presented the accurate conservative level set (ACLS) method which combined the CLS proposed by [START_REF] Olsson | A conservative level set method for two phase flow[END_REF] and the ghost fluid method. Also, [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF] improved the computation of the interface normals by using the reconstructed level set function obtained by an FMM. [START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF] modified the reinitialization equation of the ACLS to remain conservative and reduce interface deformation. A review of some advances and application of the level set-method can be found in [START_REF] Gibou | A review of level-set methods and some recent applications[END_REF].

Front Tracking methods

In the Front tracking method, the interface is represented by connected markers advected with an imposed velocity. When using markers for the interface representation, the deformations, breakups, and coalescent are not handled explicitly as in the other methods. According to [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF], the implementation of this kind of method requires overcoming several challenges, such as 1) data structure used to describe the front, 2) the maintenance of the front as it is deformed by the flow, 3) the front interaction with the fixed grid, where the Navier-Stokes equation are solved and 4) the topology changes (interaction between two interfaces).

One of the first to develop this method for flows with shocks and interface was [START_REF] Glimm | A numerical method for two phase flow with an unstable interface[END_REF]. In this method, a sharp interface representation is achieved by extrapolating the variables at one side of the front to the grid point at the other side. Fig. 1.5 shows the interface representation with the front tracking method. [45] developed a method to simulate multiphase flows with a front for identifying each fluid, where a single system of governing equations is used for the whole field. [START_REF] Salih | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF] use a stationary grid for the flow field and an unstructured grid for the interface that moves through the fixed grid. The discontinuous variables are smoothed across the interface by introducing an indicator function. According to the authors, the main drawback of this method is to handle the surface force term when interfaces are close to each other. This method was fatter developed by [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]. For improving the representation of the front, [START_REF] Holavanahalli S Udaykumar | Multiphase dynamics in arbitrary geometries on fixed cartesian grids[END_REF] used the 'cut cell' method for multiphase flows simulations, where the main idea is to modify the grid near the interface in a way that the grid lines coincide with the interface.

Hybrid methods

Several hybrid methods can be found in the literature to benefit from a specific feature of more than one method for interface representation simultaneously. For example, as a solution for the mass conservation problem when using the level set methods, [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF] proposed the HPLSM (Hybrid Particle Level set Method). [START_REF] Enright | A hybrid particle level set method for improved interface capturing[END_REF] implanted Lagrangian marker particles in the interface and used them to correct the level set function when the particles cross the interface. This way, maintaining the geometric properties of the level set method while improving the mass conservation errors in regions of high curvature. [START_REF] Aulisa | Interface reconstruction with least-squares fit and split advection in three-dimensional cartesian geometry[END_REF] proposed a combined Eulerian-Lagrangian scheme for the advection of the interface for the VOF method. [START_REF] Liu | A coupled phase-field and volume-of-fluid method for accurate representation of limiting water wave deformation[END_REF] proposed a coupled phase-field and volume of fluid method for studying water waves during breaking. This method takes advantage of the phase-field method for the accurate representation of the dynamic processes and reduces the overall computational cost of the VOF method.

In the literature, a popular hybrid method is the coupled level set volume of fluid (CLSVOF) process first proposed by [START_REF] Bourlioux | A coupled level-set volume-of-fluid algorithm for tracking material interfaces[END_REF][START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF]. This method owes its success because it takes advantage of both strategies. With the CLSVOF, the sharp representation of the interface and the exemplary description of the geometrical properties of the level set method are combined with the mass conservation properties of the VOF method. For coupling the level set method with the volume of fluid method, first, the volume fraction must be defined in terms of a level set function. For this, a reconstructed level set is defined. Since the volume fraction of a grid cell must be the same, regardless of the method used to estimate it. The volume fraction value obtained with the level set function is corrected with the volume fraction obtained with the VOF method. For the advection of the interface, for both VOF and LS functions, the coupled second-order conservative operator split advection scheme, described in [START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF], is used. Some applications for the CLSVOF method in the literature are: [START_REF] Tomar | Numerical simulation of bubble growth in film boiling using a coupled level-set and volume-of-fluid method[END_REF] used a CLSVOF based method for studying the bubble growth in film boiling; [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF] implemented the CLSVOF method proposed by [START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF] for simulating the primary breakup of a turbulent liquid jet; and [START_REF] Wang | A coupled level set and volume-of-fluid method for sharp interface simulation of plunging breaking waves[END_REF] used a CLSVOF method for the simulation of plunging breaking waves.In [START_REF] Boniou | Comparison of interface capturing methods for the simulation of two-phase flow in a unified lowmach framework[END_REF] a comparison of CLSVOF with 3 other methods for interface capturing is performed, concluding that CLSVOF is the most promising method due to its mass-conserving properties.

More recently, a hybrid moment of fluid/level set (HyMOFLS) approach for the interface reconstruction was proposed by [START_REF] Asuri Mukundan | A hybrid moment of fluid-level set framework for simulating primary atomization[END_REF]. In the HyMOFLS method, the moment of fluid (MOF) method initially proposed by [START_REF] Dyadechko | Reconstruction of multimaterial interfaces from moment data[END_REF] is used for capturing the under-resolved liquid structures, and the CLSVOF method is used for capturing the resolved liquid structures. This is possible because the MOF method uses the centroid in each cell to improve the interface reconstruction, increasing its accuracy in under-resolved interface areas. To identify the under-resolved liquid structures, [START_REF] Asuri Mukundan | A hybrid moment of fluid-level set framework for simulating primary atomization[END_REF] use the interface resolution quality (IRQ) approach developed by [START_REF] Canu | Curvature-based interface resolution quality (irq) indicator to assess simulation accuracy[END_REF]. Then, [START_REF] Asuri Mukundan | Detailed numerical simulations of primary atomization of airblasted liquid sheet[END_REF] use the HyMOFLS method for studying the primary atomization of air blasted liquid sheet. Among the analysis presented by [START_REF] Asuri Mukundan | Detailed numerical simulations of primary atomization of airblasted liquid sheet[END_REF], we can find the effect of the velocity field on the atomization characteristics, the occurrence of the secondary atomization, coalescences, and a detailed study of the information on the near-field atomization characteristics.

Phase change 1.3.1 Experimentation

In general, experiments on heat transfer, evaporation and combustion of droplets are mainly focused on stationary droplets supported by fibers ( [START_REF] Mikami | Microgravity experiments on flame spread along fuel-droplet arrays using a new droplet-generation technique[END_REF][START_REF] Pa Strizhak | Heating and evaporation of suspended water droplets: Experimental studies and modelling[END_REF][START_REF] Chauveau | An analysis of the d2-law departure during droplet evaporation in microgravity[END_REF]) or droplets in tandem ( [START_REF] Kristyadi | Monodisperse monocomponent fuel droplet heating and evaporation[END_REF][START_REF] Maqua | Monodisperse droplet heating and evaporation: experimental study and modelling[END_REF]). In the first case, one of the main difficulties in estimating the role of the fiber in the temperature distribution of the droplet. This problem is analyzed by [START_REF] Chauveau | An analysis of the droplet support fiber effect on the evaporation process[END_REF], which observed that most of the fiber techniques used in the literature, which use relatively large fibers, enhance the evaporation rate of droplets due to the increased heat transfer by conduction through the fiber. On the topic of combustion, [START_REF] Mikami | Microgravity experiments on flame spread along fuel-droplet arrays using a new droplet-generation technique[END_REF] proposed a droplet-array generation technique to study flame propagation in microgravity and presented results to demonstrate the effects of the fiber during the experiments. Moreover, to access the validity of the D 2 law, [START_REF] Pa Strizhak | Heating and evaporation of suspended water droplets: Experimental studies and modelling[END_REF] implemented a technique to reduce the impact of the fiber on the heat transfer between the droplet and its environment and performed suspended droplet experiments in stagnant hot gas environments. Also, [START_REF] Chauveau | An analysis of the d2-law departure during droplet evaporation in microgravity[END_REF] analyzed numerically and experimentally the heating and evaporation of droplets in suspension and proposed an analytical model for the droplet temperature distribution. Another interesting application is the explosive breakup of an oil-water emulsion which is also studied with a suspended droplet [START_REF] Strizhak | Evaporation, boiling and explosive breakup of oil-water emulsion drops under intense radiant heating[END_REF][START_REF] Antonov | Impact of holder materials on the heating and explosive breakup of two-component droplets[END_REF][START_REF] Dv Antonov | Breakup and explosion of droplets of two immiscible fluids and emulsions[END_REF]. One example is the work done by [START_REF] Strizhak | Evaporation, boiling and explosive breakup of oil-water emulsion drops under intense radiant heating[END_REF], which studied the vaporization of oil-water drops at high temperatures. Here, the vaporization rate is studied as a function of the gas temperature, and the effects of the concentration of organic compounds on the droplet are also presented. The conditions for the boiling and subsequent explosion of the droplet are also discussed. In addition, the effects of a turbulent flow on the evaporation rate of a single suspended droplet have been studied by [START_REF] Birouk | Current status of droplet evaporation in turbulent flows[END_REF][START_REF] Verwey | Experimental investigation of the effect of droplet size on the vaporization process in ambient turbulence[END_REF][START_REF] Verwey | Fuel vaporization: Effect of droplet size and turbulence at elevated temperature and pressure[END_REF].

In the case of droplets in tandem, [START_REF] Castanet | Evaporation of closely-spaced interacting droplets arranged in a single row[END_REF] studied the evaporation of a line of droplets in a high-temperature chamber. They observed a strong influence of the distance between droplets on the evaporation rate; this shows that the vapor generated by one drop can impact the evaporation of a subsequent drop if there is not enough space for the vapor to dissipate. This scenario is common in dense areas of a spray. [START_REF] Deprédurand | Heat and mass transfer in evaporating droplets in interaction: Influence of the fuel[END_REF] used a similar configuration of monodisperse droplets to study the effects of the nature of the fuel on the Nusselt and Sherwood numbers. In addition, [START_REF] Castanet | Heat convection within evaporating droplets in strong aerodynamic interactions[END_REF] studied heat transfer in a monodisperse droplet configuration in a diffusion flame fueled by the vapor generated by droplet evaporation.

Another interesting approach for the experimental study of evaporation droplets is the Lagrangian measurement of the trajectories and diameters of the droplets. To this end, [START_REF] Chareyron | Testing an in-line digital holography 'inverse method'for the lagrangian tracking of evaporating droplets in homogeneous nearly isotropic turbulence[END_REF] implemented an in-line digital holography technique to obtain Lagrangian statistics of evaporating droplets in a turbulent flow. This technique was used by [START_REF] Marié | Lagrangian measurements of the fast evaporation of falling diethyl ether droplets using in-line digital holography and a high-speed camera[END_REF] to study the relatively fast evaporation of diethyl ether in a falling droplet configuration. The authors state that this simplified configuration allows the results of droplet velocities and evaporation to be compared with existing theoretical models, demonstrating the accuracy of this technique. In addition, [START_REF] Marié | Digital holographic measurement of the lagrangian evaporation rate of droplets dis-persing in a homogeneous isotropic turbulence[END_REF] applied the same technique to study the evaporation of the diethyl ether in a homogeneous and nearly isotropic turbulence flow. They observed an enhancement of the evaporation rate in certain droplets with different trajectories, indicating the influence of the velocity fluctuations in the droplet evaporation time. More recently, [START_REF] Méès | Statistical lagrangian evaporation rate of droplets released in a homogeneous quasiisotropic turbulence[END_REF] presented Lagrangian statistical results for the evaporation rate in a similar configuration. And thus, supporting and quantifying the previous observations made by [START_REF] Marié | Digital holographic measurement of the lagrangian evaporation rate of droplets dis-persing in a homogeneous isotropic turbulence[END_REF].

Even though significant advancements regarding heat and mass transfer measurement techniques have been made in the literature, e. g., Interferometric Particle Imaging [START_REF] Av Bilsky | Interferometric technique for measurement of droplet diameter[END_REF], Laser-Induced Fluorescence [START_REF] Yerbol K Akhmetbekov | Planar fluorescence for round bubble imaging and its application for the study of an axisymmetric two-phase jet[END_REF] and Phase Rainbow Refractometry [START_REF] Wu | Phase rainbow refractometry for accurate droplet variation characterization[END_REF], it is still very challenging to quantify the vaporization rate in a dense spray where many liquid structures interact between each other. Several numerical studies have been developed in the past years to describe further the phase change process in dense spray or turbulent liquid jets.

Modeling

To simulate the evaporation process of a droplet, numerous mathematical models exist. Following [START_REF] Sazhin | Advanced models of fuel droplet heating and evaporation[END_REF], we can subdivide the models into three different groups: 1) hydrodynamic models, which only describe the vapor transport from the immediate proximity of the interface to the ambient gas, assuming that the vapor around the interface is saturated; 2) Kinetic models, that take into account the detachment of liquid molecules from the interface based on the kinetic Boltzmann equations; 3) Molecular dynamics models, also consider the molecular interaction at the interface but the use a molecular model instead. In this work, we are only going to analyze the hydrodynamics models.

The basis for the theory droplets evaporation in a gaseous medium date from 1877 and was introduced by Maxwell [START_REF] Sazhin | Advanced models of fuel droplet heating and evaporation[END_REF]. In its model, the vaporization rate of and static droplet is only governed by the diffusion process:

ṁd = -4πR d D v (ρ v s -ρ v ∞ ) (1.1)
where D v and ρ v is the diffusion coefficient and the density of the vapor. Ignoring the effect of the flow around the droplet due to the difference between the densities of each phase (Stefan flow). Later, improvement to this approach were proposed to add the effect of the Stefan flow to the model, under the assumption that the total density of the gas in the vicinity of the droplet remains constant

(ρ t = ρ v + ρ g = const): ṁd = -4πR d D g ρ t ln (1 + B M ) (1.2)
where B M = ρvs -ρv∞ ρgs is the Spalding mass number.

Another alternative is to estimate the vaporization rate based on the energy balance across the interface of the droplet

ṁd = -4πk g R d c pv ln (1 + B T ) (1.3)
where B T = cpv (Tg-T Γ )

L(T Γ )+ q l ṁd
, T Γ is the interface temperature and c pv is the heat capacity of the vapor. Even though this kind of models have been implemented in several CFD codes there is still room for its further improvement [START_REF] Ss Sazhin | New approaches to numerical modelling of droplet transient heating and evaporation[END_REF].

A well-known improvement of the classical model was made by [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF]. With the interest of developing a simple but efficient model that can be used in spray combustion models. To take into account the convective transport caused by the relative motion of the droplet, [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF] use the 'film theory'; where the mass and heat transfer at the interface is modeled by introducing the concept of gas films with a constant thickness:

δ T 0 = 2R d Nu o -2 (1.4) δ M0 = 2R d Sh o -2 (1.5)
where Sh 0 and Nu 0 are the Sherwood and Nusselt numbers, respectively. An option for estimate Sh 0 and Nu 0 are the Frossling correlations. It is important to highlight that another correlations, probably more suitable for the case considered, can be used for this end.

In [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF], two models for considering the internal recirculation of the moving droplet in the temperature evolution of the interface are compared: 1) the 'extended model', where the instantaneous velocity film inside the moving droplet is approximated by the Hill spherical vortex solution and 2) 'effective conductivity model', which use an effective thermal conductivity factor to account for the variation in the heat transfer due to the internal recirculation of the droplet. For both approaches, the interface temperature is considered uniform but time varying. In the extended film model, the droplet vaporization rate can be estimated with:

ṁd = 2πR d ρ g D g Sh * ln (1 + B M ) (1.6) ṁd = 2πR d k g c pv Nu * ln (1 + B T ) (1.7)
where Sh * and Nu * are the modified Sherwood and Nusselt numbers, determined by:

Sh * = 2 + Sh 0 -2 F M (1.8)
Nu * = 2 + Nu 0 -2 F T (1.9)
and F T and F M are correlation factors which depend of heat and mass transfer properties of the system, respectively.

Here, the gas phase is considered ideal, and the Raoult law estimates the molar and mass fraction at the interface. Additionally, other assumptions for developing the model are: the gas phase heat and mass transfer are in a quasisteady state, the pressure drop in the gas is negligible, and the thermophysical properties of both phases are constant and computed following the '1/3 rule', recommended by [START_REF] Hubbard | Droplet evaporation: Effects of transients and variable properties[END_REF]. [START_REF] Abramzon | Droplet vaporization model in the presence of thermal radiation[END_REF][START_REF] Abramzon | Convective vaporization of a fuel droplet with thermal radiation absorption[END_REF] extended this model by taking into account the contribution of the thermal radiation and the temperature dependence of the liquid properties.

Direct numerical simulation of phase change

Incompressible formalism

We consider that the simulation is a DNS when the vaporization rate is computed directly from the temperature or the species fields without using any model. Depending on the physical phenomena of interest, e.g., for simulating the evapora-tion process, the vaporization rate is based on the gradient of the species near the interface. In CFD codes, a common strategy found in the literature is implementing an interface capturing method and modifying it to consider the phase change. Usually, this is done implicitly by computing an interface regression speed based on the vaporization rate and the density of the liquid and adding it to the velocity field used for the interface advection scheme. Or explicitly, using a volumetric sink term in the advection equation for the liquid volume fraction.

Furthermore, a discontinuity in the velocity field occurs across the interface. Note that most numerical methods dedicated to the convective transport of the scalars in DNS of two-phases flows are developed initially to handle continuous and divergence-free velocity fields. Hence, they cannot use this velocity field directly. Otherwise, unphysical deformation of the interface will occur due to an oscillatory behavior of the numerical scheme, which will lead to errors in the mass conservation. Several strategies have been proposed to overcome this issue. Most numerical approaches share the same idea: create a divergence-free velocity field dedicated to the convective transport of the interface capturing methods. This section presents a discussion about the simulation methods of vaporization simulation.

Beginning with the volume of fluid method, [START_REF] Calimez | Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs[END_REF] is one of the first authors that coupled the phase change with the VOF method for DNS of reacting twophases flows. Even though its work was a significant advancement for the field at the moment, the relatively strong assumptions made during the development of the method limited it to idealized configurations. This is probably due to the computational constraints of time. Then, in the domain of boiling, [START_REF] Samuel | A volume of fluid based method for fluid flows with phase change[END_REF] proposed a method for phase change using a VOF method for interface capturing. One of the contributions of [START_REF] Samuel | A volume of fluid based method for fluid flows with phase change[END_REF] was improving the computation of the discontinuous heat fluxes at the interface, highlighting the fact that heat fluxes near the interface are of great importance for the mass transfer calculation. The heat and mass transfer of a droplet impinging onto a hot surface was studied by various authors using this approach [START_REF] Nikolopoulos | A numerical investigation of the evaporation process of a liquid droplet impinging onto a hot substrate[END_REF][START_REF] Strotos | Numerical investigation on the evaporation of droplets depositing on heated surfaces at low weber numbers[END_REF]. For these studies, an important issue is the treatment of the velocity discontinuity at the interface.

A method to compute the liquid and gas velocities near the interface was proposed by [START_REF] Schlottke | Direct numerical simulation of evaporating droplets[END_REF]. They emphasized the importance of these velocities in configurations where the Stefan flow is relatively high. Here, the liquid and gas velocity fields are computed considering the volume source term in the continuity equation dedicated to the Stefan flow. An iterative algorithm is implemented to compute and distribute the error in the divergence of the continuity equation to the neighboring cell. This error is due to the difference between the volume and mass averaged velocities at the interface. Another important point of this work is the computation of the vaporization rate using the concentration of the vapor at the interface, allowing us to consider the phase change in temperature below the boiling temperature.

Recently, [START_REF] Reutzsch | A consistent method for direct numerical simulation of droplet evaporation[END_REF] improved the method of [START_REF] Schlottke | Direct numerical simulation of evaporating droplets[END_REF] by adding a correction term for the volume source term in the continuity equation which takes into account the difference between the volume and mass averaged velocities. Consequently, the iterative algorithm or the distribution of the error is no longer needed. A different approach was proposed by [START_REF] Malan | Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change[END_REF]. In their method, [START_REF] Malan | Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change[END_REF] used the liquid velocity (u l ) for the interface advection, which is computed from the velocity field (u) and an intermediary velocity ( ũ) as follow:

u l = u - ũ
ũ is computed in a similar form than u, using an intermediary variable (p):

∇ • ∆t ρ n+1 ∇p = ṁ 1 ρ g - 1 ρ g δ Γ (1.10)
where the RHS of eq. 1.10 is similar to the phase change source term in the Poisson equation for the pressure but with the opposite sign. u l is divergence-free at the interface by construction, allowing the use of regular advection schemes for the VOF variable.

A similar idea has been proposed by [START_REF] Scapin | A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows[END_REF]. In their paper, [START_REF] Scapin | A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows[END_REF] used an FFTbased two-fluid Navier-Stokes solver, which results in a reduction of the computation time. An algebraic VOF method for interface capturing. Also, they obtained a u l from a potential computed in a similar way than p in eq. 1.10. Since they use an FFT-based solver for the constant-coefficient Poisson equation, the density difference between the phases is not considered during the computation of the LHS of eq. 1.10. Another interesting approach was proposed by [START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF]. The liquid velocity field in [START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF] is computed by extending the velocity in the liquid phase into the gas phase using Aslam's extension method ( [START_REF] Tariq D Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF]). As mentioned by the authors, the extended velocity field (u * l ) is not divergence-free, leading to erroneous variations in the liquid volume. To overcome this problem, u * l is projected into its divergence-free part using:

u l = u * l + ∇ • W (1.11)
where W is a potential derived from the Helmholtz-type equation,

aW + ∇ • W = ∇ • u * l
where a is zero for all liquid-containing cells and cells that are within three grid cells of the interface. A drawback of this method is that the Aslam's extension implies an important increase in the computational cost.

Among the advancements done in this matter, we can find the method of [START_REF] Bures | Direct numerical simulation of phase change in the presence of non-condensable gases[END_REF]. Which is devoted to the direct simulation of phase change, taking into account the effect of the mixture vapor/non-condensible-gas on the heat and mass transfer processes. This last method was improved by [START_REF] Bures | Direct numerical simulation of evaporation and condensation with the geometric vof method and a sharp-interface phase-change model[END_REF]. They incorporated the method of [START_REF] Malan | Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change[END_REF] for building the liquid velocity field for advecting the interface.

Continuing with the level-set methods for the simulation of phase change. [START_REF] Son | Numerical simulation of film boiling near critical pressures with a level set method[END_REF] is one of the first employing the level-set method to simulate thin film boiling simulations. Here, the discontinuities at the interface are smeared out by considering the interface to have a finite thickness. Similar simulations were presented by [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible navierstokes equations with phase change[END_REF], which proposed a level set based sharp interface method for the direct simulation of phase change. For treating the discontinuities at the interface, [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible navierstokes equations with phase change[END_REF] used a ghost fluid method, and the projection method proposed by [START_REF] Duc Q Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF]; devoted initially to propagating a flame front in an incompressible fluid. Here, ghost velocities are defined in both phases, and the velocity jump condition due to the phase change is used for this purpose, following:

u ghost l = u g + s d ρ l -ρ g ρ g (1.
12)

u ghost g = u l -s d ρ l -ρ g ρ g (1.13)
where u ghost l and u ghost g are the liquid and gas velocity extension, respectively. And s d is the interface regression speed.

For improving mass conservation, [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] proposed a modified version of this method. [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] imposed the divergence-free condition to the liquid velocity field and its extension in the gas phase; this is done by solving an additional Poisson equation for a potential function. Even though this modification added additional cost to the simulations, nonnegligible improvement was presented in the mass conservation. Also, this method can consider more physical phenomena, such as evaporation and non-homogeneous vaporization rate. Later on, [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF] modify the method of [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] for simulating the phase change to study the Leiden- frost effect. The process of [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF] can consider both evaporation and boiling in the same liquid structure while keeping the temperature continuous across the interface. [START_REF] Urbano | Direct numerical simulation of nucleate boiling in zero gravity conditions[END_REF] used the method for boiling presented by [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF] for studying nucleate boiling of droplets in zero gravity conditions. [START_REF] Popescu | On the influence of liquid/vapor phase change onto the nusselt number of a laminar superheated or subcooled vapor flow[END_REF] used a level-set/Ghost fluid method for studying the effect of phase change in the heat transfer. Recently, [START_REF] Boniou | On the numerical simulation of evaporating two-phase flows using sharp interface capturing methods[END_REF] preformed simulation of two-phase flows with phase change using the level set method and the VoF method, to compare their accuracy in estimating the vaporization rate. In fig. 1.6 are shown examples of simulation of vaporizing isolated droplets using some of the methods mentioned above.

Concerning the front tracking methods, [START_REF] Juric | Computations of boiling flows[END_REF] proposed a front tracking approach for the simulation of film boiling. In [START_REF] Juric | Computations of boiling flows[END_REF], a delta function formulation is used to add the source terms associated with the surface tension, mass transfer, and latent heat. This method was later improved by [START_REF] Esmaeeli | Computations of film boiling. part ii: multi-mode film boiling[END_REF]. Then, the latter method was used in [START_REF] Esmaeeli | Computations of film boiling. part ii: multi-mode film boiling[END_REF] for studying the heat and mass transfer of film boiling on a horizontal periodic surface. [START_REF] Esmaeeli | A front tracking method for computations of boiling in complex geometries[END_REF] implemented a front tracking/finite difference method for simulations of boiling in complex geometries. [START_REF] Irfan | A front tracking method for direct numerical simulation of evaporation process in a multiphase system[END_REF] studied both the temperature and vapor gradient driven vaporization rate with a front tracking method.

Compressible formalism

Until now, the studies mentioned above use an incompressible formalism to resolve the Navier-Stokes equations. This assumption is no longer suitable for simulating multiphase flows with phase change in an enclosed environment or with encapsulated gas structures such as atomization configurations. In these kinds of studies, an outflow boundary condition is used to evacuate the divergence created at every iteration in the gas. Consequently, the gas density and the thermodynamic pressure remain constant, leading to an overestimation of the Stefan flux in the case of n-decane vaporization in air. A weakly compressible approach should be developed to account for the increase in gas density and pressure due to vaporization and/or thermal dilatation.

Studies capable of describing this type of problem are scarce in the literature. An example is the work of [START_REF] Wang | Vaporization of liquid droplet with large deformation and high mass transfer rate, ii: Variable-density, variable-property case[END_REF] which proposed a method for modeling the evaporation of two-phase flows while solving the Navier-Stokes equations under a low Mach number limit. In this model, the mass and heat transfer occurs in a sink layer in the liquid and a source layer in the gas phase. More recently, [START_REF] Michael S Dodd | Analysis of droplet evaporation in isotropic turbulence through dropletresolved dns[END_REF] studied the evaporation of an isolated droplet and the effects of the initial liquid volume fraction on the evaporation rate in a homogeneous isotropic turbulence configuration. A low Mach number formulation is used to solve the Navier-Stokes equations in the gas phase while maintaining the incompressible assumption in the liquid phase. Also, DNS results were validated with the measurement of nheptane droplets evaporating at high temperature and pressure obtained by [START_REF] Verwey | Fuel vaporization: Effect of droplet size and turbulence at elevated temperature and pressure[END_REF]. Among the observations of [START_REF] Michael S Dodd | Analysis of droplet evaporation in isotropic turbulence through dropletresolved dns[END_REF], the evaporation rate increase with the Reynolds number. One possible explanation is that turbulent flow removes vapor in the region near the interface, increasing the vapor fraction gradient. Also, the evaporation rate decreases with increasing droplet diameter, which can be explained by the interaction of the droplets with vapor produced by neighboring droplets.

[10] proposed a low Mach number formulation for studying the evaporation of droplets in homogeneous shear turbulence. Here, the evaporation framework based on the arithmetic VOF method for the interface capturing proposed by [START_REF] Scapin | A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows[END_REF] is extended to account for the gas density changes with the low-mach formalism proposed by [START_REF] Barba | An interface capturing method for liquid-gas flows at low-mach number[END_REF]. In their study, [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] observed, among other things, an increase in the evaporation rate in the turbulent environment compared to the evaporation rate obtained in a stagnant gas. Another approach was proposed by [START_REF] Poblador | Volume of fluid method for low-mach-number compressible supercritical liquid jet[END_REF] to solve the low-mach formulation of the Navier-Stokes equation for studying supercritical liquid jets. In this study, the liquid velocity for advecting the interface is constructed using Aslam's extrapolation method in a similar to [START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF]. In [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF], a pressure-based method was developed for the simulation of compressible and turbulent two-phase flows with phase change. This formalism described atomization processes (high Weber and Reynolds numbers), allowing an accurate representation of encapsulated gas structures with their own thermodynamic pressure. A Coupled Level Set/Volume of Fluid (CLSVOF) interface capturing method was implemented where the VOf equation was adapted to consider the compressibility effects and vaporization:

∂α l ∂t + ∇ • (α l u) = α l (1 -α l ) D + α l ∇ • u -α l ṁ 1 ρ l - 1 ρ g + ṁ ρ l (1.14)
whehe α l is the liquid volume fraction, D is a compressibility term and ṁ is the volumetric vaporization rate. The LHS of eq. 1.14 is similar to the VOf equation used in the incompressible formalisms and first three terms of the RHS are dedicated to the representation of the compressibility effects cause by the changes in the pressure and the vaporization. The last term of eq. 1.14 is a volumetric sink term that is in charge of removing the quantity of liquid dictated by the vaporization rate and it is also responsible of the interface regression. Additionally, a pressure equation was developed to account for the compressibility effects and the Stefan flow in the Navier-Stokes equation. Among the simulation results presented by [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF], there is a compressible Homogeneous Isotropic Turbulent (HIT) where the interface regression and the Stefan flow cause by the phase chance are taken into account. These results are significant because they open the door to the study of numerous physical phenomena that are difficult or impossible to analyze by experiment. However, coupling between the energy transfer and the species concentration with the vaporization rate was not yet implemented.

Turbulent mixing

The influence of turbulent flows in the evolution of the mixing between two or more species is of great importance in nature and engineering, e. g., combustion. After the fuel atomization and evaporation in a combustion chamber, the turbulent flow enhances the fuel vapor and air mixing at a molecular level, improving the combustion. Turbulent mixing can be defined as a three-stage process: entrainment, dispersion or stirring, and diffusion of one or more species into another [START_REF] Paul E Dimotakis | Turbulent mixing[END_REF]. In the most simple cases, the mixing is passive, meaning that the mixture of the different species does not affect the configuration's dynamic.

Examples of this kind of mixture are the mixture of two fluids with the same density, the dispersion and mixing of particles with a small Stokes number in airflow, and the dilution of small concentration ink in a liquid. This can be characterized as Level-1 mixing and is mainly driven by the turbulent flow.

In level-2, flows dynamics are affected by the mixing. This kind of mixing can be observed in the mixture of two fluids with different densities under the influence of gravity and in the case of mixing due to a temperature gradient (natural convection). In level-3, the mixing is coupled with the flow dynamic and also changes the properties and composition of each species. A clear example of this type of mixing is combustion, where the flow and composition are dominated by the reaction rate, which directly depends on the level of mixture between the fuel and the oxygen. Understanding the turbulent mixing processes is especially important to the turbulent reacting flows community. This is because mixing among the reactants is enhanced by turbulence and may thus significantly increase the reaction rate ([111]).

Turbulent mixing in single-phase flows

The studies dedicated to turbulent mixing focus mostly on the Level-1 mixing with high Reynolds numbers of a few canonical cases, e.g., grid/isotropic turbulence, channel and pipe flow, jets, etc. [START_REF] Paul E Dimotakis | Turbulent mixing[END_REF]. Generally, there are two different ways to analyze the turbulent mixing of a passive scalar. First, we can focus on the spatial and temporal variation of the integral-scale properties (e. g., variance and flux) to determine the mixing properties of a passive scalar. And the second is to focus on the inertial and dissipation scales.

Experimental measurements of the single-point probability density function (pdf) of a passive scalar in homogenous turbulent shear flow made by [START_REF] Tavoularis | Experiments in nearly homogenous turbulent shear flow with a uniform mean temperature gradient. part 1[END_REF] indicated an almost normal distribution for the scalar pdf. Later, various experimental and theoretical works show that, under certain conditions, the scalar pdf presents an exponential tail due to anomalous mixing ( [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]). [START_REF] Jayesh | Probability distribution, conditional dissipation, and transport of passive temperature fluctuations in grid-generated turbulence[END_REF] studied the evolution of the scalar pdf and dissipation rate for passive temperature fluctuations in decaying grid-generated turbulence; concluding that the scalar pdf has a pronounced exponential tail for turbulence Reynolds number Re λ greater than 70 but below this values the pdfs can be modeled by a normal distribution.

Regarding the analysis of the small-scale behavior of the passive scalar fluctuations. [START_REF] Katepalli R Sreenivasan | The passive scalar spectrum and the obukhovcorrsin constant[END_REF] compiled and analyzed the slope of the scalar spectrum for various shear flows. They observed that the expected slope of -5/3 in the scalar energy spectrum is reached for Re λ > 1000, suggesting the possibility of a "universal" Obukhov-Corrsin constant. Analysis of the third-order statistics (skewness) of the scalar pdf has shown the anisotropy at the dissipation and inertial scales ( [START_REF] Katepalli R Sreenivasan | On local isotropy of passive scalars in turbulent shear flows[END_REF]). This phenomenon was numerically studied by [START_REF] Pumir | A numerical study of the mixing of a passive scalar in three dimensions in the presence of a mean gradient[END_REF]. They used a DNS of passive scalar mixing in the presence of a mean gradient, obtaining a skewness of order 1.

Since experimental measurements of mixing in high Reynolds numbers are difficult to obtain, the reliance on DNS simulations increases to study the mixing. [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF] used DNS for studying the turbulent mixing of a passive scalar in an isotropic turbulent flow. For the resolution of the Navier-Stokes equation, a pseudospectral method was used, and to maintain a statistically stationary velocity field, the simulations were forced. This mean that energy is added to the simulation in a specific range of wavenumber at the same rate that is dissipated. Among the observations of [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF], we can find the effect of the initial conditions on the mixing process and the evolution of the probability density function of the scalar and scalar-dissipation rate during the mixing. Also, the scalar pdf's tendency to become self-similar in the later stages of the simulations. As an extension of the work done by [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] used DNS to study the mixing of two passive scalars in stationary homogeneous, isotropic turbulence. This work was focused on the evolution of the joint pdf of the two passive scalars to improve the models of multi-scalar mixing.

For high Sc numbers, the scalar fluctuations are dissipated at higher wavenumbers than the velocity fluctuations ( [START_REF] Warhaft | Passive scalars in turbulent flows[END_REF]). The effects of the Schmidt number on the turbulent transport of a passive scalar were studied by [START_REF] Yeung | Schmidt number effects on turbulent transport with uniform mean scalar gradient[END_REF], using DNS of forced isotropic and homogeneous turbulence. They observed a modest inertialconvective range consistent with Bachelor's predictions in the scalar spectrum. More recently, [START_REF] Yeung | High-reynolds-number simulation of turbulent mixing[END_REF] studied a 2048 3 DNS of the mixing of passive scalars in turbulent flows, observing a departure from small scale isotropy. [START_REF] Yeung | Spectrum of passive scalars of high molecular diffusivity in turbulent mixing[END_REF][START_REF] Yeung | Direct numerical simulation of turbulent mixing at very low schmidt number with a uniform mean gradient[END_REF] used DNS to study the passive scalar mixing with a very low Schmidt number, and the aim was to confirm the existence of the -17/3 power of the scalar spectral density in the inertial-diffusive range.

Another interesting phenomenon that can affect turbulent mixing is the presence of a shock wave; this is particularly important for supersonic propulsion engines. The scalar mixing and velocity dynamics in the presence of a shock wave have been recently studied by [START_REF] Boukharfane | Evolution of scalar and velocity dynamics in planar shock-turbulence interaction[END_REF] using DNS of compressible singlephase multi-components flows. The results of [START_REF] Boukharfane | Evolution of scalar and velocity dynamics in planar shock-turbulence interaction[END_REF] shows the enhancement of the scalar mixing when a shock wave is present. Moreover, [START_REF] Buttay | Turbulent mixing and molecular transport in highly under-expanded hydrogen jets[END_REF] studied the turbulent mixing and molecular transport in a supersonic hydrogen jet. They also compared the performance of the Hirschfelder and Curtiss approximation for the molecular diffusion (which describes the transport by molecular diffusion as a function of the diffusion coefficient and the vapor mass fraction gradient) with the results obtained with a detailed multi-component transport representation.

Turbulent mixing in two-phase flows

The aforementioned studies of the passive scalar mixture in a turbulent flow use a single-phase configuration. Where the source of the passive scalar is defined at the initialization of the simulation as one or various blobs of vapor arbitrarily distributed throughout the domain. Moreover, the turbulent mixing of a passive scalar has been studied with a Lagrangian point-particle approach. An example is the work of [START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF], where the effect of the spatial distribution of the droplets on the evaporation and turbulent mixing is studied. They use a spectral DNS solver for the evolution of the turbulent carrier phase and a Lagrangian model for the dispersed phase, where the droplets are assumed as spherical with a diameter much smaller than the grid space. Here, the vaporizing droplets act as vapor point sources. Various publications focus in the study of the turbulent mixing in similar scenarios ([126, 127, 128]). The main interest of this kind of work is to understand the effect of the vaporization source term in the scalar advection/diffusion equation on the evolution of the mean vapor mass fraction, fluctuations and dissipation rate.

A limitation of the Lagrangian point-particle approximation is that the sim-

L. Germes Martínez

DNS of phase change June 28, 2022 ulations are restricted to particles smaller than the Kolmogorov scale. To simulate particles larger than the smaller length scales of the turbulent carrier flows, a fully resolved DNS must be performed ( [START_REF] Balachandar | Turbulent dispersed multiphase flow[END_REF]). This kind of configuration was adopted by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. They studied the mixing process in a two-phase HIT configuration with a vaporizing liquid. Statistical and spectral analyses were presented to demonstrate the effects of the liquid volume fraction in the velocity and scalar fields. In this work, the coupling between the vaporization rate and the flow dynamics via the velocity jump condition at the interface is not considered. Consequently, it is limited to the study of low vaporization rates, where the phase change process barely affects the interface and flow dynamics. [START_REF] Bouali | Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows[END_REF] used the same configuration of [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF] for studying the influence of the liquid/gas interface on the scalar dissipation rate. A numerical contribution made by [START_REF] Bouali | Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows[END_REF] is the implementation of Aslam-Chiu's extension for imposing the vapor mass fraction at the interface.

Another approach that has shown very promising results for studying turbulent mixing is the Large Eddy Simulations (LES), augmented by Sub-Grid Scale (SGS) modeling [START_REF] Paul E Dimotakis | Turbulent mixing[END_REF]. An application of this kind of simulation was presented by [START_REF] Wegner | Comparative study of turbulent mixing in jet in cross-flow configurations using les[END_REF], which used LES to study the effect of the angle between the jet and the oncoming cross-flows. Concluding that it is favorable to inject the jet at an angle inclined against the oncoming main flow.

Chap. 2 | Incompressible two-phase flow simulations with phase change

Introduction

In this chapter, and throughout this manuscript, a system composed of two immiscible fluids (e.g., a single-component liquid and a mixture of gases) is studied. First, the governing equations of our incompressible formalism dedicated to the DNS of two-phase flows with phase change are presented. Since the one-fluid formalism is adopted to represent each phase, special attention must be paid to the jump conditions at the interface. For this reason, the next part of the chapter is a discussion on the standard jump conditions used for inert two-phase flow simulations and the additional jump conditions required for phase change. This is followed by the numerical methods implemented in the in-house ARCHER code [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]. This part discusses modifications made to the interface capturing method to account for the evaporated liquid. In addition, three different methods to handle the velocity jump in the projection method are presented. We then show results in a 2D static cylinder configuration to validate and compare the three methods. Finally, application results are presented to show the accuracy and robustness of our incompressible formalism.

Governing equations

In this first section, the Navier-Stokes equations are presented under the assumption of incompressibility. Then, the interface jump conditions necessary to respect the conservation of mass, momentum, and species when applying the one-fluid formulation and to consider the phase change in the system are discussed.

Incompressible two-phase flow simulations with phase change

Navier-Stokes equations

The continuity and Navier-Stokes equations for incompressible flows can be written as follows:

     ∇ • u = 0 ∂u ∂t + (u • ∇) u = - ∇P ρ + ∇ • (2µ ¯ ) ρ + f v ol (2.1a) (2.1b)
where ρ, u and P are the density (kg • m -3 ), velocity (m • s -1 ) and pressure (P a), respectively. µ is the dynamic viscosity

(kg • m -1 • s -1 ) , f v ol is the volume forces (m • s -2
) and ¯ = 1 2 ∇u + ∇u T is the viscous stress tensor. The vaporization rate is estimated as a function of the local temperature and the normal gradient of the vapor mass fraction. As a result, the energy and species conservation equations are added to the system:

ρc p ∂T ∂t + u • ∇T = ∇ • (λ∇T ) (2.2)
where T is the temperature (K), c p is the heat capacity

(J • kg -1 • K -1 ) and λ is the thermal conductivity (W • m -1 • K -1 ).
Given that the liquid phase is considered as monocomponent, the species conservation equation is only described in the gas phase:

∂Y v ∂t + ∇ • (uY v ) = ∇ • (D m ∇Y v ) (2.3)
where Y v is the vapor mass fraction and D m is the mass diffusivity (m 2 • s -1 ).

It is worth mentioning that even though the velocity field is not divergencefree at the interface when the phase change is considered, the assumption of incompressibility still stands in each phase. It has been done in previous work in the literature ( [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Schlottke | Direct numerical simulation of evaporating droplets[END_REF][START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible navierstokes equations with phase change[END_REF]). But this is possible with the condition that the velocity divergence produced by the Stefan flow at the interface must be evacuated of the domain, e. g., with an outflow boundary condition.

Interface jump conditions

Only a single governing equation system is solved for all phases involved when the one-fluid formalism is used. To satisfy the conservation of mass, energy, momentum, and species, jump conditions are needed when there is an interface. In other words, the jump conditions allow the phases to interact. The jump operator ([•] Γ ) across the interface Γ is defined as [ξ] Γ = ξ l -ξ g where ξ a given scalar.

Commonly, for the cases without phase change, the following jump conditions are used:

[ρ] Γ = ρ l -ρ g (2.4) [µ] Γ = µ l -µ g (2.5) [P -n • (2µ ¯ ) • n] Γ = σκ (2.6)
where σ as the the surface tension (N • m -1 ), κ is the curvature (m -1 ) and n is the vector normal to the interface pointing from the liquid to the gas phase.

Mass conservation

The mass jump condition is obtained by performing a mass balance on a control volume (dV ) containing a portion of the interface between the phases. A representation of the velocities defined in the interface, the control volume, and the part of the interface (dΓ) is shown in the figure 2.1. Assuming that the interface is infinitely thin, mass accumulation is not possible. This means that the mass flux to the interface on the liquid side has to be the same as the mass flux in the gas phase. This can be written as [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF][START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF]:

ρ l u Γ l -u Γ • n -ρ g u Γ g -u Γ • n = 0 (2.7)
where u Γ is the interface velocity (m • s -1 ), u Γ l and u Γ g are the velocities (m • s -1 ) at the liquid and gas side at a distance infinitely small to the interface, respectively. Because we assume that the liquid and gas density are constant in time and space in the incompressible formalism, there is not need to specify that ρ l and ρ g has to be defined close to the interface. The first term on the LHS of eq. 2.7 represent the mass flux toward the interface from the liquid phase which represent the local vaporization rate ω. In addition, the second term on the LHS of eq. 2.7 represent the vapor mass flux from the interface to the surrounding air:

ρ l u Γ l -u Γ • n = ρ g u Γ g -u Γ • n = ω (2.8)
For cases without phase change, the equality of the normal component of the velocities on both sides of the interface can be deduced from eq. 2.7: Furthermore, assuming the no-slip condition at the interface, the tangential components of the velocities must also be equal. Therefore, the velocity jump condition is:

u Γ • n = u Γ l • n = u Γ g • n (2.9)
[u] Γ = 0 (2.10)
Then, if mass transfer is considered, another velocity can be defined [START_REF] Calimez | Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs[END_REF]:

s d = u Γ -u Γ l (2.11)
where s d is the interface regression speed (m • s -1 ). s d is collinear to the normal of the interface because the no-slip condition at the interface is preserved, and its definition simplifies to the difference of the normal components of u Γ and u Γ l :

s d = u Γ -u Γ l • n • n (2.12)
from eq. 2.8, we have:

s d = - ω ρ l n (2.13)
Depending on which side of the interface we are on, u Γ can be defined in two different ways. From the liquid side, u Γ is the sum of the liquid velocity at the interface (u Γ l ) and the interface regression speed (s d ):

u Γ = u Γ l + s d (2.14)
And following the same analysis, u Γ can be also be defined from the gas phase L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022 as:

u Γ = u Γ g + ρ l ρ g s d (2.15)
The difference between equations 2.15 and 2.14 is the velocity jump across the interface:

[u] Γ = ρ l -ρ g ρ g s d (2.16)
Based on eq. 2.16 we can see the two physical phenomena responsible for the velocity discontinuities at the interface. First, the density difference between the liquid and gas phase (ρ l -ρ g ). The second is the vaporization rate at the interface, represented in eq. 2.16 by the interface the regression speed. In the incompressible formalism, the temporal variation of the velocity jumps depends only on the interface temperature and the vapor concentration near the interface because the densities remain constant.

The jump condition for momentum conservation is also affected by mass transfer due to the interface regression velocity. Adding the evaporation-related term to eq. 2.6, the pressure jump condition is:

[P -n • (2µ ¯ ) • n] Γ = σκ -ω2 1 ρ Γ (2.

17)

Energy and species conservation

To solve the energy conservation equation in the one-fluid formalism, thermodynamic equilibrium is assumed at the interface, i.e., the temperature of each phase is the same at the interface:

T Γ = T Γ l = T Γ g
This assumption greatly simplifies the estimation of the spatial derivative when solving the diffusion term of eq. 2.2 since

[T ] Γ = 0
However, discontinuities in physical properties, such as thermal conduction and heat capacity, are present:

[λ] Γ = λ l -λ g (2.18) [c p ] Γ = c pl -c pg (2.19)
In addition, the latent heat must be considered for cases with phase change, which is the amount of energy needed by the molecules to pass from one state to another. Consequently, the heat flux across the interface is discontinuous, and an additional jump condition must be considered to respect the energy conservation [START_REF] Calimez | Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs[END_REF]:

λ l ∇T Γ l • n -λ g ∇T Γ g • n = -h lv ω (2.20)
where h lv is the latent heat (J • kg -1 ). In the case of a liquid evaporating in a gas at high temperature, the first term on the LHS of eq. 2.20 represents the heat flux from the interface to the liquid and the second is the heat flux from the gas in the direction to the interface. It is important to mention that eq. 2.20 is derived under the assumption that the energy source term due to friction and the viscous dissipation are negligible.

Finally, to respect the conservation of species, for a single component liquid, the interface jump condition for Y v is

ρ g D m ∇Y v • n| Γ g = ω(1 -Y v s ) (2.21)
where Y v s is the vapor mass fraction at the interface. To estimate Y v s , thermodynamic equilibrium is assumed at the interface, allowing a direct relation between the energy and Y v . To this end, the Clausius-Clapeyron relation is used:

X = exp - h lv M v ap R 1 T Γ - 1 T B (2.22)
where X is the vapor mole fraction at the interface, M v ap is the molar mass of the vapor (kg • mol -1 ), R is the ideal gas constant

(J • K -1 • mol -1
), T Γ is the temperature (K) of the interface and T B is the boiling temperature (K). From the mass definition, Y v s can be calculated as:

Y v s = XM v ap XM v ap + (1 -X)M g (2.23)
where M g is the molar mass of the inert gas. And, the denominator of the RHS of eq. 2.23 represents the mean molar mass of the gas mixture.

Interface jump conditions of the different variables reveal the strong coupling between heat/mass transfer and the velocity jump. Handling this coupling is challenging and requires robust and accurate numerical methods dedicated to DNS of two-phase flow with phase change.

Numerical Methods

This section discusses numerical strategies used to perform DNS of multiphase flows with phase change. Numerous difficulties appear when dealing with these simulations due to the diversity of physical processes simultaneously. In addition, special care is required when implementing the jump conditions at the interface described in the previous section. We will start with the description of the interface capturing method, adapted to include the interface regression due to phase change. Then, we describe the projection method for the time evolution of the Navier-Stokes equations. And finally, we present the numerical methods implemented to follow the temporal evolution of the temperature and species fields.

Interface capturing

The Coupled Level-set/Volume of Fluid (CLSVOF) method implemented in the ARCHER code [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF] is extended to handle multiphase flows with phase change. To this end, a critical step is modifying the VOF equation since it controls the conservation of mass in our CLSVOF method.

For demonstrating the VOF equation, the first step is to define the density as a weighted average of the density of each phase:

ρ = α l ρ l + α g ρ g
where α l and α g are the local liquid and gas volume fraction, respectively. Then, this definition of ρ is introduced into the continuity equation (eq. 2.1a). To account for the variation of the liquid volume in the presence of the phase change, a sink term is added:

∂α l ∂t + ∇ • (α l u l ) = - ṁ ρ l (2.24)
where ṁ represent the amount of evaporated liquid volume in each cell containing an interface. ṁ is estimated by:

ṁ = ρ l Σ s d (2.25)
where Σ = S Γ V cell is the surface density (m -1 ) in a mixed cell, S Γ is the surface obtained in the PLIC reconstruction and V cell is the cell volume.

The convective term of eq. 2.24 is discretized using a modified version of the "coupled" second-order conservative operator split advection scheme defined in [START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF]. For the interface convection, a divergence-free liquid velocity field (u l ) is used. This velocity field is continuous at the interface, and its construction is explained in the following section. The first step is to advance the liquid volume fraction in each direction without considering the evaporation term:

                   αl -α n l ∆t + (α n l u l ) i+ 1 2 -(α n l u l ) i-1 2 ∆x = αl u l,i+ 1 2 -u l,i-1 2 ∆x αl -αl ∆t + ( αl v l ) j+ 1 2 -( αv l ) j-1 2 ∆y = αl v l,j+ 1 2 -v l,j-1 2 ∆y ᾱl -αl ∆t + ( αl w l ) k+ 1 2 -( αw l ) k-1 2 ∆z = ᾱl w l,k+ 1 2 -w l,k-1 2 ∆z (2.26a) (2.26b) (2.26c)
The convection in each direction is coupled in a final equation as a second step. Here, the sink term is added to obtain the following equation:

α n+1 l -ᾱl ∆t = -ᾱl w l,k+ 1 2 -w l,k-1 2 ∆z -αl v l,j+ 1 2 -v l,j-1 2 ∆y -αl u l,i+ 1 2 -u l,i-1 2 ∆x - ṁ ρ l (2.27)
Since liquid evaporation is explicitly accounted for by adding a sink term on the RHS of eq. 2.27; local undershoot of the liquid volume fraction (α l < 0) may occur. This undershoot may be caused by a discrepancy between the available liquid in a mixed cell and the magnitude of the sink term. To respect mass conservation and ensure accurate interface regression due to evaporation, the remaining fraction of evaporated liquid is distributed to the neighboring mesh cell located in the direction opposite to the interface normal computed during the PLIC reconstruction.

In a VOF/LS coupling, it is necessary to solve both of their transport equations simultaneously, the Level Set transport equation is defined as:

∂φ ∂t + ∇ • (φu l ) = 0 (2.28)
where φ is the signed distance function between a grid point and the interface, it takes positive values for the nodes in the liquid phase and negative for those in the gas. The discretization of the convective term in the level set equation is made following the original split advection scheme defined in [START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF][START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]. The signed distance function is used to obtain the geometrical information of the interface, e.g., the vector normal to the interface (n) and the curvature (κ):

n = ∇φ |∇φ| κ = -∇ • n (2.29) L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022
A redistancing algorithm is applied following the same principle of previous work ( [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF], [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF], [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]).

Projection method

As discussed in the previous chapter, the treatment of interfacial velocity discontinuities is one of the significant challenges when simulating two-phase flows with phase change. We also mentioned some of the strategies found in the literature adapted to this kind of velocity field. This work will compare two of these methods: the implicit method and the explicit method. The main reason for this comparison is to choose the method that best suits the main objective of this work, which is to develop a compressible formalism to simulate atomization regimes with phase change.

A common point of both methods is that the velocity field is calculated by solving the momentum equation with a projection method. Moreover, this equation is discretized on a staggered Eulerian mesh. The scalar variables (P , T , α l , etc.) are defined at the center of the cell, and the velocities are located at the faces of the cell.

Implicit method

The first step of the projection method is to calculate an intermediate velocity field (u * ). To this end, eq. 2.1b is temporally discretized as follows:

u n+1 -u * ∆t + u * -u n ∆t + (u n • ∇) u n = ∇ • (2µ ¯ ) ρ n+1 + f v ol I - ∇P ρ n+1 (2.30)
Then, we split eq. 2.30 and take the terms defined by I to calculate u * :

u * = u n -∆t (u n • ∇) u n - ∇ • (2µ ¯ ) ρ n+1 + f v ol (2.31)
In the eq. 2.31, the convective terms are calculated with a 5th order WENO scheme, and the viscous terms are solved with the method proposed by [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. It is essential to mention that the surface tension force is not considered in calculating the intermediate velocity. What is left of eq. 2.30 is used for the computation of the pressure field:

u n+1 = ∆t ∇P ρ n+1 + u * (2.32)
The next step is to apply the divergence operator in each side of eq. 2.32:

∇ • u n+1 = ∇ • ∆t ∇P ρ n+1 + ∇ • u * (2.33)
Generally, the velocity field in the next time step is divergence-free when a incompressible formalism is implemented, i.e., ∇ • u n+1 = 0. But, when phase change is considered, there is a velocity divergence at the interface, which magnitude can be expressed using the velocity jump condition (eq. 2.16):

∇ • u n+1 = ρ l -ρ g ρ g s d δ Γ (2.34)
where δ Γ is a Dirac function; used to identify the cells of the mesh that contain an interface.

Then, eq. 2.34 is injected in eq. 2.33 to obtain the Poisson equation for the pressure:

-∇ • ∆t ∇P ρ n+1 = ∇ • u * - ρ l -ρ g ρ g s d δ Γ (2.35)
Here, the Poisson solver used for this Helmholtz type equation (eq. 2.35) consists of an MGCG (MultiGrid preconditioned Conjugate Gradient) method coupled with the Gauss-Seidel "Red-Black" iterative scheme.

Finally, the velocity field in the next time step is computed with eq. 2.32.

Dirat function

There are several approaches in the literature to describe the Dirac function.

Here, we will use two of them:

• Continuum approach:

δ Γ = |∇α l | (2.36)
A 2nd order central scheme is used for the discretization of eq. 2.36. This method smears the velocity jump over a few cells around the interface. Some authors have observed the formation of spurious currents in the onefluid velocity when implementing the velocity jump using this approach [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Soo | Direct numerical simulation of incompressible multiphase flow with phase change[END_REF]. This can further complicate the construction of a liquid velocity field for interface advection. • Sharp approach:

δ Γ = S Γ V cell (2.37)
This approach is the same used in this work to estimate the volumetric

L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022
evaporation rate in the sink term of the VOF equation. The local surface S Γ is given by the PLIC reconstruction. It allows a more accurate treatment of the velocity jump, the affected cells being only those containing an interface. In addition, it considers the amount of surface area in the grid cell, which could improve the estimation of Stefan flux.

Liquid velocity construction

When an implicit treatment of the velocity discontinuity is implemented, the Stefan flow is automatically accounted for in the one-fluid formulation of the velocity field. This means that all three velocities (interface, liquid, and gas) at the interface defined in the section 2.2.2 are present in the same mesh cell. Consequently, the velocity field generated from the projection method is not suitable for interface advection. Several techniques can construct a velocity field for liquid phase convection in the literature. Most of them are based on the idea of extracting the velocity of the liquid from the one-fluid velocity [START_REF] Palmore | A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows[END_REF][START_REF] Lc Malan | A geometric vof method for interface resolved phase change and conservative thermal energy advection[END_REF][START_REF] Scapin | A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows[END_REF].

The strategies used in this method consist of extending the velocity in the liquid phase into the gas phase by eliminating the velocity jump induced by the phase change. Therefore, a continuous velocity field that is divergence-free by construction is obtained. The first step of this approach is to construct the Stefan flow velocity (u S ) using a velocity potential:

       ∇ • ∇ϕ ρ n+1 = ρ l -ρ g ρ g s d δ Γ u S = ϕ ρ n+1 (2.38a) (2.38b)
Equation 2.38a is solved using the same numerical schemes used for eq. 2.35. Next step is to subtract the Stefan velocity from the one-fluid velocity to obtain the liquid velocity u l :

u l = u -u S (2.39)
u l is then used for the advection of the VOF and the LS functions.

Explicit method

For the explicit method, we will use the method proposed by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] dedicated to handling discontinuous velocity fields at the interface. This method follows two separate velocity fields, one for each phase. Both velocity fields are coupled in the Poisson equation for the pressure and the extension of the gas velocity field in the liquid phase. Also, the velocity extensions are calculated following the strategies presented in equations 1.12, which allow having continuous liquid/gas velocity fields at the interface. A Ghost fluid method (GFM) populates the ghost cells on both interface sides.

The first step of the projection method is the same as for the implicit method. The temporal discretization of the momentum equation described in eq. 2.30 is used for the computation of the intermediate velocity (u * ). This step is similar for both liquid and gas velocity fields:

         u * l = u n l -∆t (u n l • ∇) u n l - ∇ • (2µ ¯ l ) ρ n+1 + f l,v ol u * g = u n g -∆t u n g • ∇ u n g - ∇ • (2µ ¯ g ) ρ n+1 + f g,v ol (2.40a) (2.40b)
For the equations 2.40, the convective terms are calculated with a 5th order WENO scheme and the viscous terms are solved with the method proposed by [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. Once again, the remainder of eq. 2.30 is used to calculate the pressure and velocity fields. For the pressure, we use eq. 2.33 and assume at first that the velocity field at the next time step is divergence-free (∇ • u n+1 = 0):

-∇ • ∆t ∇P ρ n+1 = f (2.41)
where

f =    ∇ • u * l if φ > 0, ∇ • u * g if φ ≤ 0,
This is a crucial part of the method because it is the point of coupling between the phases during the solution of the Poisson equation. As for eq. 2.35, the Poisson solver used consists of an MGCG (MultiGrid preconditioned Conjugate Gradient) method coupled with the Gauss-Seidel "Red-Black" iterative scheme. Then, the real part of the velocity field of each phase is obtained with:

   u n+1 l = u * l -∆t ∇P n+1 ρ n+1 if φ > 0, u n+1 g = u * g -∆t ∇P n+1 ρ n+1 if φ ≤ 0, (2.42)
Now, to use the GFM, the liquid and gas velocity extensions must be estimated. For the liquid, a ghost pressure (P ghost ) is computed following:

-∇ • ∆t ∇P ghost ρ n+1 = ∇ • u * l (2.43) L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022
Eq. 2.43 is solved following the same methods used for eq. 2.35. And, the liquid velocity extension (u ghost l ) is computed from P ghost and u * l :

u ghost l = u * l -∆t ∇P ghost ρ n+1 (2.44)
The velocity field resulting from the combination of the liquid velocity field and its extension in the gas phase is not divergence-free. To enforce this condition, [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] proposed to extend this velocity field onto its divergence-free part using an intermediate velocity field (W ):

   W = u n+1 l if φ > 0, W = u ghost l if φ ≤ 0, (2.45)
This extra step, compared to the method proposed by [START_REF] Duc Q Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF], improves mass conservation by reducing the deformation of the interface due to numerical instabilities during its convection ( [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]). Then, the following Poisson equation is solved for a potential function (ϕ):

∇ 2 • ϕ = ∇ • W (2.46)
And u n+1 l is obtain by:

u n+1 l = W -∇ϕ (2.47)
The velocity extension for the gas is computed by adding the interface velocity jump to the liquid velocity:

u ghost g = u n+1 l -s d ρ l -ρ g ρ g (2.48)
Finally, the gas velocity field is obtained:

u n+1 g =    u ghost g if φ > 0, u n+1 g if φ < 0, (2.49)
The gas velocity and its extension are continuous and not divergence-free at the interface. In this case, the resulting gas velocity field remains unchanged for the following parts of the formalism. This is because the gas velocity plays no role in the interface convection.

With this formalism , two continuous velocity fields have been obtained in

   u n+1 = u n+1 l if φ > 0, u n+1 = u n+1 g if φ ≤ 0, (2.50)

Energy transport

For the temporal evolution of the temperature field, the energy equation (eq. 2.

2) is solved. The convective term is discretized using a 5th order WENO scheme. Large temperature gradients can occur at the interface, which is a source of error and can lead to numerical heating across the interface. To avoid this, the temperature field of each phase is extrapolated using Aslam's extension method [START_REF] Tariq D Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF]. The idea is to extend the scalar of one phase linearly, following the normal direction of the interface, into the other, populating the ghost grid points. Then a 5th order WENO scheme is used for convection of the liquid temperature and its ghost values (u • ∇T l ) and the same procedure is applied for the gas side (u • ∇T g ). Finally, these terms are introduced into eq. 2.2 depending on the sign of the Level Set function.

A 2nd order central difference scheme is used for the discretization of the diffusion term. Since mass transfer is considered, the energy jump at the interface (eq. 2.20) must be taken into account in eq. 2.2. Consequently, the discretization scheme for the diffusion term is modified at the points near the interface using a GFM. As an example, the discretization of the diffusion term for the 1D case when φ i > 0 and φ i+1 < 0 illustrated in fig. 2.2 is presented. First, eq. 2.20 is discretized with the help of

θ i = |φ| i |φ| i + |φ| i+1
where θ i represent the distance of the point i from the interface with respect to the point in the other side of the interface:

λ l T Γ -T i θ∆x -λ g T i+1 -T Γ (1 -θ)∆x = -ωh lv (2.51)
Then, from eq. 2.51 the interface temperature (T Γ ) is obtained:

T Γ = λ g θT i+1 + λ l (1 -θ)T i -ωh lv θ(1 -θ)∆x λ l (1 -θ) + λ g θ (2.52)
We only need to determine the temperature derivative in the liquid phase for this case. For this reason, we solve for the first term of the LHS of eq. 2.51 and replace the interface temperature for its value given by eq. 2.52. After some manipulation and the definition of a ghost thermal conductivity coefficient (λ GH ) which depend on the level set function:

λ GH = λ l λ g λ l (1 -θ) + λ g θ (2.53)
An equation for the temperature derivative in the liquid side is obtained:

λ l T Γ -T i θ∆x = λ GH T i+1 -T i ∆x - ωh lv (1 -θ) λ g (2.54)
An advantage of this method is that the energy jump condition is added to the system without further modifying the numerical scheme. In addition, the latent heat is weighted at each point as a function of its distance from the interface. For the case studied in fig. 2.2, the discretized diffusion term in the energy equation has the following form:

∂ ∂x λ ∂T ∂x i = λ GH T i+1 -T i ∆x -ωh lv (1-θ) λg -λ l T i -T i-1 ∆x ∆x (2.55)
There are three other cases depending on the phase of the point considered and the phase of the point across the interface. The discretized diffusion terms Incompressible two-phase flow simulations with phase change for each case are summarized in the table 2.1

Cases ∂ ∂x λ ∂T ∂x i φ i > 0 and φ i+1 < 0 λ GH T i+1 -T i ∆x - ωh lv (1-θ) λg -λ l T i -T i-1 ∆x ∆x φ i < 0 and φ i+1 > 0 λ GH T i+1 -T i ∆x + ωh lv (1-θ) λ l -λg T i -T i-1 ∆x ∆x φ i > 0 and φ i-1 < 0 λ l T i+1 -T i ∆x -λ GH T i -T i-1 ∆x - ωh lv (1-θ) λg ∆x φ i < 0 and φ i-1 > 0 λg T i+1 -T i ∆x -λ GH T i -T i-1 ∆x + ωh lv (1-θ) λ l ∆x
Table 2.1: Discretized energy diffusion terms for the different cases.

Vapor mass fraction transport

To transport the vapor mass fraction (Y v ), eq. 2.3 is solved only in the gas phase as for eq. 2.2, the diffusion term is discretized using a 2nd order central difference scheme, and the convective term is computed with a 5th order WENO scheme.

To enforce the Dirichlet boundary condition on the interface (Y v s ), a 1st order Aslam-Chiu extension method is implemented. In this method, the main idea is to introduce a mirror point in the gas phase, which is the image of the ghost point of Y v in the liquid, initially estimated with the Aslam extension method. Then, the mirror point is used to correct the estimation of the Aslam extension, following the value of the boundary condition given by eq. 2. [START_REF] David | An interface tracking method for a 3d eulerian hydrodynamics code[END_REF]. More details on this method are available in [START_REF] Bouali | Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows[END_REF]. The estimation of the ghost point of the vapor fraction (Y Gh v ) in the liquid is computed as:

Y Gh v = 2Y v s -Y ext v -2φ d Y v dn (2.56)
where

Y ext v
is the first estimation of the ghost point, d Yv dn is extracted from the linear Aslam's extension done before the convection step of Y v and Y v s is obtained with eq. 2.21.

The interface temperature (T Γ ) needed for eq. 2.22 is estimated with (in 1D):

T Γ = T l,i |φ i+1 | + T ghost l,i+1 |φ i | |φ i+1 | + |φ i | (2.57)
where T ghost l is the extended liquid temperature on the gas phase. Eq. 2.57 can be generalized into two and three dimensions. If the liquid and the gas temperatures are used instead of the liquid temperature and its extension to interpolate L. Germes Martínez For both scalar fields (T and Y v ), a correction step is implemented when a grid node passes from being in one phase to the other due to the evolution of the interface. During this step, Y v and T of a 'switching' point i are replaced with its ghost values given by Aslam's extension method during the previous time step. After integrating the interface capturing method, the correction aims to avoid artificial heating/cooling in regions near the interface.

Summary of methods

This work implements three different methods to deal with the various numerical challenges that arise when considering phase change. The principles and mathematical basis were explained in the previous sections. Here, we want to recapitulate each method and highlight their differences.

Numerical challenge Method A Method B Method C [u] Γ Implicit Implicit Explicit δ Γ |∇α l | S Γ V cell - u l
Stefan flow cancellation Stefan flow cancellation 2 velocities + GFM T and Y v convection u u u l , u g Table 2.2: Summary of the three methods main characteristics.

For method A, the velocity jump at the interface is treated implicitly with the help of a source term in the Poisson equation and a Dirac function estimated with the gradient of the liquid volume fraction. At the same time, the Stefan flow is removed from the one-fluid velocity to construct the liquid velocity. In addition, the convection of temperature and vapor mass fraction is performed with the one-fluid velocity. Method B differs from method A only in calculating the Dirac function, which is done with the surface area calculated during PLIC reconstruction. In method C, the velocity jump is treated explicitly by adding it directly to the gas velocity extension. In addition, the liquid velocity field is constructed during the resolution of the projection method. Therefore, no additional velocity field needs to be considered. Finally, for the convection of scalars in method C, either the liquid or the gas field is used depending on the sign of the level set function.

Simulation results

This section presents three academic cases to validate our incompressible formalism. First, a 2D static cylinder configuration with an imposed evaporation rate is studied to demonstrate the accuracy of our sink term in the VOF equation, combined with the CLSVOF method for capturing reacting interfaces. This configuration is performed using the numerical methods described in the table 2.2. Then, a D 2 law configuration is investigated to validate the implementation of the thermodynamics of the formalism and the coupling between the temperature, the vapor mass fraction, and the evaporation rate. Finally, a convected evaporation cylinder is studied to show the ability of the method to handle a moving interface with an inhomogeneous velocity jump.

The physical properties for the air, water and n-decane are listed in table 2.3. The mass diffusion coefficient for the vapour in air is D m = 2.1 × 10 -5 m 2 .s -1 and D m = 1 × 10 -5 m 2 .s -1 for the n-decane vapour in air. The surface tension for water/air interface is σ = 0.07 N.m -1 and σ = 0.0135 N.m -1 for the n-decane/air interface. This section has two main objectives:

Fluid ρ kg m 3 µ kg ms λ W mK C p J kgK M kg mol h lv J kg T B (K) Air 1.
• Validate the estimation of the sink term in the VOF equation (eq. 2.24) responsible for the interface regression due to phase change. • To test the efficiency of the methods for treating the liquid velocity extension mentioned in section 2.3.2. This allows us to choose the most appropriate method for this work.

To this end, a 2D static water cylinder with a constant evaporation rate is analyzed. Since the evaporation rate is imposed, the energy and species conservation equations are not considered for this configuration. However, the Navier-Stokes equations are solved.

A difficulty in simulating two-phase flows with evaporation is the spurious currents created in the velocity field of the liquid phase closed to the interface.

A probable cause of the formation of these spurious currents in the numerical treatment of velocity discontinuities at the interface and surface forces ( [START_REF] Popinet | A front-tracking algorithm for accurate representation of surface tension[END_REF]). If the magnitude of the spurious velocities is significant, unrealistic interface deformation will occur, leading to problems in mass conservation. The velocity of the liquid phase is assumed to be zero for the static configuration of the cylinder (u l = 0), with the regression velocity (s d ) being the only velocity to influence the interface evolution. On the other hand, the one-fluid velocity must contain the Stefan flow around the interface.

Initially, the velocity of both phases is zero, and outflow boundary conditions are used in all directions. The evolution of the mass of the liquid is relatively easy to predict since a constant evaporation rate and a circular shape are assumed and is computed as follows:

m th (t) = ρ l π (R D -s d t) 2 (2.58)
Here, we choose to investigate two different vaporization regimes. The first one corresponds to a small Stefan flow, and the second represents a large Stefan flow.

Small Stefan Flow

To represent a small Stefan flow, the configuration presented in [START_REF] Malan | Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change[END_REF] is reproduced. Here, the initial cylinder of radius R D = 0.23 m is at the center of a square domain of length l x = l y = 1 m. The interface regression velocity is held at a constant value of s d = 0.05 m • s -1 and the density ratio is ρ l ρg = 2. Simulations are carried out for all 3 methods with a mesh of 128 2 and a constant time step of dt = 0.0005 s. The velocity jump for this configuration is:

s d ρ l ρ g -1 = 0.05 m • s -1 (2.59)
The magnitude of the velocity jump can be considered the maximum value of the Stefan flow closed to the interface. for our three methods. We observed an excellent agreement between theoretical and simulation results for all methods. These results reveal that the circular shape is maintained throughout the life of the cylinder. In addition, it shows that the spurious currents in the liquid velocity field are maintained at a level that allows good interface regression.

Large Stefan Flow

To test our three methods with a large Stefan Flow configuration, we reproduce the configuration presented in [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. Here, for each simulation, the initial cylinder radius is R D = 200 µm and the domain is a square with a side length of l x = l y = 6R D . Simulations are made for five mesh sizes: 32 2 , 48 2 , 64 2 , 96 2 and 128 2 and a CFL number of 0.2. The constant interface regression speed is s d = 10 -3 m.s -1 . Also, the density ratio is ρ l ρg = 815.7. This configuration represents a Stefan flow 16.3 times stronger compared to the previous configuration.

The time evolution of the normalized liquid mass for the simulation with a mesh of 128 2 is shown in fig. 2.4. Here again, good agreement is obtained between the theoretical and simulation results for the three methods. Given that the points plotted in figure 2.4 are chosen at every same number of iterations, and the CFL number is constant. A smaller time step is observed in method A, which indicates that the velocity jump is larger than in the other methods.

Then, a mesh convergence study is carried out to analyze the liquid mass error further. For this end, the L p norms (defined with the absolute error of the normalized liquid mass) are plotted in fig. 2.5, which are computed as:

L ∞ = max |m * th (t) -m * (t)| (2.60) L 1 = |m * th (t) -m * (t)| N (2.61) L 2 = |m * th (t) -m * (t)| 2 N (2.62)
Order of convergence between 1 and 2 is obtained for the three methods. These results show an improvement in the magnitude of the errors compared to the results obtained by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] although the order of convergence is almost the same. This difference can be explained by using a conservative method to capture the interface (CLSVOF method). In addition, the magnitude of the errors for the least resolved configurations is smaller for method C, followed by method B. For example, in fig. 2.5a, a radical increase of L ∞ is observed for methods A and B in the simulations with a mesh of 32 2 . Here, in the last moments of the simulation, the spurious current pushes the cylinder away from the center of the domain, causing this increase of the maximum error. Hence, the spurious currents are more prominent in the less resolved simulations for the A and B methods. In addition, with method C, more consistent results are obtained for each mesh.

Results obtained with method C are explained by using velocities without discontinuities for the convective term of the Navier-Stokes and the interface transport equations. In addition, the improvement in method B compared with method C is explained by the sharp definition of the Dirac function. Furthermore, fig. 2.5 shows that all three methods converge to a similar level of error for the simulation with the more refined grid.

The one-fluid velocity field for the three methods for the mesh of 128 2 is shown in fig. 2.6. The Stefan velocity, normal to the interface, can be seen in all three figures. Since we are in an incompressible formalism and the velocity of the liquid is assumed to be null, the maximal magnitude of the Stefan flow can be easily estimated by:

Sd 2 ρ l ρ g -1 = 0.81 m • s -1 (2.63)
This maximum value is well reproduced in each method, with the most significant deviation found in method A. In the color field of fig. 2.6, we can see that the one-fluid velocities constructed with the implicit velocity jump treatment (method A and B) are smooth at the interface compared to method C, where the velocity jump is sharper. The one-fluid velocity for method C is constructed from the liquid and gas velocity following the sign of the level set function.

In Figures 2.7, the liquid velocity field for the three methods is presented. Since the velocity of the cylinder is supposed to be zero in the static configuration, the velocity field shown in figures 2.7 is the spurious velocities. These velocities can influence the transport of the interface, deforming it and causing mass conservation errors. The smallest magnitudes of the spurious velocities are found in method C. Also, the small spurious velocities could be because the liquid velocity field for method C is extended into the gas by a velocity field u ghost l constructed using only the information from the liquid phase. In addition, for methods B and C, spurious currents are formed away from the interface because u S is defined over the entire domain. These spurious currents do not interfere with the advection of the interface since the velocities used are around the mixed cells.

Finally, the gas velocity field and its extension in the liquid phase are shown in fig. 2.8. Here, we can see that the gas velocity is continuous, and the velocity jump is added to the ghost points in the liquid phase following the normal direction of the interface.

Throughout the rest of this manuscript, method C is chosen for our simulations of two-phase flows with phase change. This decision was made following the results obtained for the static vaporizing cylinder configuration with large Stefan flows. Moreover, the nature of the construction of the liquid velocity field in the method A and B could lead to further complications when extending the technique to a compressible formalism, where the source of divergence in the velocity field is not only due to the phase change, e.g., when thermal expansion and spatial and temporal density variations are considered. Since the construction of the liquid velocity field in these methods is based on the cancellation of the velocity divergence at the interface.

D 2 law
A 3D static droplet configuration is compared with the well-known D 2 law to validate the implementation of thermodynamics in the incompressible formalism.

According to this law, the square of the droplet diameter decreases linearly with time due to heat and mass diffusion in the gas [START_REF] William | Fluid dynamics and transport of droplets and sprays[END_REF]. To obtain the normalized squared radius equation, one must solve the continuity equation in spherical coordinates assuming incompressibility and a constant liquid mass flow from the interface to the air ( ṁ = 4πr 2 ρu).

Among the assumptions necessary for the D 2 law, the isolated droplet must be in an infinite continuous medium. Reproducing this condition in a numerical domain is too computationally expensive. It would imply that the droplet interface is far enough away from the domain boundary not to perturb the spatial profile of the velocity, vapor mass fraction, or temperature. To be representative of the finite-domain simulation, we integrate the species and energy conservation equations using the values of each variable at the boundary of the numerical

L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022
domain, yielding:

Y v (r ) = 1 + (Y v BC -1) exp - ṁ 4πρD m 1 r - 1 r BC (2.64) ln T (r ) -h lv cp -ϕ ṁ T BC -h lv cp -ϕ ṁ = - ṁ 4πρD t g 1 r - 1 r BC (2.65)
where Y v BC , T BC are the vapor mass fraction and temperature at the boundary of the domain, respectively. r BC is the distance from the center of the droplet to the limit of the domain. D m is the mass diffusivity. D t g = λ ρgcpg is the thermal diffusivity. ṁ is computed by evaluating eq. 2.64 at the interface

(Y v (R D ) = Y v s ): ṁ = -4πρ r r BC r BC -r ln Y v s -1 Y v ∞ -1 (2.66)
And ϕ can be interpreted as the heat flow entering the liquid. ϕ is obtained by evaluating eq. 2.65 at the interface (T (R D ) = T Γ ):

ϕ = - ṁ h lv + exp -ṁ 4πρD t g 1 r -1 r BC (T BC c p -h lv ) -T Γ c p 1 -exp -ṁ 4πρD t g 1 r -1 r BC (2.67)
This equation considers the length of the numerical domain in our theoretical calculation, which makes the results more comparable.

This configuration is complicated because we are trying to reproduce the theoretical results given by solving the conservation equations in spherical coordinates using a Cartesian grid. To compare the vapor mass fraction and temperature profiles calculated by eq. 2.64 and 2.65 with the simulation results, it is necessary to wait long enough for the configuration to transition from initial conditions to a steady-state. In other words, the vapor mass fraction and the temporal derivative of the temperature must be small enough to be considered negligible. To reduce the computational cost and to avoid droplet motion from the center of the domain, u l = 0. s d is the only velocity responsible for the interface evolution.

A water droplet with an initial radius of R D = 150µm is considered. The initial temperature of the drop is T l = 353K, and the initial temperature of the gas is T g = 573K. The domain is a cube with a side length of 8R D . The simulations are performed with three mesh sizes: 32 3 , 48 3 , and 64 3 . And outflow boundary conditions are used in all directions. On Fig. 2.9, the temperature (2.9a) and the vapor mass fraction (2.9b) are presented at t = 0.08 s for the finest mesh (64 3 ). A spherical symmetric property around the droplet is observed for both fields. Due to the relatively large difference between the liquid and the gas temperature, the temperature gradients in the liquid phase are not noticeable in fig. 2.9a. The temporal evolution of the droplet's normalized diameter is shown in fig. 2.10. The interface regression speed is higher in the early stages of the simulations. The large vapor mass fraction gradients can explain this behavior due to the assumption of dry air at the initial time. In addition, the initial temperature of the liquid is higher than the temperature at the equilibrium state, which means that the evaporation rate will decrease with the temperature temporal evolution. As the temperature and vapor mass fraction stabilize, we can observe that the evolution of the normalized squared diameter becomes linear as predicted by the D 2 law. The same trend appears in the three simulations with different meshes. Here, the difference between the curves is related to capturing the temperature and vapor mass fraction boundary layer. These results show the ability of our method to handle strong temperature gradients at the interface and the coupling between vaporization rate and fluid dynamics efficiently. Fig. 2.11 and fig. 2.12 show the temperature and vapor mass fraction profile in the gas phase obtained by solving Eq. 2.65 and 2.64 and the profiles obtained with three simulations (32 3 , 48 3 and 64 3 ) as a function of the normalized radius r * = r R D . Good agreement with the theory is found in all simulations.

Convected evaporating cylinder

In this section, a simulation of a convected cylinder evaporating is presented. This setup illustrates the ability of the formalism to handle the convective term of the Navier-Stokes equation in the presence of discontinuities in the velocity field combined with the CLSVOF interface capturing method. The initial radius of the cylinder is R D = 300 µm and the initial velocity and temperature of the liquid and gas are

u l = 1 m • s -1 , u g = 0 m • s -1 , T l = 323 K and T g = 873 K, respectively.
Outflow boundary conditions are used in all directions. For temperature and vapor mass fraction, Dirichlet boundary conditions are imposed: T = 873K and Y v = 0. A 2D rectangular domain with length in the x-direction l x = 8R D and in the y-direction l y = 4l x is considered. The simulations are carried out with a grid size of 128 × 512. Fig. 2.13 shows the temperature and vapor mass fraction fields at t = 6 ms. The numerical methods implemented in this work to treat different jump conditions and temperature gradients at the interface allow us to capture various phenomena. For example, we can identify the thermal and inhomogeneous species boundary layers that form at the interface. The presence of strong temperature gradients and vapor mass fraction gradients at the front implies the formation of a recirculation zone due to the moderate Reynolds number investigated (Re D = ρgDu init µg = 40). This recirculation zone causes temperature and vapor mass fraction homogenization behind the cylinder, forming a plume in both fields. In addition, the variation of the vapor mass fraction around the cylinder and of the interface temperature induces an inhomogeneous evaporation rate (see fig.

2.14). As expected, the highest evaporation rates is at the front of the cylinder. This result highlights the importance of considering the local interface temperature and the local surface density to compute the evaporation rate.

In fig. 2.15 the streamlines of the gas flow behind the cylinder colored with the temperature and the isocontour of vu l = 0 in the gas phase (black line) for two simulations are shown: on the right, the evaporating falling cylinder; on the left, we have the same configuration, but this time, without evaporation (heating only). Here, u l is the volume averaged velocity in the liquid phase. A nearly axisymmetric vortex and a slight detachment and enlargement of the recirculation zone are observed in the back of the evaporating cylinder. These differences are explained by the presence of the Stefan flow, which creates an envelope around the interface and separates the airflow from the cylinder. have been compared with those found in [START_REF] Dennis | Numerical solutions for steady flow past a circular cylinder at reynolds numbers up to 100[END_REF] for a solid cylinder. In fig. 2.16 the values of v -u l over a vertical line in the center of the domain (solid line), and the location of the interface (dash-dotted lines) for the water cylinder are shown. These kind of curves allow us to compute the recirculation zone length, since it is the distance from the location of the interface at the back of the cylinder to the place where the velocity changes sign. The case without evaporation is represented in fig. 2.16a and the case with evaporation is shown in fig. 2.16b. Here, we observe a small relative enlargement of the recirculation zone in the case with evaporation. Additionally, good agreement is found for the non-evaporating water cylinder with a 1.4% difference with the reference, while the evaporating configuration has a difference of 4.7%.

To further demonstrate the influence of Stefan flow on the recirculation zone characteristics, two additional simulations are performed for a n-decane cylinder with the same Reynolds number. Also, the same temperature difference between the initial liquid temperature and its boiling temperature (∆T = T B -T l ) is used in each simulation. However, a larger Stefan flow is obtained for n-decane due to its thermophysical properties. Fig. 2.17 shows the streamlines behind the n-decane cylinder: on the right, the evaporating configuration, and on the left, the non-evaporating one. Here we can see that the detachment of the recirculation zone is more pronounced. Additionally, fig. 2.18 shows the values of vu l over a vertical line in the center of the domain (solid line), and the location of the interface (dash-dotted lines) for the n-decane cylinder. Fig. 2.18a correspond to the case without evaporation (see fig. 2.17 left). In this case, the beginning of the recirculation zone is attached to the back of the cylinder (v -u l = 0). On the contrary, in fig. 2.18b, a peak on the Figure 2.17: Velocity streamlines behind the n-decane cylinder and visualisation of vu l = 0 iso-contour. Right: evaporating drop, left: non-evaporating drop. velocity in the opposite direction of the cylinder velocity is observed, shifting the place where the recirculation zone starts and its size. Moreover, there is another velocity peak in the front of the cylinder with the same direction of the cylinder velocity. This changes in the direction of the gas velocity correspond to the Stefan flow. This time, a difference with the reference value of almost 32% is found for the evaporating n-decane cylinder. The results are summarized in the table 2.4. It is worth mentioning that the same phenomenon of the detachment of the recirculation zone from the backside of the interface and its enlargement was also observed by [START_REF] Thamali R Jayawickrama | The effect of stefan flow on the drag coefficient of spherical particles in a gas flow[END_REF] and [START_REF] Thamali R Jayawickrama | The effect of stefan flow on nusselt number and drag coefficient of spherical particles in non-isothermal gas flow[END_REF]. They studied the effect of Stefan flow on the drag coefficient and Nusselt number of a spherical particle in a steady airflow.

Conclusion

An incompressible formalism for simulating two-phase flows with phase change is coupled with a mass conservative interface capturing method (CLSVOF). A volume sink term in the VOF equation, based on surface density, is used to account for the evaporation process. The use of accurate extensions for discontinuous variables and the ghost fluid method allows the description of additional jump conditions at the interface due to heat and mass transfer. These techniques can handle strong temperature gradients while avoiding numerical diffusion and artificial heating within the fluid.

To deal with the velocity jump at the interface, three different methods were described and tested in a configuration with a small Stefan flow and a large Stefan flow. The first two methods consist of an implicit formulation of the velocity jump, and the third method is based on an explicit formulation. Similar observations to those of [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] were made, with a level set method, for the magnitude of spurious currents in the liquid velocity field when the velocity jump term in the Poisson equation accounted implicitly, and the Dirac function is smeared around the interface. These results are improved by including a sharp treatment of the velocity jump term. Then, judging by the results obtained in the static configuration and its potential to be extended to a compressible formalism, the method with the explicit treatment is chosen for the following section, dedicated to compressible simulations.

Finally, several validation cases illustrate the mass conservation properties of our method and the quality of the coupling between vaporization rate and flow dynamics.

Chap. 3 | Development of compressible formalism with heat and mass transfer

Introduction

As we mentioned in the previous chapter, incompressible formalism dedicated to simulations of two-phase flows with phase change is limited to non-confined environments, i.e., at least one outflow boundary condition must be considered.

The velocity divergence created at the interface due to the velocity jump must be evacuated from the domain to maintain the incompressibility assumption. Additionally, to be mass conservative in a closed environment simulation with a vaporizing liquid, the gas density should vary with time depending on the vapor molecular mass. This latter limitation has been addressed in the literature for the incompressible formalism by predicting the gas density at the next time step using the amount of vaporized liquid at the current time step ( [START_REF] Calimez | Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs[END_REF]). Nevertheless, this approach is unable to predict the spatial variation of pressure or density. Therefore, a weakly compressible formalism for simulations of two-phase flows in vaporization in a closed environment is considered. This way, the density and pressure changes induced by the vapor generated at the interface are taken into account, allowing simulations two handle many liquid and gas structures with independent thermodynamic conditions. This chapter presents the governing equations for compressible two-phase flow simulations with heat and mass transfer. Moreover, the numerical schemes implemented to solve them are discussed. As in the previous chapter, the onefluid formalism is considered, where a single system of equations is solved for all phases in the domain. Also, we use configurations in which a single-component liquid evaporates into a mixture of ideal gases. In this case, the two phases are considered compressible, and an equation of state adapted to each phase is chosen to close our system of equations. A mass conservative method of capturing L. Germes Martínez DNS of phase change June 28, 2022 the interface (CLSVOF) is used. This method was adapted to account for the evolution of the liquid mass due to phase change and density variation. Next, the interface capturing method is coupled to a modified projection method for the time evolution of the momentum equation based on two velocities.

A validation case of an evaporating cylinder in an enclosed environment is presented for the compressible formalism. Finally, a 3D Homogeneous Isotropic Turbulence (HIT) configuration is studied to illustrate the ability of the method to handle compressible turbulent two-phase flows with phase change.

Governing equations 3.2.1 Navier-Stokes equation

The Navier-Stokes and the continuity equations for a fully compressible flow can be written as:

       ∂ρ ∂t + ∇ • (ρu) = 0 ∂ρu ∂t + ∇ • (ρu ⊗ u) = ∇ • Ω + ρf v ol (3.1a) (3.1b)
where ρ, u and P are the density, velocity and pressure, respectively. f v ol is the volume forces, Ω = -P Ī + τ is the total stress tensor and Ī the identity matrix. The viscous stress tensor is defined as:

τ = 2µ ¯ - 2 3 µ (∇ • u) Ī (3.2)
where ¯ = 1 2 ∇u + ∇u T . c is the speed of sound, µ the dynamic viscosity. We consider here that both fluids are compressible. Consequently, the system of equations 3.1 must be closed by two equations of state, one for each phase. For the liquid, a modified version of the Tait equation is chosen:

ρ l = ρ 0 P -P 0 B + 1
1 γ l (3.3)
where ρ 0 is a density at a reference state, P 0 is a reference pressure, B is a constant, and function of the isothermal compressibility of the fluid, and γ l is a pressure-independent parameter also inherent to the liquid.

Concerning the gas phase, for isentropic cases, a Laplace's law is chosen:

ρ g = P C γ 1 γg (3.4) L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022
where γ g is the heat capacity ratio and C γ is a constant dependent of the initial condition of the ideal gas. Moreover, for configurations with heat and mass transfer a perfect gas Equation of state is used :

ρ g = P r T (3.5)
where r is the specific ideal gas constant.

Energy and species conservation equation

The energy conservation equation can be expressed in term of the temperature [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]:

ρc p ∂T ∂t + u • ∇T = ∇ • (λ∇T ) + α T T ∂P ∂t + u • ∇P + τ : ∇u + Q (3.6)
where c p is the heat capacity, α T is the coefficient of thermal expansion, T is the temperature, λ is the thermal conductivity, and Q is an energy source term.

Since the liquid phase is considered as mono-component, the species conservation equation is only described for the gas phase:

∂ρ g Y v ∂t + ∇ • (ρ g uY v ) = ∇ • (ρ g D m ∇Y v ) (3.7)
where Y v is the vapor mass fraction and D m is the mass diffusivity.

Interface jump conditions

The governing equations do not directly include phase change terms because the one-fluid formalism adopted in the previous chapter is maintained for the compressible formalism. This means that adapted interface jump conditions must be implemented to respect the conservation of different quantities, such as mass, momentum, energy, and species. The classical jump conditions presented in equations 2.4 remain applicable, with the only difference that the density of each phase varies in time and space. Therefore, in the equations 2.4, 2.8 and 2.13, we have to specify that the densities are at an infinitely small distance from the interface:

[ρ] Γ = ρ Γ l -ρ Γ g (3.8) ρ Γ l u Γ l -u Γ • n = ρ Γ g u Γ g -u Γ • n = ω (3.9) s d = - ω ρ Γ l n (3.10)
The remaining jumps conditions stay unchanged and are summarized in table 3.1.

Quantity

Jump condition

Density ρ Γ l -ρ Γ g Viscosity µ l -µ g Thermal conductivity λ l -λ g Pressure [P -n • (2µ ¯ ) • n] Γ = σκ -ω2 1 ρ Γ Velocity ρ Γ l -ρ Γ g ρ Γ g s d
Heat flux -h lv ω Table 3.1: Summary of the jump conditions for the different quantities.

Pressure equation

When compressible flows with heat and mass transfer are involved, the total pressure depends on several factors, such as thermal expansion, phase change, etc. To simulate this kind of flow, following the pressure variations in space and time is essential to conserve mass, especially in confined environments. This is because the gas density is proportional to the pressure, as we can see from the equations of state (equations 3.4 and 3.5). To obtain an equation for the pressure evolution, [START_REF] Caltagirone | A multiphase compressible model for the simulation of multiphase flows[END_REF] presented a method that consists of reformulating the continuity equation by developing the time derivative of pressure as a function of temperature and density. A similar approach can be found in [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF] and [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF] for isentropic configurations. More recently, a similar formulation of the pressure equation that takes thermal dilatation into account is applied in [START_REF] Urbano | A semi implicit compressible solver for two-phase flows of real fluids[END_REF]. The first step in the development of our pressure equation is to write the derivative of pressure as a function of temperature and density:

DP Dt = ∂P ∂T ρ DT Dt + ∂P ∂ρ T Dρ Dt (3.11)
The first term on the RHS represents the pressure change due to thermal expansion, and the second represents the pressure change due to the compress-ibility effects and phase change. To develop this equation, the definition of several thermodynamic relations is used:

• Isothermal compressibility coefficient (χ T ):

χ T = - 1 V ∂V ∂P T (3.12)
where V is the specific volume (m 3 ). Since the density is inversely proportional to the volume, it is relatively easy to demonstrate that we can write χ T as:

χ T = 1 ρ ∂ρ ∂P T (3.13)
• Isobaric Thermal Expansivity (α T ):

α T = 1 V ∂V ∂T P = - 1 ρ ∂ρ ∂T P (3.14)
• Heat capacity ratio (γ) and Reech's relation:

γ = c p c v = χ T χ S (3.15)
where χ S = 1 ρ ∂ρ ∂P S is the isentropic compressibility coefficient. • General Mayer's relation:

c p -c v = T α 2 T ρχ T (3.16)
• Speed of sound (c) and its relation with the isentropic compressibility coefficient:

c 2 = ∂P ∂ρ S = 1 ρχ S (3.17)
Now, we are going to use the equations 3.13 and 3.14 to determine the coefficient of the second term of eq. 3.11:

α T χ T = -1 ρ ∂ρ ∂T P 1 ρ ∂ρ ∂P T (3.18)
In the following, we use some properties of the partial derivative to obtain a final expression for the coefficient:

α T χ T = -∂ρ ∂T P ∂ρ ∂P T (3.19)

L. Germes Martínez

DNS of phase change June 28, 2022

α T χ T = - ∂ρ ∂T P ∂P ∂ρ T (3.20) α T χ T = ∂P ∂T ρ (3.21)
For the coefficient of the first term in the RHS of eq. 3.11, we take eq. 3.13 and if ∂ρ = 0:

∂P ∂ρ T = 1 χ T ρ (3.22)
To continue the development of our pressure equation, knowing from 3.1a that Dρ Dt = -ρ∇ • u we have:

DP Dt = α T χ T DT Dt - 1 χ T ∇ • u (3.23)
Then, substituting eq. 3.6 into eq. 3.23:

1 - α 2 T T ρχ T c p DP Dt = α T ρχ T c p ∇ • (λ∇T ) + τ : ∇u + Q - 1 χ T ∇ • u (3.24)
Using the general Mayer's relation and the heat capacity ratio (equations 3.16 and 3.15, respectively), we can obtain:

1 - α 2 T T ρχ T c p = 1 γ (3.25)
Finally, with the relation between the speed of sound and the isentropic compressibility coefficient and the Reech's relation, we find

χ T = γ ρc 2
and replacing this into eq. 3.24, the final form of the pressure equation is:

DP Dt = α T c 2 c p ∇ • (λ∇T ) + τ : ∇u + Q -ρc 2 ∇ • u (3.26)

Numerical Methods

The numerical methods implemented for the compressible formalism are discussed in this section. This method can be considered an extension of the incompressible formalism presented in the previous chapter. Several terms are added to the discretized equations to account for the compressible effects of temperature and pressure variations. First, modifications of the interface capturing method are discussed. Second, the projection method presented in the previous chapter, based on the method proposed by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] is modified with the help of the pressure equation (eq. 3.26). Finally, the strategies were implemented for solving the energy and species conservation equations.

Interface capturing

In two-phase flow simulations, the interface is the boundary between the phases. The numerical method used to capture the interface is coupled to the solution of the Navier-Stokes equations by the definition of discontinuous fluid properties, such as density and viscosity. In addition, the identification of each fluid in the domain will also depend on the interface capturing method used. For our compressible formalism, the Coupled Level-set/Volume of Fluid (CLSVOF) method implemented in the ARCHER code [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF] is extended, this time, to handle multiphase flows with phase change and compressibility effects.

In this section, we begin by presenting the VOF equation for the compressible formalism with phase change. With this approach, density variations due to temperature, pressure gradients, and evaporation in closed environments can be represented. To demonstrate the VOF equation, the first step is to define the density as a weighted average between the density of each phase:

ρ = α l ρ l + α g ρ g (3.27)
where α l = V l V cell and α g = Vg V cell are the local liquid and gas volume fraction, respectively. Then, the definition of ρ is introduced into the continuity equation (eq. 3.1a):

∂ (α l ρ l + α g ρ g ) ∂t + ∇ • ((α l ρ l + α g ρ g ) u) = 0 (3.28)
The next step is to separate this equation in two parts, ones for each phase:

       ∂α l ρ l ∂t + ∇ • (α l ρ l u) = - ṁ ∂α g ρ g ∂t + ∇ • (α g ρ g u) = ṁ (3.29a) (3.29b)
By simplifying equations 3.29a and 3.29b and using the conventional nota-tion for the material derivative

D[•] Dt = ∂[•] ∂t + u • ∇[•] , we obtain:        ∂α l ∂t + ∇ • (α l u) = - α l ρ l Dρ l Dt - ṁ ρ l ∂α g ∂t + ∇ • (α g u) = - α g ρ g Dρ g Dt + ṁ ρ g (3.30a) (3.30b)
where ṁ is the volumetric evaporation rate (Kg • m -3 • s -1 ) and represent the amount of liquid evaporated in each cell containing an interface, which explained the negative sign in the RHS in eq. 3.29a. Also, ṁ, in eq. 3.29b represent the amount of gas produced due to the phase change. ṁ is estimated by:

ṁ = ρ l Σ s d (3.31)
where Σ = S V cell is the surface density in a mixed cell, S is the surface, obtained in the PLIC reconstruction. The material derivative of the density can be obtained by writing eq. 3.1a on a non-conservative form and using eq. 3.26 for the velocity divergence:

Dρ l Dt = 1 c 2 l DP Dt - α T l c pl ∇ • (λ∇T ) (3.32)
The pressure and temperature terms of eq. 3.32 are solved after the solution of the Poisson equation for the pressure and during the temporal evolution of the temperature, respectively. Both will be detailed in the following sections.

In our method, only α l is followed by solving eq. 3.30a and α g is calculated using the property:

α l + α g = 1
The convective term of eq. 3.30a is discretized, as for the incompressible formalism, using a modified version of the split advection scheme of the conservative second-order "coupled" operator defined in [START_REF] Sussman | A coupled level set and volumeof-fluid method for computing 3d and axisymmetric incompressible twophase flows[END_REF]. The interface is convected with a liquid velocity field (u l ). This velocity field is continuous and allows us to transport the interface without affecting the velocity jump at the interface. The construction of the liquid velocity field is explained in the next section. The first step of the convection of α l is to advance the liquid volume fraction in each direction without considering the phase change or the compressibility terms:

                   αl -α n l ∆t + (α n l u l ) i+ 1 2 -(α n l u l ) i-1 2 ∆x = αl u l,i+ 1 2 -u l,i-1 2 ∆x αl -αl ∆t + ( αl v l ) j+ 1 2 -( αv l ) j-1 2 ∆y = αl v l,j+ 1 2 -v l,j-1 2 ∆y ᾱl -αl ∆t + ( αl w l ) k+ 1 2 -( αw l ) k-1 2 ∆z = ᾱl w l,k+ 1 2 -w l,k-1 2 ∆z (3.33a) (3.33b) (3.33c)
The computation of the liquid flux across each cell face is detailed in [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]. In a second step, the convection along every direction is coupled in a final equation where the sink and compressibility terms are added:

α n+1 l -ᾱl ∆x = -ᾱl w l,k+ 1 2 -w l,k-1 2 ∆z -αl v l,j+ 1 2 -v l,j-1 2 ∆y -αl u l,i+ 1 2 -u l,i-1 2 ∆x - α n+1 l ρ l Dρ l Dt - ṁ ρ l (3.34)
As in the previous chapter, the signed distance function is used to obtain the geometrical information of the interface, e.g., the vector normal to the interface (n) and the curvature (κ):

n = ∇φ |∇φ| κ = -∇ • n (3.35)
A redistancing algorithm is applied following the same principle of previous work ( [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF], [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF], [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]).

An equation for the velocity divergence can be obtained by adding the equation system 3.30:

∇ • u = - α l ρ l Dρ l Dt - α g ρ g Dρ g Dt + ṁ 1 ρ g - 1 ρ l (3.36)
From eq. 3.36, we can conclude that the velocity divergence depends on several phenomena. For example, in a grid cell located entirely on the liquid or gas side, the velocity divergence is governed by the spatial and temporal changes of the density, which are themselves controlled by the local thermodynamic state, as we can see in eq. 3.32. Furthermore, in the mixed cell, the velocity divergence also depends on the magnitude of the vaporization rate and the density ratio, both of which are represented in the phase change term ṁ 1 ρg -1 ρ l . Moreover, according to eq. 3.36, in the mixed cell, the velocity divergence can be estimated by a weighted average of the compressibility term of each phase following the amount of volume in the cell. A similar approach for the estimation of the velocity divergence in two-phase flows with phase change is found in [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF][START_REF] Poblador | Volume of fluid method for low-mach-number compressible supercritical liquid jet[END_REF].

Projection method

A projection method is used for the time evolution of the momentum equation (eq. 3.1b). Using a projection method decouples the velocity field and the pressure computation. This method consists of an adaptation of the method proposed by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] (dedicated to discontinuous velocity fields at the interface) to account for pressure variation and compressibility effects. An Eulerian staggered mesh is used, where the velocities are defined on the faces and the other variables (e.g., P, α l , T , and Y v ) in the center of the cell. The main idea of the original method is the construction of two velocity fields, the liquid velocity field and its extension in the gas phase and the gas velocity field and its extension in the liquid phase. This allows avoiding the numerical instabilities created when using conventional numerical schemes with a velocity discontinuity at the interface induced by the phase change.

Computation of the intermediate velocities

The first step of the projection method is to compute the intermediate velocities. This is done by solving the momentum equation (eq. 3.1b) without the pressure term. After splitting the temporal derivative of eq. 3.1b, we obtain:

u * = u n -∆t   (u n • ∇) u n - ∇ • 2µ ¯ -2 3 µ∇ • u n Ī ρ n+1   (3.37)
where ρ n+1 is the density computed at the face of the staggered mesh after the resolution of the interface. For the resolution of the viscous term, the method proposed by [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] is implemented. Since this method was originally conceive for incompressible flows simulations, the term 2 3 µ∇ • u n Ī is added to the formalism. Therefore, the new viscous tensor, in matrix notation is: velocity, the final result of the term ∇ • τ has to be located at the faces of the mesh. To this end, attention must be paid when computing the different components of the viscous tensor. For example, in a two-dimensional case, the spatial discretization of the viscous term is:

τ = 2µ ¯ - 2 3 µ∇ • u Ī = 2µux -2 3 µ(ux
∇ • τ =        ( 2µux -2 3 µ(ux +vy +wz ) ) i+1,j -( 2µux -2 3 µ(ux +vy +wz ) ) i,j ∆x + (µ(uy +vx )) i+ 1 2 ,j+ 1 2 -(µ(uy +vx )) i+ 1 2 ,j-1 2 ∆y   i+ 1 2 ,j   (µ(uy +vx )) i+ 1 2 ,j+ 1 2 -(µ(uy +vx )) i-1 2 ,j+ 1 2 ∆x + ( 2µvy -2 3 µ(ux +vy +wz ) ) i,j+1 -( 2µvy -2 3 µ(ux +vy +wz ) ) i,j ∆y   i,j+ 1 2      (3.39)
A harmonic mean that depends on the sign and magnitude of the level-set function is used to compute the viscosity at the cell center and cell faces of the mesh when an interface is present. Moreover, the convective terms of eq. 3.37 are computed with a 5th order WENO scheme.

In our compressible formalism, an intermediate velocity u * is calculated for both liquid and gas fields u * l , u * g by repeating the operation for each phase:

                 u * l = u n l -∆t   (u n l • ∇) u n l - ∇ • 2µ ¯ l -2 3 µ∇ • u n l Ī ρ n+1   u * g = u n g -∆t   u n g • ∇ u n g - ∇ • 2µ ¯ g -2 3 µ∇ • u n g Ī ρ n+1   (3.40a) (3.40b) Pressure solver
The next step of the projection method is to determine the pressure field. For this purpose, a Poisson equation is derived from the following equation:

u n+1 = u * -∆t ∇P n+1 ρ n+1 (3.41)
The divergence operator is then applied to both sides of eq. 3.41:

∇ • u n+1 = ∇ • u * -∇ • ∆t ∇P n+1 ρ n+1 (3.42)
Then, we insert eq. 3.42 into the pressure equation (eq. 3.26), resulting in the following Helmholtz equation for the pressure:

-∇ • ∇P n+1 ρ n+1 + 1 ρc 2 ∆t 2 P n+1 = 1 ρc 2 ∆t 2 (P n -∆tu • ∇P n ) + ∇ • u * ∆t + α T ρc p ∆t ∇ • (λ∇T ) (3.43)
This implicit temporal discretization of the acoustic term removes the acoustic CFL restriction, allowing larger time steps. Furthermore, the convective term of the pressure (u • ∇P ) is discretized with a 5th order WENO scheme. In addition, the Poisson solver used for this Helmholtz type equation consists of an MGCG (MultiGrid preconditioned Conjugate Gradient) method coupled with the Gauss-Seidel "Red-Black" iterative scheme.

The acoustic terms of the equation 3.43 are calculated as a weighted average depending on the volume of each fluid in the cell, as suggested by [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF]:

1 ρc 2 ∆t 2 P n+1 = α l ρ l c 2 l ∆t 2 + α g ρ g c 2 g ∆t 2 P n+1
(3.44)

1 ρc 2 ∆t 2 (P n -∆tu • ∇P n ) = α l ρ l c 2 l ∆t 2 + α g ρ g c 2 g ∆t 2 P n + α l ρ l c 2 l ∆t u l • ∇P n - α g ρ g c 2 g ∆t u g • ∇P n (3.45)
In eq. 3.45, we can see that the liquid and gas velocities are used depending on the magnitude of α l . Since u l and u g are continuous at the interface, numerical instabilities are avoided when implementing the 5th order WENO scheme due to the velocity discontinuity. In addition, the second RHS term of eq. 3.43 is estimated, as the incompressible formalism presented in the previous chapter, following the sign of the level set equation:

∇ • u * ∆t =    ∇•u * l ∆t if φ > 0, ∇•u * g ∆t if φ ≤ 0, (3.46)
The energy term on the RHS of eq. 3.43 is also estimated following the sign of the level-set function, and ∇ • (λ∇T ) is taken from the computation of the energy equation diffusive term. This formulation allows us to be consistent with the implementation of the energy jump due to the phase change; which is explained in the following section:

α T ρc p ∆t ∇ • (λ∇T ) =    α T l ρ l c pl ∆t ∇ • (λ∇T ) if φ > 0, α T g ρgcpg∆t ∇ • (λ∇T ) if φ ≤ 0, (3.47)
For computing the velocity field in the next time step (n + 1), we use eq. 3.41. And repeat the same operation for the liquid and gas velocity using the total pressure field:

   u n+1 l = u * l -∆t ∇P n+1 ρ n+1 if φ > 0, u n+1 g = u * g -∆t ∇P n+1 ρ n+1 if φ ≤ 0, (3.48)
Velocities and its extensions

To use the GFM, the next step is to compute a liquid and gas velocity extensions.

Starting with the liquid velocity extension, a ghost pressure (P ghost ) is computed:

-∇ • ∇P ghost ρ n+1 + 1 ρc 2 ∆t 2 P ghost I = ∇ • u * l ∆t + α T ρc p ∆t ∇ • (λ∇T ) II (3.49)
Eq. 3.49 is solved following the same methods used for eq. 3.43. Concerning the incompressible formalism, the new term I is related to flow acoustics, and the new term II is the contribution of thermal dilatation. These two terms have been added compared to the original projection method proposed by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. The liquid velocity extension (u ghost l ) is computed from P ghost :

u ghost l = u * l -∆t ∇P ghost ρ n+1 (3.50)
Since both phases are considered to be compressible, for this reason, the steps dedicated to enforce the divergence-free condition in the liquid velocity field and its extension are skipped. It is worth mention that the liquid velocity extension u ghost l is continuous but non divergence-free. And the liquid velocity in n + 1 is obtained from the liquid velocity field (eq. 3.48) and its extension (eq. 3.50): we add the interface velocity jump to the liquid velocity in the liquid phase:

u n+1 l =    u n+1 l if φ > 0, u ghost l if φ ≤ 0,
u ghost g = u n+1 l -s d ρ l -ρ g ρ g (3.52)
Finally, the gas velocity field is obtained, following the sign of the level set function:

u n+1 g =    u ghost g if φ > 0, u n+1 g if φ ≤ 0, (3.53)

Energy transport

For the temporal evolution of the temperature field, in the energy equation (eq. 3.6), the term responsible for the energy generation by friction ( τ : ∇u) and the term of volumetric energy generation Q are neglected because the production or consumption of energy by a chemical reaction in the system is not considered, yielding:

ρc p ∂T ∂t + u • ∇T = ∇ • (λ∇T ) + α T T ∂P ∂t + u • ∇P (3.54)
Eq. 3.54 is solved in two-step. First, an intermediary temperature (T * ) is computed containing the convective and diffusive terms. Then, the pressure term is added after the resolution of the Poisson equation(eq. 3.43).

The convective terms of eq. 3.54 are discretized using a 5th order WENO scheme. Large temperature gradients can occur at the interface, which is a source of error and can lead to numerical heating. To avoid this, the temperature field of each phase is extrapolated using Aslam's extension method [START_REF] Tariq D Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF]. The idea is to extend the scalar of one phase linearly, following the normal direction of the interface, into the other, populating the ghost grid points. Then a 5th order WENO scheme is used for convection of the liquid temperature and its ghost values (u l • ∇T l ) and the same procedure is applied for the gas side (u g • ∇T g ). Finally, these terms are introduced into eq. 3.54 depending on the sign of the Level Set function.

A 2nd order central difference scheme is used to discretize the diffusion term. When the mass transfer is considered, the energy jump at the interface (eq. 2.20) must be taken into account in eq. 3.54. For this reason, the numerical scheme for the diffusion term is modified in the grid points close to the interface using a GFM. The GFM used for the compressible formalism is similar to the one presented in section 2.2.2.

Vapor mass fraction transport

For the transport of the vapor mass fraction (Y v ), eq. 3.7 is solved only in the gas phase. Eq. 3.7 reduces to:

∂Y v ∂t + (u g • ∇) Y v = ∇ • (ρ g D m ∇Y v ) ρ g (3.55)
As for eq. 3.54, the diffusion term is discretized using a 2nd order central difference scheme, and the convective term is calculated with a 5th order WENO scheme. To apply the Dirichlet boundary condition at the interface (Y v s ), a 1st order Aslam-Chiu extension method is implemented. More details about this method are available in [START_REF] Bouali | Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows[END_REF]. The interface temperature (T Γ ) needed for eq. 2.22 is estimated with eq. 2.57.

Regarding the calculation of the time step, a CFL condition similar to [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF] is used. To summarize the steps of the numerical methods, fig. 3.1 shows the algorithm used in an Euler time step for the compressible two-phase flow solver.

Simulation results

This section presents and discusses the results from simulations with our compressible formalism for two-phase flows. The first simulations are devoted to validating mass conservation in a closed environment. We then show the ability of our formalism to handle atomization regimes using a two-phase homogeneous and isotropic turbulence configuration.

Validation: total mass conservation in enclosed environment

The conservation of total mass in a closed environment in the presence of a vaporizing liquid can be validated by using the same configuration used in section 2.4.1 (static water cylinder configuration). In this kind of configuration, the amount of liquid mass will decrease with time following the magnitude of the evaporation rate. As a result, the gas mass has to increase because there is no mass transfer between the domain and its surroundings, as in the incompressible two-phase simulations presented in the previous chapter. To quantify the mass transferred in a single time step due to evaporation, a mass balance in the domain after the evaporation step is performed: where m total = m l + m g is the total mass of the domain, which should remain constant over time. And

m total = ρ l (α l -α ev ap ) dV + ρ g (1 -α l + α ev ap ) dV (3.56)
α ev ap = V ev ap V cell = ∆tΣ s d
is the liquid volume fraction vaporized during a time step in a mixed cell. Then, comparing eq.3.56 with a mass balance carried out before the vaporization step:

m total = ρ l α l dV + ρ g (1 -α l ) dV (3.57)
we can obtain the amount of mass transferred between the phases:

∆m = (ρ g -ρ l ) α ev ap dV (3.58)
From eq. 3.58, we can conclude that the amount of liquid mass transferred to the gas is directly related to the density difference between the phases. In our compressible formalism, the increase in mass in the gas phase is accounted for by adding the velocity jump in the ghost cells of the gas velocity field, which will cause an increase in pressure during the solution of Poisson's equation (eq. 3.43). This will increase the gas density, estimated in this case by eq. 3.4. There is also an increase in the gas phase volume due to the liquid's vaporization. Moreover, this increase is taken into account through the sink term in the advection equation VOF (eq. 3.30a) and the fact that α g = 1 -α l .

To simulate an enclosed environment, periodic boundary conditions are employed in all directions. This configuration cannot be performed with an incompressible solver because the Stefan flow generated at the interface should be evacuated by at least one outflow boundary condition. The initial cylinder radius is R D = 200 µm. A 2D rectangular domain is considered with a side length of l x = 24R D . Simulations are performed with a 512 2 grid resolution. In this configuration, the vaporization is imposed at ω = 1 kg m 2 s . Consequently, this subsection does not resolve the temperature and species transport equation.

The temporal derivative of the total gas mass can be estimated ∂mg ∂t = ωS. Consequently, an estimation of the gas mass at a given instant n + 1 (m n+1 g ), including the gas mass increase induced by the vaporization is written as :

m n+1 g = m n g + ωS n ∆t (3.59)
where m n g is the total gas mass at the previous time step of the simulation. Note that the result is given by eq. 3.59 are not theoretical because the local interface area S Γ from the PLIC reconstruction is extracted in each cell and then summed to obtain the total surface S. Fig. 3.2 shows the evolution of the normalized gas mass (m * g = mg(t) m 0 g ) obtained in the simulation where m 0 g is the initial total gas mass, and compared with the reference results given by eq. 3.59. A good agreement is found between the simulation results and the reference, which validates the coupling between the velocity jump term in the Poisson pressure equation and the phase change sink term in the VOF equation. For completeness, Fig. 3.3 shows the simulation results for the evolution of the liquid mass and the reference results computed in a similar way with:

m n+1 l = m n l -ωS n ∆t (3.60)
Furthermore, as previously specified, the total mass in an enclosed environment should remain constant during the simulation to respect continuity. In this configuration, the total mass reduces by less than 0.02% at the end of the simulation. The same simulation is repeated for 3 others mesh resolutions (See Table 3.2). This result demonstrates the excellent mass conservation properties of our CLSVOF interface capturing method.

L. Germes Martínez

Compressible HIT configuration with phase change

Similar to previous works, the idea is to study the influence of the evaporation phenomenon inside a forced HIT two-phase flows. In [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], a passive scalar was used to represent the evaporation and mixing process. Among the limitations of this procedure, the influence of the Stefan flow in the velocity field and the pressure increase due to the evaporation process was not considered. A significant improvement has been proposed by [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF] where a phase change source term was directly introduced in the continuity and pressure equations. The results obtained illustrated good mass conservation properties. However, the influence of the temperature in the vapor mass fraction at the interface was not implemented yet to the general formalism: the vaporization rate was directly imposed. Contrary to [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF], here, the velocity jump is introduced implicitly to the Poisson equation for the pressure in the ghost points of the gas velocity field. The coupling between the energy, the continuity, and momentum equations is added by assuming thermodynamic equilibrium at the interface.

The thermodynamic properties are those of n-decane (table 2.3). The initial gas density is ρ g = 25 kg • m -3 . The pressure is P = 4.11 × 10 6 P a and for the Tait equation, the parameters used are B = 10 9 P a, P 0 = 10 5 P a, ρ 0 = 750 kg • m -3 and γ l = 1.215. The initial liquid volume fraction is Φ = 10% and the targeted mean kinetic energy is k = 3.6 m 2 •s -2 . The initial liquid temperature is T l = 340 K and the initial gas temperatures is T g = 573 K. The domain is a cube with a side length of 1.5 × 10 -4 m and a 256 3 mesh has been retained for this simulation. Fig. 3.4 shows an instantaneous temperature field and the iso-contour of the level-set function, which represents the liquid/gas interface. As well as in other compressible and incompressible HIT configurations presented in previous works, many breakup and coalescence events are observed. Also, small gas structures are observed inside the liquid; these structures are characterized by an increase in pressure due to the creation of vapor and a decrease in temperature because of the energy consumption induced by the phase change. Additionally, the influence of the convection due to the turbulent velocity field in the temperature field is observed. Similar phenomena are found in the vapor mass fraction field (Fig. 3.5).

The velocity divergence is shown in fig. 3.6 (top) to observe the compressibility and dilatation effects on the velocity field. By comparing Fig. 3.6 (top) and Fig. 3.4, an increase in the velocity divergence has been observed in the presence of temperature gradients in the gas phase. This phenomenon can be explained by rewriting Eq. 3.26 to obtain the equation of the velocity divergence:

∇ • u = - 1 ρc 2 DP Dt + α T ρc p (∇ • (λ∇T )) (3.61)
where the second term is thermal dilatation. This term is dominant in the velocity divergence budget in the largest gas structure, meaning that the acoustic/compressible effects (first term) have a lower order of magnitude in this structure. However, the divergence inside encapsulated gas structure is generated differently: compressible effects have a major role in this scenario. The Stefan flow compresses the gas in a very small volume, leading to an increase in the divergence and gas density. This effect is illustrated in Fig. 3.6 (bottom): the divergence inside the encapsulated gas structure is almost one order of magnitude larger than the maximum divergence in the main gas structure. The temperature inside this structure is quasi uniform. Peaks of velocity divergence are also observed close to the interface due to the presence of the Stefan flow (velocity jump).

Conclusion

The development of a pressure-based method for compressible multiphase flows with phase change is presented. Previous numerical methods implemented for incompressible formalism are extended to a weakly compressible formalism. In our formalism, the compressible Navier-Stokes equations and additional terms related to the compressibility of the flow, such as thermal expansion or the presence of a velocity divergence, are considered. In addition, the VOF equation is modified to consider liquid compressibility, variable density, and the evaporation process. A validation case is investigated in a static cylinder configuration evaporating in a closed environment. Then, a comparison is made between the results obtained with the compressible framework and a reference equation representing the mass balance in a closed environment, showing the accuracy of the formalism.

Finally, the latter formalism is used to simulate a 3D compressible HIT configuration and demonstrate the ability of our method to handle strong interface curvature, variable density, thermal expansion, encapsulation of the gas structure with its thermodynamic pressure, collisions, and breakups.

Chap. 4 | Towards a complete description of mixing and evaporation processes

Introduction

In the literature, the scalar turbulent mixing is studied in several scenarios, such as single-phase flow, dispersed and dense two-phase flow [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF][START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF][START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF][START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]]. In the first example, the temporal evolution of the scalar variance is characterized by an exponential decay due to the scalar dissipation rate [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF]; when a mean scalar gradient is not imposed. This behavior changes drastically when a twophase flow configuration with phase change is considered. In this case, an initial increase of the scalar variance is observed until a maximum value is reached, followed by a quasi-exponential decay is observed. [START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF] explained this behavior by analyzing the scalar variance transport equation, where different source terms appear when there is a vaporizing liquid in the domain. Here, the temporal evolution of the scalar variance results from an equilibrium between the scalar dissipation rate and the vaporizing source term. In the work of [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], the scalar turbulent mixing in a two-phase HIT configuration with a constant vapor mass fraction at the interface is investigated. In this chapter, a continuation of this study is performed using our compressible formalism, where the coupling between vapor mass fraction and temperature is solved. Thus, a non-homogeneous representation of the vapor mass fraction at the interface is possible. First, several improvements to the governing equations and the numerical method described in 3 are presented and discussed. For instance, the influence of the vapor molecular mass is accounted for in the computation of the gas density and heat capacity. Then, the numerical methods devoted to solving the governing equations with the corresponding interface jump conditions and how to compute surface averaged quantities are presented.

Here, scalar mixing and evaporation in an atomization environment are analyzed using the compressible two-phase HIT configuration. To this end, we perform two simulations with different evaporation regimes. Furthermore, the evaporation regime effects on the temporal evolution of the mean vapor mass fraction and the mean liquid and gas temperatures are discussed. Then, we analyze the scalar turbulent mixing by considering the temporal evolution of the vapor mass fraction main statistics. Finally, the behavior of the evaporation rate and the vapor mass fraction at the interface is studied.

Governing equations

This section considers a two-phase flow composed of a liquid and a mixture of ideal gases. Both phases are fully compressible and are separated by an infinitely thin interface. In addition, heat and mass transfer between the phases through the interface is possible. The gas phase is initially described as dry air until conditions necessary for liquid evaporation are met. At that time, the air begins to mix with the vapor produced at the interface. The liquid is considered a singlecomponent fluid with physical properties similar to n-decane. The following system of equations is solved for both phases:

                       ∂ρ ∂t + ∇ • (ρu) = 0 ∂ρu ∂t + ∇ • (ρu ⊗ u) = ∇ • Ω + ρf v ol ρc p DT Dt = ∇ • (λ∇T ) + α T T DP Dt ∂ρ g Y v ∂t + ∇ • (ρ g uY v ) = ∇ • (ρ g D m ∇Y v ) (4.1a) (4.1b) (4.1c) (4.1d)
A forcing method maintains the mean turbulent kinetic energy near a prescribed value. To implement it, the method described by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF] is used. Here, the volume forces term in eq. 4.1 (f v ol ) is used as a source term defined as

f v ol = Au
where A is a forcing coefficient and u represent the velocity fluctuations.

The thermal conductivity (λ), thermal dilatation coefficient (α T ), and the mass diffusivity (D m ) are considered as constants in each phase. However, the heat capacity is considered constant in the liquid, and it is a function of the vapor L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022 mass fraction in the gas, which is defined as a weighted average:

c pg = k=a,v Y k c pk (4.2)
Since both of phases are compressible, two equations of state are needed for closing the system 4.1, i. e., one for each phase. For the liquid, a modified version of the Tait's equation is chosen:

ρ l = ρ 0 P -P 0 B + 1 1 γ l (4.3)
where ρ 0 is a density at a reference state, P 0 is a reference pressure, B is a constant function of the isothermal compressibility of the fluid, and γ l is a pressureindependent parameter, inherent to the fluid. And for the gas, the ideal gas equation of state is used:

ρ g = P M RT (4.4)
where R is the ideal gas constant. In this work, a significant influence of evaporation on the gas properties is expected since the initial volume fraction φ = V l V box is equal to 10 % and the initial temperature gradient between the phases is relatively large. For these reasons, the molar weight of the gas mixture (M) defined

as M = k=a,v Y k M k -1
allows the vapor composition to be considered in the gas density prediction. For the evolution of the pressure, we use eq. 3.26, and the energy jump condition summarized in table 3.1 are used.

Numerical methods

Interface representation and projection method

The system of equations 4.1 is similar to the equations presented in the section 3.2. Hence, similar numerical methods are used to solve them. However, some changes are introduced compared to the previous chapter. For instance, when the liquid begins to evaporate, the heat capacity of the gas phase depends on the vapor mass fraction and the heat capacity of each species. Also, the gas molecular mass present in the EOS (eq. 4.4) undergoes similar changes. This way, the physical properties of the vapor are taken into account, bringing us closer to a more realistic model of the evaporation phenomenon.

Surface averaged quantities

This section is devoted to describing the estimation of the variables at the liquid/gas interface. 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 local surface area and mean curvature are extracted. To this end, we use the "sur f ace_oper ator s" 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 gas-liquid interface using geometrical variables: toward a unified description of the disperse and separated phase flows[END_REF] and [START_REF] Di Battista | Post-processing of two-phase dns simulations exploiting geometrical features and topological invariants to extract flow statistics: application to canonical objects and the collision of two droplets[END_REF] 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 iso-surface, respectively. Finally, the surface averages quantities are calculated as follows:

ψ = 1 S r ψdS r (4.5)
where ψ = ω, Y v and S r is the surface area reconstructed from the zero level-set surface.

Numerical configuration

This section is devoted to describing the numerical configuration used for the study of evaporation and turbulent scalar mixing. Since one objective is to investigate the influence of the evaporation rate on the time evolution of the mentioned processes. First, the values of the equilibrium temperature and the vapor mass fraction of the system are estimated. It allows us to predict whether there is a significant difference in the final state of the system when an initial condition is changed. Thus, this helps us to choose the initial parameters of the configurations. Because we consider that when the final state of the system is not the same, the time evolution of the different quantities also changes. For this purpose, an enthalpy balance is solved by an iterative algorithm. In the following, we describe the two compressible two-phase HIT configurations and the simulation initialization procedure.

Equilibrium Conditions

Phase equilibrium is achieved when there is no mass transfer between the phases [START_REF] Yunus | Thermodynamics: an engineering approach[END_REF]. Consequently, the system has reached phase equilibrium with no vapor mass fraction gradient, which is the driving force for mass transfer. In this work, the word equilibrium refers to the phase equilibrium, and it is reached when the vapor mass fraction in each point of the gas phase is equal to the vapor mass fraction at the interface. Since the vapor mole fraction at the interface is related to the temperature via the Clausius-Clapeyron relation (eq. 2.22), thermal equilibrium is reached when the system reaches phase equilibrium.

An enthalpy balance is performed to estimate the simulation conditions at the equilibrium state. If we considered an isolated system where three components (air, vapor, and liquid) can coexist at the same time, we can write:

∆Q l + ∆Q g + ∆Q v = Q lv (4.6)
where ∆Q l , ∆Q g and ∆Q v represent the enthalpy change (J) between the initial and the equilibrium state of the liquid, gas and vapor phase, respectively. When phase change occurs, the enthalpy change in the system is caused by the latent heat, which is accounted in the eq. 4.6 by the term Q lv . The temperature at the equilibrium (T eq ) can be determined from eq. 4.6, by using the definition of the sensible enthalpy (Q = mc p ∆T ), the latent heat (Q = mh v l ), and the following assumptions:

• Heat capacities of the components 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 ev ap h lv m l,eq c pl + m 0 a c pa + m v ,eq c pv (4.7)
where 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 . The amount of evaporated mass can be determined by m ev ap = m v ,eq -m 0 v .

To determine the equilibrium conditions, the liquid mass at the equilibrium (m l,eq ) is varied until eq. 4.7 is satisfied. To this end, similar to [START_REF] Palanti | An implicit formulation to model the evaporation process in the eulerian-lagrangian spray atomization (elsa) framework[END_REF], an iterative algorithm is implemented.

Fig. 4.1 shows the algorithm used to determine T eq . First, the initial variables (T 0 l , T 0 g , m 0 l , m 0 g and m 0 v ) are introduced and a test is done to verify that the evaporation is possible at this condition. By assuming zero evaporation (m l,eq =

X eq = n v ,eq n a + n v ,eq
where n v ,eq = mv,eq Mv and n v ,eq = m 0 a Ma are the numbers of vapor and air moles at the equilibrium, respectively. Also, the vapor mole fraction at the saturation (X th = f (T eq )) is computed using equation 2.22. If X eq > X th , the system is saturated, thus evaporation is no possible given this initial condition. As a result, only heat transfer can occur in the domain and a corresponding T eq is computed.

Otherwise, if X eq < X th , a second test is done for the total evaporation case (m l,eq = 0). Here, if X eq > X th , the initial conditions allow the complete evaporation of the liquid mass and T eq is the final gas temperature. In the contrary, partial evaporation of the liquid is possible. In this case, T eq is determined using a bisection method where m l,eq is varied in the interval [0, m 0 l ] until the relative error e = |Xeq-X th | X th is less than an arbitrary tolerance tol = 10 -10 . Finally, the vapor mass fraction at the equilibrium (Y v ,eq ) is computed with:

Y v ,eq = X eq M v ap X eq M v ap + (1 -X)M g (4.8)

Parameters choice

The main objective of this chapter is to study scalar mixing in a turbulent environment with the presence of an evaporating liquid. Here, the interface acts as an inhomogeneous vapor mass fraction source and, at the same time, acts as a sink for the temperature. Similar to previous works ( [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF][START_REF] Martinez | A new dns formalism dedicated to turbulent two-phase flows with phase change[END_REF]), a compressible two-phase HIT configuration with phase change is considered. Simulations are carried out in a three-dimensional box with a length size l x = 1.5 × 10 -4 m and a mesh of 256 3 . In addition, periodic boundary conditions are considered in each direction to simulate an enclosed environment. A forcing method is implemented to keep the mean kinetic energy (k) at a value close to 3.6 m 2 • s -2 . The liquid has similar physical properties as n-decane, and the gas is initially considered dry air. Then, when the evaporation of the liquid 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 4.1. The parameters for the Tait's equation are B = 10 9 P a, P 0 = 10 To investigate the influence of the evaporation regime on scalar turbulent mixing, two simulations are performed in which only the initial liquid temperature is changed. The reason is that the vapor mass fraction at the interface depends on the difference between the interface and boiling temperatures (see equation 2.22 and 2.23). The initial liquid temperatures of the simulations are T 0 l = 400 K, 340 K and the initial gas temperature for both cases is T 0 g = 573 K. The first simulation can be considered a large evaporation case compared to the second simulation since the liquid temperature is closer to the boiling point. The equilibrium state and temporal evolution of temperature and vapor mass fraction are expected to be different for both simulations due to the variation of the evaporation regime. Using the algorithm presented in section 4.4.1, the variables at the equilibrium are calculated for each case. The results are summarized in the table 4.3. In a hypothetical case where evaporation is not considered in the system (only heating), the temperature at equilibrium is somewhere between the initial liquid and gas temperatures. In the table 4.3, we can observe that the temperature at equilibrium for case 1 is lower than the initial temperature of the liquid phase. This is explained by the large evaporation rate, where the liquid's latent heat consumes an important part of the initial enthalpy of the system. Consequently, the system has less enthalpy than the initial state at the equilibrium state. Different results are observed for case 2, where the equilibrium temperature is higher than the initial temperature of the liquid. Since case 2 has a lower evaporation rate, the system's enthalpy in the initial state is sufficient to evaporate the liquid and increase the liquid temperature. However, the increase in liquid temperature is lower than in the case without evaporation because some enthalpy is still consumed.

Cases T 0 l (K) T 0 g (K) T eq (K) Y v ,
Moreover, the impact of the evaporation regime can be seen in the amount of liquid evaporated. In case 1, the vapor mass fraction at equilibrium is 2.4 times the vapor mass fraction of case 2. This indicates a more significant amount of evaporated liquid in case 1, caused by the increased evaporation rate. Another parameter we can analyze is the liquid mass at equilibrium normalized by the initial liquid mass. 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 ) of the simulations. Thus, the comparison between the different cases investigated is more pertinent. In addition, to analyze different stages of the turbulent mixing 4 times are defined (t 1 , t 2 , t 3 and t 4 ), each correspond to a specific value of the normalized mean vapor mass fraction (Y v * = 0.18, 0.32, 0.54 and 0.72). It is important to note that the physical time represented by each defined instant is not the same for the two simulations due to their different paths towards the equilibrium state.

In addition, the values of Y v * are chosen taking into account the time required to reach the equilibrium state, and the computational cost involved. Nevertheless, at least 20 eddy turnover time is performed in 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 flows[END_REF][START_REF] Bouali | Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows[END_REF]. Here, the level-set function and velocity fields are initialized in a cubic domain with periodic boundary conditions in each direction, based on the work of [START_REF] Luret | Modeling collision outcome in moderately dense sprays[END_REF]. It consists of eight drops, equally spaced in the domain. The sum of their volume is equal to the initial liquid volume fraction. Each droplet has a rotational velocity corresponding to the average kinetic energy.

At the beginning of the simulation, the linear forcing, species, and energy transport equations are deactivated; and the initial liquid and gas temperatures are kept constant by the temperature correction step. Two additional tempera-ture fields are created in this step, one for the liquid temperature and its ghost points in the gas and the other for the gas temperature and its ghost point in the liquid. Aslam's extension method is used to populate the ghost points of each temperature field. Then, the temperature correction step identifies when a mesh point shifts from one phase to the other and assigns its ghost value from the previous iteration.

Once the domain is destabilized, the linear forcing is activated, and the turbulent kinetic energy reach its prescribed level. Afterward, diffusive sub-iterations 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 (4.9a) (4.9b)
where k is the current sub-iteration.

The number of diffusive sub-iterations (k end ) is computed by:

k end = t BL ∆t dif f (4.10)
where t BL = δ Y v s dif f is the time necessary to generate a boundary layer close to the interface. More precisely, this procedure move a vapor mass fraction isocontour 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 diffusion restrictions:

V cf l = 6v isc dx 2
where v isc = max λg ρgcpg , λ l ρ l c pl , D m . This procedure avoids numerical instabilities created by high temperature and vapor mass fraction gradients in the initial fields. Finally, the energy and vapor mass fraction transport equations (equations 4.1c, 4.1d) are activated and the simulation statistics are collected from this point.

Results

Heat and mass transfer

The instantaneous temperature fields at the domain boundaries for the two cases with interface visualization (zero level-set isocontour) at t 2 are presented in L. Germes Martínez At this point in the simulations, the effects of the turbulent velocity field are easily appreciated in the temperature transport. As in the previous chapter, many interactions between the liquid and gas structures induced by the velocity field are observed. The opacity of the interface visualization is reduced to show the bubbles enclosed within the larger liquid structures. Due to the reduced gas volume, these bubbles are characterized by reaching equilibrium states at shorter times than the rest of the domain. The influence of the gas volume on the time required to reach equilibrium is investigated in more detail by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. In addition, the smallest temperatures are found near the interface, where a wide range of thermal boundary layer thickness is observed. These variations in thickness depend on the deformation of the interface and the temperature of the surrounding gas. For example, as noted in section 2.4.2, 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.

By comparing the figures 4.2, a decrease in the maximum temperature in case 2 is observed. Also, there is a larger homogenization of the gas temperature far from the interface in fig 4 .2b. This indicates a higher level of mixing for the same value of Y * v in case 2 concerning case 1. Similar phenomena are observed in fig. 4.3 for the instantaneous vapor mass fraction fields at the domain boundaries at t 2 . Compared to the temperature fields, larger gradients are present in the gas phase away from the interface. This is because the mass diffusivity (D m = 3.82 ×10 -7 m 2 •s -1 ) is smaller than the thermal diffusivity (D th = 1.70 × 10 -6 m 2 • s -1 ) resulting in a greater influence of the convection term in the vapor mass fraction transport equation compared to the diffusion term. Furthermore, since the vapor mass fraction at the interface is temperature dependent, the evaporation rate is also inhomogeneous. Recently, similar observations have been made for a single evaporating droplet in a homogeneous and isotropic turbulent environment ( [START_REF] Michael S Dodd | Analysis of droplet evaporation in isotropic turbulence through dropletresolved 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 have shown the variations of the evaporation rate due to small deformations caused by the turbulent velocity field. In addition, their results also have shown an increase in the evaporation rate due to the interaction of the vapor mass fraction with the turbulent velocity field. These results highlight the limitation of classical vaporization models based on spherical droplets with a constant evaporation rate. por mass fraction and the turbulent velocity field within the domain. Also, the vapor accumulation is observed with the tendency to be completely homogeneous when the system reaches equilibrium. Another interesting phenomenon is the vapor wake behind the smaller liquid structures as they are convected through the domain. The same behavior is observed in the simulations presented in the section 2.4.2. There are several possible outcomes for the lifetime of these smaller liquid structures, e.g., they may remain convected by the velocity field until they evaporate completely or collapse 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, the distance to another interface, etc.

The temporal evolution of the volume averaged liquid (T l ) and interface T Γ temperature normalized by the equilibrium temperature for both cases are presented in figures 4.5. t * = t τ is the time normalized by the eddy turnover time defined as:

τ = k ξ
where (•) is the volume averaging operator, and ξ is the mean kinetic energy dissipation rate. As predicted, the temperature of the liquid for case 1 decreases as it approaches the equilibrium state (see fig. 4.5a). Also, 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 GFM implemented to account for the energy jump condition and the coupling between the liquid and gas temperature. At the beginning of the simulation, the interface and liquid temperatures have different magnitudes due to the initial diffusion sub-iterations. Then, after evaporation is activated, a sudden drop in temperature is observed before starting to stabilize after two eddy turnover times. This is because the magnitude of the evaporation rate reaches a maximum value early in the simulation. Indeed, the mass boundary layer is not fully formed yet, causing an increase in the vapor mass fraction gradient at the interface. Subsequently, the vapor generated at the interface increases the saturation level, 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 thermal equilibrium has not yet been reached. An opposite behavior is observed for the liquid temperature in case 2. Fig. 4.5b shows an increase in the liquid temperature caused by the initial temperature gradient between the liquid and gas temperature. Similar to case 1, the liquid and interface temperatures have different values due to the initial subiterations. Then, the temperature rise slows down as the heat flux from the gas to the liquid decreases. Resulting in stabilizing the temperature, represented by a plateau after about 15 eddy turnover time. By comparing the two images in fig. 4.5, we can notice the effects of the evaporation regime on the characteristic time required to reach thermal equilibrium.

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 fig. 4.6. If we compare the cases, a significant difference between the initial magnitudes of ωΓ is observed, which is explained by initial liquid temperatures. Moreover, both curves have the same general behavior; a maximum value at the beginning of the simulation followed by a decrease as the system approaches the equilibrium state. However, each case takes very different paths to the equilibrium state. For example, a sharper decrease in the evaporation rate in case 1 is observed, indicating a quicker vapor accumulation near the interface. Furthermore, the curve of the temporal evolution of the evaporation rate for case 1 shows more variations than the curve corresponding to case 2. Speculatively, these variations are suspected to be related to the surface density, which depends on the evolution of the mass of the liquid and the rate of energy dissipation.

L
The normalized mean gas temperature (T * g = T g Teq ) as a function of time is shown in fig. 4.7. Similar to the liquid and the interface temperature, the gas temperatures in fig. 4.7 are not the same because the equilibrium temperatures differ. Besides, a more pronounced temperature decrease is observed in case 2. The effect of the energy jump and the Stefan flow on the heat transfer across the interface could explain this behavior. For example, the temperature decrease at points near the interface caused by latent heat results in a local decrease in the temperature gradient. In other words, the increased evaporation rate slows down the heat transfer. This also explains the increase in the characteristic time required to reach the equilibrium state. Additionally, the heat transfer decrease due to the Stefan flow has been observed by [START_REF] Thamali R Jayawickrama | The effect of stefan flow on the drag coefficient of spherical particles in a gas flow[END_REF][START_REF] Thamali R Jayawickrama | The effect of stefan flow on nusselt number and drag coefficient of spherical particles in non-isothermal gas flow[END_REF] in previous works. Another way to investigate the influence of the evaporation regime on the evolution of a system from the initial condition to the equilibrium state is to study the time evolution of the normalized vapor mass fraction (Y * v = Yv Yv,eq ), where Y v is the volume averaged vapor mass fraction in the gas phase. This is shown in fig. 4.8 for the two cases studied. Here, we can observe a deceleration of Y * v with a decrease in the evaporation rate. These results are consistent with our previous analysis of the temporal evolution of the mean liquid and gas temperatures and evaporation rate. For example, higher evaporation rates indicate faster vapor generation at the interface, contributing to the rapid increase in the normalized mean vapor mass fraction in case 1. In addition, the increased Stefan flow pushes the vapor in a direction normal to the interface, increasing the rate of change of Y * v . 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 field and the interaction between the liquid and gas structures. This last part is investigated in the next section by analyzing the evolution of the second moment and the PDFs of Y * v . From figures 4. However, the Lewis number is not constant in our simulations due to the decrease of the thermal diffusion coefficient D th = λ ρcp induced by the increase of the gas density and heat capacity. Consequently, the evaporation regime influences the difference 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, which limits the assumption that the Lewis number is constant and equal to one; normally used in combustion models.

Scalar turbulent mixing analysis Evolution of the scalar variance

This section is dedicated to studying scalar mixing using the second moment of the normalized vapor mass fraction. This has been previously done for singlephase ( [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 ([126, 149, 125]) and dense ( [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]) two-phase flows. An exponential decay characterizes the temporal evolution of the scalar variance in single-phase flows. The dissipation rate term in the transport equation for the variance explains this behavior. In DNS of dispersed two-phase flows, Lagrangian techniques are implemented to track the temporal evolution of droplet diameter, velocity, and position. These droplets create vapor pockets as the turbulent velocity field convects them. This method provides good insight into the effects of a scalar source on the temporal evolution of the scalar variance. However, it

L. Germes Martínez

DNS of phase change June 28, 2022 has some limitations, 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 flow of liquid structures are solved, similar to [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. Moreover, temperature, vapor mass fraction, and fluid dynamics are coupled in our formalism through the evaporation rate. Consequently, the vapor mass fraction is not a passive scalar, and the mixing depends on the mass and thermal boundary layers simultaneously.

* v = (Y * v -Y * v ) 2
and its theoretical maximum value

Y * v ,max = (Y * v (1 -Y * v )
of the normalized vapor mass fraction for the two cases studied. In general, the curves in fig. 4.10 have similar behavior. The dispersion increases when the simulation starts until it reaches a maximum value after approximately an eddy turnover time. Then, the dispersion starts to decay to zero, close to the saturation state. This behavior is very different from the results obtained with the single-phase flow but similar to the results obtained with the two-phase dispersed flow 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 variance 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 influential in the system than vaporization, resulting in a decay of the scalar variance. In addition, [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF] observed a similar behavior in dense two-phase flow simulations. In this case, the initial increase in dispersion is caused 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 variance and a sharper decay, resulting in a shorter characteristic time to reach equilibrium.

In a more detailed analysis of the curves presented in fig. 4.10, we observe the influence of the evaporation regime on the temporal evolution of the scalar variance. 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 regime. 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 flows. 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 influence of the Stefan flow 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 rapid decay of the scalar variance in case 1 is also induced by the larger Stefan flows, which push the vapor away from the interface. In addition, the non-homogeneity of the evaporation rate considered in the present work dramatically affects the temporal and spatial evolution of scalar dissipation rate and, consequently, the scalar variance.

Y * v,max (b)
In the following sections, a more detailed analysis of scalar generation and mixing is performed 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 confined two-phase environment with a vaporizing liquid surrounded by dry air, the scalar PDF starts as a Dirac function corresponding to Y * v = 0 when evaporation has not begun. Then, once the equilibrium state is reached, the scalar PDF is another Dirac function at Y * v = 1, where vapor generation is no longer possible due to the absence of vapor gradient.

In this work, we are interested in studying the temporal evolution of the scalar PDF shapes. For this purpose, the PDF of the normalized vapor mass fraction for four successive instants for case 1 is plotted in fig. 4.10. Each PDF corresponds to a specific value of the mean vapor mass fraction: Y * v = 0.18, 0.36, 0.54, 0.72. The general tendency described above is easily seen. An interesting behavior that is observed for the curves presented in 4.11a 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 final value. This behavior has at least two explanations; first, from the table 4.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. Second, there is a wide range of liquid structures with different characteristic lengths, resulting in the different temporal evolution of the temperature profiles of 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 result of solving the temperature and internal flow 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.

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.2 0.4 0.6 0.8 P (Y * v ) (b)
In the analysis of passive scalar mixing in single-phase flows, values of the normalized vapor mass fraction greater than unity are generally not observed. This is because the total amount of vapor is usually present at the initial condition of the simulation. In the simulation performed by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], the vapor fraction of 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 is equal to the imposed vapor mass fraction at the interface.

The first PDF presented in 4.11a correspond to a normalized time of t * 1 = 0.60. At this time, the RMS is almost at its maximum (Y * v = 0.299) which is a 130% of its initial value (see fig. The following PDF corresponds to t * 2 = 1.95 and an RMS of Y * v = 0.256, which means that the decay has already started; thus, the contribution of the scalar dissipation rate is more significant than the contribution of the vaporization source term in the variance transport equation. In this state, the effects of the turbulent velocity field on the scalar mixture are more critical. As a result, the vapor is transported into the gas field away from the interface, resulting in a rightward shift and a decrease in the PDF peak. In addition, the tail of the PDF has been reduced, indicating a reduction in the vapor mass fraction at the interface. Also, the PDF shape exhibits several changes, describing various scenarios, showing a large interface temperature dispersion and a wide range of characteristic lengths of the liquid and gas structures. Subsequently, at the third PDF, which corresponds to a time of t * 3 = 4.3, 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.156). 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 L. Germes Martínez Finally, the shrinking continues in the fourth PDF (t * 4 = 8), consistent with the asymptotic behavior described above. At this point, the RMS has decreased to 43% of its initial value (Y * v = 0.099). 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 fig. 4.11. 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 change of behavior in the last PDF values that change at every instant. This could be interpreted as the result of the mass boundary layer thickness evolution with increasing mean vapor mass fraction. 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 instants than for case 1 (see fig. 4.11) are presented in fig 4 .12. In this case, it can also be found values of Y * v > 1 but with a smaller probability than in case 1. This is because the equilibrium temperature is higher than the initial liquid temperature (see table 4.3). So, a probable reason for having values of Y * v > 1 is that the smaller liquid structures are heated at temperatures higher than the equilibrium temperature. This occurs especially in the first moments of the simulation. The probabilities for the values of Y * v > 1 decrease with time until they become nonexistent in the last PDF.

-4 -2 0 2 4 6 (Y * v -Y * v ) Y * v 0.0 0.2 0.4 0.6 0.8 1.0 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 P (Y * v ) Y * v (b)
For the same values of Y * v , an influence of the evaporation regime in the form of the PDFs can be appreciated. For example, at the first instant, the peak of the PDF has decreased considerably, and the time needed to arrive at the same

Y * v (t * 1 = 7.66
) is almost 13 times bigger than in case 1. For the RMS (Y * v = 0.182), a similar increase with respect to the initial value is observed (125%). The reduced evaporation can explain the time to arrive at the same Y * v , since less vapor is generated at the interface, and the influence of the Stefan flow is smaller. Consequently, the turbulence velocity field has more time to act on the vapor mass fraction field, which explains the reduction of the PDF peak. The same phenomenon occurs in the others PDFs. The second PDF correspond to a time of t * 2 = 10.7 and a Y * v = 0.137 representing 94% its initial value. At this time, the peak of the PDF continues to decrease. In the third and fourth PDFs (t * 3 = 17.51 and t * 4 = 41.65), the width has been reduced considerably. This indicates that a better turbulent mixing is achieved compared to case 1.

Self-similarity

* v (P (Y * v )Y * v with respect to Y * v -Y * v Y * v
) for the same time instants than fig. 4.12 are presented. Similar to the results obtained by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], a significant 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 fig. 4.13, it can be observed that in the first instant of the simulation, the shape of the PDF is affected by the initial conditions. However, the PDF adopts similar shapes once the boundary layers are formed. Nonethe- less, around the peak, we can notice a difference 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 influences the PDF shape. Furthermore, it is observed that the non-homogeneity of the vapor mass fraction at the interface affects the shape and length of the tail on the right side of the PDFs.

The form observed in fig. 4.13 is closer to a log-normal distribution confirming the observations made by [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. Consequently, the PDF shape is not properly represented by a β-distribution, mainly because of a secondary peak at the highest values of the β-PDF form, which is not present in our results.

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. Fig. 4.14 shows the PDFs of the normalized vapor mass fraction at the interface (Y * v s = Yvs Yv,eq ) at t * 1 , t * 2 , t * 3 and t * 4 for both cases. For the interface variables, the asymptotic behavior of the PDF is similar to the gas phase variables. At first, the PDF is spread around the value that corresponds to the initial conditions, Y * v s = 1.24 for the case in fig. 4.14a and Y * v s = 0.5 for the case in fig. 4.14b. Then, the PDF is similar to a Dirac function shape when the system is close to the equilibrium state.

Large dispersion of the vapor mass fraction is observed in the first PDFs. This indicates a wide range of thermal boundary layers thickness since Y * v s depends mainly on the interface temperature. Then, the mean vapor mass fraction and a small peak at the left of the mean value ( ωΓ = 5.6) are presented. For the successive instants, the decrease of the mean evaporation rate and the tail to the right can be observed, causing an increase in the PDF peak. Compared to the interface vapor mass fraction, the evaporation rate mostly depends on the interface temperature but also on the vapor concentration around the interface. These results also depend on the interface geometry, such as the curvature and the liquid and gas structure's characteristic lengths.

Dependence of the evaporation rate on curvature

To better understand the simulation state at the investigated instants, the interface visualization (zero-level set isocontours) colored by the evaporation rate magnitude for the same instants as the PDFs are presented in fig. 4.15a. As for figures 4.2, 4.3 and 4.4, the effect of the turbulent velocity field on the interface are observed. Also, these visualizations help us understand the curvature's influence on the evaporation rate. The upper left figure shows how the evaporation rate is suddenly reduced due to local saturation when the interface is concave. Also, in the upper-right figure, we can see some gas structures enclosed in the liquid and how the evaporation rate is low in these areas due to the same effects. In the last figure, the reduction of liquid mass and evaporation rate compared to the first figure is observed. Figure 4.17 shows 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 configuration 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 fig. 4.17 from left to right and from top to bottom. The first joint PDF (top left of fig. 4.17) 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 two-phase flow simulations in homogeneous shear turbulence. In their simulations, the liquid/gas interface and internal flow of the liquid droplets are resolved, and the evaporation rate depends on the vapor mass fraction and the temperature fields. [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF] observed a dispersion of ω due to the interface deformations and vapor concentration in-homogeneity induced by the turbulent velocity field. Although, only values on the positive side of the curvature axis (representing a convex interface) are presented because the interface deformation is relatively small, and there is no encapsulation of gas structures.

In this work, an opposite convention for the curvature sign is used. It is represented in fig. 4.18, where a deformed liquid structure is shown. Here, 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. 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 bigger 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 because of the difficulty of the vapor to be transported to areas of the gas far from the interface. In the second joint PDF (top right of fig. 4.17), there is an increase in the probability of ω in the negative values of the mean curvature. This indicates an increased saturation in the concave interfaces. Then, in the third joint PDF (bottom left of fig. 4.17), the decrease of ω is most noticeable as a result of the approximation of the simulation to the equilibrium state. Moreover, the values ω are less dispersed because of the homogenization of the vapor mass fraction and temperature fields. Finally, in the last joint PDF (bottom right of fig. 4.17), the probability of both convex and concave interfaces start to become comparable, and the tendency of ω toward zero is more evident.

In fig. 4.19, the temporal evolution of the joint PDF of the evaporation rate for case 2 is presented. In general, the same behavior is observed as in fig. 4.17. However, the effects of the evaporation regime are easily seen. the magnitude of the evaporation regime is considerably smaller, as noted in the previous section (see fig. 4.6). In the earlier moment, the shape of the joint PDFs for both cases are similar despite the difference in the magnitude of the evaporation rate. Once the turbulent mixing begins to homogenize the scalars in the gas phase, the differences in the joint PDFs become more apparent. For the same value of the mean vapor mass fraction, we can observe the flattened shape of the joint PDFs, which indicates a higher saturation level for the same value of the mean vapor mass fraction.

Conclusion

The compressible formalism implemented in chapter 3 is applied to analyze the evaporation process and scalar turbulent mixing in a dense atomization regime. However, minor changes in the formalism are introduced to consider the vapor molecular mass and heat capacity in the computation of density and temperature. The initial condition of the two simulations is chosen based on an estimation of the equilibrium temperature and vapor mass fraction. Here, the initial liquid temperature is varied to consider different evaporation regimes.

Then, the temporal evolution of the evaporation rate in both cases is studied. In case 1, the evaporation rate is larger, and the initial decay is more pronounced. Therefore, the higher values of the evaporation rate explain the liquid temperature decrease in case 1. Moreover, the vapor mass fraction increases more rapidly, resulting in a shorter characteristic time needed to reach equilibrium in case 1.

Furthermore, the turbulent mixing of the vapor mass fraction is analyzed. To this end, the temporal evolution of the variance for both cases is studied. An initial increase on the vapor mass fraction variance is observed. This is explained by vapor production at the interface in a primarily dry air environment. Afterward, the scalar variance decays to zero until it reaches equilibrium. Also, it is observed that the evaporation regime affects this behavior by increasing the maximum value when large evaporation rates are considered, and it is accompanied by a faster decay.

Moreover, the temporal evolution of the PDFs of the vapor mass fraction is studied. It is shown that the PDFs shape strongly depends on the scalar mixing stage. For example, in the initial stage, a peak at the value near zero is observed, accompanied by a long tail covering all possible values of the vapor mass fraction. These results show that the scalar mixing is mainly controlled by the diffusion of the vapor layers near the interface. An interesting behavior is observed in the right tail of the PDFs. Here, we can see the dependence of the shape of the PDFs tails on the mass boundary layer thickness. In addition, the effects of the evaporation regime on the temporal evaluation of scalar PDFs are observed by analyzing the results obtained with case 2. When the evaporation rate is reduced, a longer time is needed to reach the equilibrium value of the mean vapor mass fraction in the system. Consequently, the shapes of the PDFs show higher scalar mixing in case 2.

Finally, the behavior of the evaporation rate and vapor mass fraction at the interface is analyzed using the temporal evolution of the joint PDFs and PDF shapes. Both variables showed a more significant dispersion in the first instants of the simulation. 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.

General conclusion and Prospects

General conclusion

The principal objective of this work is to develop a weakly compressible formalism for the simulation of turbulent two-phase flows with phase change for the study of evaporation and turbulent mixing processes.

The first chapter presents examples of industrial processes where evaporation and turbulent mixing play an essential role. And how a deeper understanding of these phenomena can help optimize their design and reduce their energy consumption and pollutant emotions. Then, the state-of-art of numerical methods dedicated to the DNS of the evaporation is investigated. We first noticed the increased difficulty when heat and mass transfer are considered in two-phase flow simulations, and most publications dedicated to this kind of flow use an incompressible formalism. Important results were obtained that have contributed to a better understanding of the evaporation process, especially in isolated liquid or gas structures. Unfortunately, important phenomena, such as the density and pressure variations due to the vapor generation at the interface can not be captured. Additionally, it is not adapted to simulate an enclosed environment because outflows boundary conditions must be present to maintain the incompressible assumption.

Furthermore, publications that include the compressibility effects in their simulation used a low-mach formulation of the Navier-Stokes equations. One of the characteristics of this strategy is a spatially constant thermodynamic pressure, which could affect the results when several gas structures with different characteristic lengths are present in the domain. This motivates us to develop a weakly-compressible formalism using an all-Mach number formulation of the Navier-Stokes equations to simulate two-phase flows with phase change. Thus, the evaporation and scalar turbulent mixing processes can be studied simultaneously in a closed environment.

The in-house code ARCHER is adapted to simulate incompressible two-phase flows with phase change to achieve our objective. To this end, the mass conser-vative interface capturing method (CLSVOF) is modified to consider the liquid mass variation due to vaporization, using a sink term in the VOF transport equation. The sink term is estimated using the evaporation rate and the local surface density; computed using the PLIC reconstruction. In addition, for the treatment of the velocity discontinuities at the interface, three methods proposed in the literature ( [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Lc Malan | A geometric vof method for interface resolved phase change and conservative thermal energy advection[END_REF]) are tested and compared in a 2D static cylinder configuration with a constant evaporation rate. We choose for our formalism the projection method based on two velocities and explicit implementation of the velocity jump proposed by [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. Then, the thermodynamics of the formalism is validated with the results obtained with the D 2 law. Moreover, the length of the recirculation zone behind a convected cylinder with and without evaporation is compared with the reference results obtained with a solid cylinder. In the case with no evaporation, good agreement was obtained with the reference, and a detachment and enlargement of the recirculation zone were observed when evaporation is considered. These results show the effects of the Stefan flow on the system's dynamics.

The next step is to extend the formalism into a weakly-compressible formalism using the pressure-based method proposed by [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF]. Additionally, the CLSVOF method for the interface representation and the projection method was modified to account for the liquid evaporation and compressibility effects, such as thermal dilatation. To this end, a pressure equation is derived and implemented. The total mass conservation properties of the compressible formalism are validated by studying the temporal evolution of the liquid and gas mass in a 2D static cylinder with a constant evaporation rate. Then, a 3D compressible twophase HIT configuration with phase change is presented to show the formalism's ability to handle atomization regimes.

In the following chapter, our weakly-compressible formalism is implemented to analyze the evaporation and scalar turbulent mixing processes. For this part, small changes in the formalism have been introduced to consider the vapor molecular mass and heat capacity in the computation of density and temperature. Also, two simulations are performed with a compressible two-phase HIT configuration with phase change, varying the initial liquid temperature to have different evaporation regimes. By estimating the temperature and the vapor mass fraction at the equilibrium state, the effects of the evaporation regime are observed. In the case of a large evaporation rate, the equilibrium temperature is smaller than the initial liquid temperature. The opposite behavior is observed in the case of a low evaporation rate. This is attributed to the amount of energy consumed by the liquid during evaporation.

The evaporation process is studied by comparing the temporal evolution of the liquid and gas temperature, evaporation rate, and vapor mass fraction. First, as expected, we observed the liquid and gas temperature tends toward its equilibrium value, but the behavior is significantly affected by the evaporation rate. Then, the scalar turbulent mixing process is analyzed via the temporal evolution of the variance and PDF shapes of the vapor mass fraction. This has allowed us to observe the different states of the turbulent mixing. From the initial stage, where mixing is mainly controlled by the diffusion of the vapor layer close to the interface and its evolution toward a shape similar to a Dirac function at the equilibrium. Here, significantly different results have been obtained compared with the single-phase or dispersed two-phase configurations. We conclude that the effect of the interface and the evaporation on the scalar mixing are not negligible. Then, the normalized PDFs of the vapor mass fraction are presented, observing that the PDFs evolve to a similar shape, close to a log-normal distribution, in the later instants of the simulation. Finally, the relation between the evaporation rate and the interface curvature is analyzed using the temporal evolution of the joint PDFs. Here, a larger evaporation rate is observed in the positive curvatures, which corresponds to the concave section of the interface.

Prospects

In this work, promising results have been obtained to contribute to a deeper understanding of the evaporation and scalar turbulent mixing processes in a dense atomization regime. Considering the phenomena captured by our formalism, e.g., non-homogeneous evaporation rate and temperature at the interface, different gas structures with independent thermodynamic pressure, a gas density that depends on the local temperature and pressure, etc., the number of possible applications is vast. Although, several numerical development prospects can be investigated to improve the incompressible and compressible formalisms. This section will explore the possible prospects of this work related to numerical developments and physical applications.

Physical applications

Beginning with the incompressible formalism, a falling droplet configuration with phase change could help access the range of applicability of the classical models found in the literature. For example, a continuation of the work done by [START_REF] Duret | Evaporating droplets: comparisons between dns and modelling[END_REF],

where the effects of the Stefan flow, latent heat, and interface regression in the heat and mass transfer of a droplet are investigated, could result in more realistic correlations for the Nusselt and Sherwood numbers used in evaporation modeling. Another possible application for our formalism is the analysis of the particles captured by droplets, generally used for studying air pollution. Significant improvement could be made in this field since the models do not consider the Stefan flows, which should influence the particle capture because it pushes the particles away from the interface when evaporation occurs.

For the compressible formalism, a more detailed study conducted on the two-phase HIT configuration could be used to propose new subgrid-scale models for LES and RANS simulations. In addition, the study of the temporal evolution of the vapor mass fraction PDF showed that the form of the PDF in the last moments of the simulation is not well represented by the Gaussian or β distributions, contrary to the results obtained in single-phase or dispersed two-phase simulations. This is particularly interesting for turbulent combustion simulations, where the assumed PDF approach is commonly used to characterize the mixing and estimate the chemical reaction rate.

Another possible application for our compressible formalism is the study of 'boil off' in cryogenic fluid containers. This is important in storing and transporting liquefied fuels, such as natural gas and hydrogen. The 'boil off' consists of spontaneous vaporization when the liquid fuel is transported in a tank. Consequently, the pressure increases, and a safety valve releases some of the fuel. In the case of hydrogen, this release constitutes a significant percentage of the fuel loss.

Numerical development prospects

There are several possible improvements to our formalism, including implementing a consistent computation of mass and momentum fluxes in two-phase flows simulation. A candidate is the numerical scheme implemented in ARCHER by [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF]. But this scheme has to be modified to account for velocity discontinuities generated by the phase change. In addition, an essential issue of the simulation presented in chapter 4 is the related computational cost. Consequently, new strategies are needed to reduce the number of interface extrapolation and Poisson solvers required for the simulations. In addition, for some configurations, an adaptive mesh refinement (AMR) method could help capture boundary layers while reducing the number of mesh points in the whole domain.

Currently, the simulations with the compressible formalism are limited to L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022

an enclosed environment because a proper outflow boundary condition is not implemented yet. This is challenging because of the treatment of acoustic waves. The advantage of proper boundary conditions is that the range of applications of our formalism can be extended to liquid jets or liquid films. Moreover, new numerical strategies to consider evaporation and boiling simultaneously can significantly improve our formalism. This has already been done in the literature with a level-set method in [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF].

Annexe Synthèse générale

Introduction

Le transfert de chaleur et de masse dans les écoulements multiphasiques joue un rôle essentiel dans de nombreux processus environnementaux et industriels, tels que la formation de nuages et de gouttes de pluie, le refroidissement par atomisation, les réacteurs nucléaires, le séchage par atomisation, la combustion de carburant, etc. Ce dernier exemple est présent dans une grande diversité de moteurs à combustion interne, par exemple, le moteur d'une voiture ou un moteur pour la production d'énergie électrique et un moteur à combustion continue d'un avion ou d'une rocket.

Ces processus ont un principe commun : le carburant est injecté sous forme liquide dans la chambre de combustion. Il subit ensuite plusieurs changements morphologiques, de la formation initiale du jet à la rupture secondaire pour former de petites structures liquides. Ces structures interagissent avec un environnement gazeux à haute température où la température du liquide augmente et où l'évaporation commence. L'étape suivante est la combustion de la vapeur de carburant.

À ce stade, l'évaporation dépend fortement de l'interaction entre les structures liquides et le mélange turbulent de la vapeur de carburant. De plus, l'efficacité de la combustion est liée au mélange de la vapeur de carburant et de l'oxydant. Par conséquent, la compréhension des processus d'évaporation et de mélange turbulent pourrait aider à optimiser la conception des chambres de combustion et ainsi réduire la consommation de carburant et les émissions polluantes.

Le séchage par atomisation est un autre procédé industriel dans lequel l'évaporation et le mélange turbulent jouent un rôle essentiel. Il est largement utilisé dans l'industrie alimentaire, par exemple pour la production de lait en poudre ou de café instantané. Le séchage par atomisation est un procédé de séchage permettant de transformer un liquide en un produit séché particulaire en une seule opération [START_REF] Robert H Perry | Perry's chemical engineers' handbook[END_REF]. Il augmente la surface du liquide par atomisation et est suivi d'un contact avec de l'air à haute température. À ce stade, le liquide s'évapore et des particules se forment alors que la gouttelette est encore en suspension ; l'étape finale est la récupération du produit séché. Comme pour la combustion par pulvérisation, l'efficacité de l'évaporation du liquide dépend du mélange de la vapeur et de l'air sec et de la taille moyenne des gouttelettes.

Ce travail se concentrera sur la simulation numérique directe (DNS) du processus d'évaporation sur des régimes d'atomisation turbulents. Dans ces simulations, le changement de phase se produit à l'interface et à des températures inférieures au point d'ébullition. Il est principalement contrôlé par les gradients de densité du gaz, de température et de fraction massique de vapeur à l'interface. La vapeur générée à l'interface crée une divergence de vitesse qui dépend du rapport de densité entre les phases. Par conséquent, un flux de Stefan est produit, qui est un flux de masse dans la direction normale de l'interface.

Même si de nombreux codes dans la littérature sont dédiés à l'évaporation avec un formalism DNS, la plupart d'entre eux se concentrent sur des gouttelettes ou des bulles isolées et utilisent un formalisme incompressible. Par conséquent, les simulations sont limitées à des environnements ouverts, et les changements de densité et de pression dus à la production de vapeur et au mélange ne peuvent pas être pris en compte. Ce dernier point est important pour estimer l'influence du flux de Stefan sur la dynamique du système.

Un formalisme compressible doit être mis en oeuvre pour une représentation plus réaliste du processus d'évaporation dans un régime d'atomisation. Il existe plusieurs voies pour surmonter cette tâche, par exemple, les solveurs de Riemann et les solveurs basés sur la pression. Cependant, dans les simulations avec un faible nombre de Mach, la précision des solveurs de Riemann a tendance à diminuer. Par conséquent, ce travail se concentre sur un solveur basé sur la pression couplé à une méthode de capture d'interface. Cette méthode présente plusieurs avantages, notamment sa compatibilité avec les méthodes de projection standard généralement utilisées dans le formalisme incompressible. De plus, une formulation implicite des termes acoustiques est permise, ce qui réduit les restrictions sur le pas de temps. Comme la méthode est basée sur un formalisme entièrement compressible pour les écoulements multiphasiques, les changements de pression spatiaux et temporels sont capturés. Cela permet aux simulations de traiter de nombreuses structures liquides et gazeuses, chacune ayant sa propre densité et pression thermodynamique. Après la vaporisation, la vapeur est transportée et mélangée avec la phase gazeuse. Dans la plupart des applications industrielles, la phase gazeuse est soumise à un champ de vitesse turbulent qui affecte de manière significative le mélange vapeur-air. Ce mélange turbulent peut être défini comme un processus en trois étapes : entraînement, agitation et mélange ( [START_REF] Paul E Dimotakis | Turbulent mixing[END_REF]). La première étape correspond à la génération de vapeur, et le mélange est contrôlé par la diffusion des couches de vapeur proches de l'interface vers la phase gazeuse. Dans la deuxième étape, le mélange est régi principalement par le champ de vitesse turbulent et est caractérisé par la déformation de l'interface due à l'interaction avec les plus grands tourbillons de l'écoulement. Enfin, lorsque les fluides deviennent plus homogènes dans l'état final, les gradients de concentration diminuent jusqu'à ce qu'il n'y ait plus de transfert de chaleur ou de masse. La plupart des recherches consacrées à ce phénomène se déroulent dans un environnement purement gazeux ou dans des écoulements diphasiques dispersés, laissant de côté l'influence des interactions vapeur/interface. Ceci nous motive à utiliser notre formalisme compressible pour contribuer à une meilleure compréhension du mélange turbulent.

Le travail de cette thèse de doctorat est basé sur le code interne ARCHER pour les écoulements diphasiques incompressibles et compressibles 2D et 3D ( [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]). Il a été développé dans le laboratoire Coria au cours des deux dernières décennies, et il a été validé et appliqué dans plusieurs configurations liées au processus d'atomisation.

Simulations d'écoulements diphasiques incompressibles avec changement de phase

Dans ce chapitre, un système composé de deux fluides non miscibles (par exemple, un liquide monocomposant et un mélange de gaz) est étudié. Tout d'abord, les équations résolues dans notre formalisme incompressible dédié au DNS des écoulements diphasiques avec changement de phase sont présentées. Puisque le formalisme one-fluid est adopté pour représenter chaque phase, une atten-tion particulière doit être accordée aux conditions de saut à l'interface. Pour cette raison, la partie suivante du chapitre est une discussion sur les conditions de saut standard utilisées pour les simulations d'écoulements diphasiques inertes et les conditions de saut supplémentaires requises pour le changement de phase. Cette partie est suivie par les méthodes numériques mises en oeuvre dans le code ARCHER interne : [START_REF] Ménard | Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary breakup of a liquid jet[END_REF]. Cette partie traite des modifications apportées à la méthode de capture d'interface pour tenir compte du liquide évaporé. En outre, trois méthodes différentes pour gérer le saut de vitesse dans la méthode de projection sont présentées. Nous montrons ensuite des résultats dans une configuration de cylindre statique 2D pour valider et comparer les trois méthodes. Enfin, des résultats d'application sont présentés pour montrer la précision et la robustesse de notre formalisme incompressible.

Résultats

Régression de l'interface due à un taux de vaporisation constant

• Valider l'estimation du terme puit dans l'équation VOF (eq. 2.24) responsable de la régression de l'interface due au changement de phase.

• Tester l'efficacité des méthodes responsables de la génération de la vitesse du liquide mentionnées dans la section 2.3.2. Cela nous permet de choisir la méthode la plus appropriée pour ce travail. À cette fin, un cylindre d'eau statique 2D avec un taux d'évaporation constant est analysé. Puisque le taux d'évaporation est imposé, les équations de conservation de l'énergie et des espèces ne sont pas considérées pour cette configuration. Cependant, les équations de Navier-Stokes sont résolues.

D 2 law
Une configuration statique de gouttelettes en 3D est comparée à la loi D 2 pour valider l'implémentation de la thermodynamique dans le formalisme incompressible. Selon cette loi, le carré du diamètre de la gouttelette diminue linéairement avec le temps en raison de la diffusion de la masse vers le gas [START_REF] William | Fluid dynamics and transport of droplets and sprays[END_REF]). Pour obtenir l'équation du rayon carré de la goutte, il faut résoudre l'équation de continuité en coordonnées sphériques en supposant l'incompressibilité et un débit massique liquide constant de l'interface vers l'air.

Sur la Fig. 4.20, la température (4.20a) et la fraction massique de vapeur (4.20b) sont présentées à t = 0.08 s pour la maille la plus fine (64 3 ). Une propriété de symétrie sphérique autour de la gouttelette est observée pour les deux champs de temperature et le champ de fraction massique de vapeur. En raison de la différence relativement importante entre la température du liquide et celle Ce chapitre présente les équations résolues pour les simulations d'écoulements diphasiques compressibles avec transfert de chaleur et de masse. De plus, les schémas numériques mis en oeuvre pour les résoudre sont discutés. Comme dans le chapitre 2, le formalisme one-fluid est considéré, où un seul système d'équations est résolu pour toutes les phases du domaine. Nous utilisons également des configurations dans lesquelles un liquide à un seul composant s'évapore dans un mélange de gaz idéaux. Dans ce cas, les deux phases sont considérées comme compressibles, et une équation d'état adaptée à chaque phase est choisie pour fermer notre système d'équations. Une méthode conservative en masse pour capturer l'interface (CLSVOF) est utilisée. Cette méthode a été adaptée pour tenir compte de l'évolution de la masse liquide due au changement de phase et à la variation de densité. Ensuite, la méthode de capture d'interface est couplée à une méthode de projection modifiée pour l'évolution temporelle de l'équation de quantité de mouvement basée sur deux vitesses.

L. Germes Martínez

Un cas de validation d'un cylindre qui s'evapore dans un environnement fermé est présenté pour le formalisme compressible. Enfin, une configuration 3D de Turbulent Homogène Isotrope (THI) diphasique est étudiée pour illustrer la capacité de la méthode à traiter des écoulements diphasiques turbulents compressibles avec changement de phase.

Résultats

Configuration THI diphasique compressible avec changement de phase Similaire aux travaux précédents, l'idée est d'étudier l'influence du phénomène d'évaporation à l'intérieur d'un écoulement diphasique turbulent. Dans [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], un scalaire passif a été utilisé pour représenter le processus d'évaporation et de mélange. Parmi les limitations de cette procédure, l'influence du flux de Stefan dans le champ de vitesse et l'augmentation de pression due au processus d'évaporation n'ont pas été considérés. Une amélioration significative a été proposée par [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF] où un terme source de changement de phase a été directement introduit dans les équations de continuité et de pression. Les résultats obtenus ont illustré de bonnes propriétés de conservation de la masse. Cependant, l'influence de la température et la fraction de masse de vapeur à l'interface n'était pas encore implémentée dans le formalisme général : le taux de vaporisation était directement imposé. Contrairement à [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF], ici, le saut de vitesse est introduit implicitement dans l'équation de Poisson pour la pression dans les points ghost du champ de vitesse du gaz. Le couplage entre les équations d'énergie, de continuité et de quantité de mouvement est ajouté en supposant l'équilibre thermodynamique à l'interface.

La Fig. 4.22 montre un champ de température instantané et l'iso-contour de la fonction level-set, qui représente l'interface liquide/gaz. Comme dans d'autres configurations THI compressibles et incompressibles présentées dans les travaux précédents, de nombreux événements de rupture et de coalescence sont observés. De plus, de petites structures gazeuses sont observées à l'intérieur du liquide ; ces structures sont caractérisées par une augmentation de la pression due à la création de vapeur et une diminution de la température en raison de la consommation d'énergie induite par le changement de phase. De plus, l'influence de la convection due au champ de vitesse turbulent dans le champ de température est observée. Des phénomènes similaires se retrouvent dans le champ de fractionnement de la masse de vapeur (Fig.

Vers une description complète des processus de mélange et d'évaporation

Dans la littérature, le mélange turbulent scalaire est étudié dans plusieurs scénarios, tels que l'écoulement monophasique, l'écoulement diphasique dispersé et dense [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF][START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF][START_REF] Reveillon | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF][START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. Dans le premier exemple, l'évolution temporelle de la variance scalaire est caractérisée par une décroissance exponentielle due au taux de dissipation scalaire [START_REF] Eswaran | Direct numerical simulations of the turbulent mixing of a passive scalar[END_REF] ; lorsqu'un gradient scalaire moyen n'est pas imposé. Ce comportement change radicalement lorsqu'une configuration d'écoulement diphasique avec changement de phase est considérée. Dans ce cas, une augmentation initiale de la variance scalaire est observée jusqu'à atteindre une valeur maximale, suivie d'une décroissance quasi-exponentielle. [START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF] a expliqué ce comportement en analysant l'équation de transport de la variance scalaire, où différents termes sources apparaissent lorsqu'il y a un liquide qui se vaporise dans le domaine. Ici, l'évolution temporelle de la variance scalaire résulte d'un équilibre entre le taux de dissipation scalaire et le terme source de vaporisation. Dans le travail de [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], le mélange turbulent scalaire dans une configuration THI diphasique avec une fraction massique de vapeur constante à l'interface est étudié. Dans ce chapitre, une continuation de cette étude est effectuée en utilisant notre formalisme compressible, où le couplage entre la fraction de masse de vapeur et la température est résolu. Ainsi, une représentation nonhomogène de la fraction de masse de vapeur à l'interface est possible. Tout d'abord, plusieurs améliorations aux équations résolues et à la méthode numérique décrite dans 3 sont présentées et discutées. Par exemple, l'influence de la masse molaire de la vapeur est prise en compte dans le calcul de la densité et de la capacité thermique du gaz. Ensuite, les méthodes numériques consacrées à la résolution des équations avec les conditions de saut d'interface correspondantes et la manière de calculer les quantités moyennées en surface sont présentées.

Ici, le mélange scalaire et l'évaporation dans un environnement d'atomisation sont analysés en utilisant la configuration THI diphasique compressible. À cette fin, nous effectuons deux simulations avec différents régimes d'évaporation. De plus, les effets du régime d'évaporation sur l'évolution temporelle de la fraction massique moyenne de la vapeur et des températures moyennes du liquide et du gaz sont discutés. Ensuite, nous analysons le mélange turbulent scalaire en considérant l'évolution temporelle des principales statistiques de la fraction de masse de vapeur. Enfin, nous étudions le comportement du taux d'évaporation et de la fraction massique de vapeur à l'interface.

Résultats Transfert de chaleur et de masse

Les champs de température instantanés aux frontières du domaine pour les deux cas avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 sont présentés dans la fig. 4.24. À ce stade des simulations, les effets du champ de vitesse turbulent sont facilement appréciés dans le transport de la température. Comme dans le chapitre 3, de nombreuses interactions entre les structures liquide et gazeuse induites par le champ de vitesse sont observées. L'opacité de la visualisation de l'interface est réduite pour montrer les bulles enfermées dans les structures liquides plus grandes. En raison du volume de gaz réduit, ces bulles sont caractérisées par l'atteinte d'états d'équilibre à des temps plus courts que le reste du domaine. L'influence du volume de gaz sur le temps nécessaire pour atteindre l'équilibre est étudiée plus en détail par [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. En outre, les températures les plus faibles sont trouvées près de l'interface, où une large gamme d'épaisseur de la couche limite thermique est observée. Ces variations d'épaisseur dépendent de la déformation de l'interface et de la température du gaz environnant. Par exemple, une structure liquide se déplaçant vers un gaz à haute température présentera une épaisseur de couche limite thermique plus faible à l'avant en raison des gradients de température plus importants.

En comparant les figures 4.24, on observe une diminution de la température maximale dans le cas 2. De plus, il y a une plus grande homogénéisation de la température du gaz loin de l'interface dans la figure 4.24b. Ceci indique un niveau de mélange plus élevé pour la même valeur de Y * v dans le cas 2 par rapport au cas 1.

Des phénomènes similaires sont observés dans la fig. 4.25 pour les champs de fraction massique de vapeur instantanée aux frontières du domaine à t 2 . Par rapport aux champs de température, de plus grands gradients sont présents dans la phase gazeuse loin de l'interface. Ceci est dû au fait que la diffusivité de la masse (D m = 3, 82 × 10 -7 m 2 • s -1 ) est plus petit que la diffusivité thermique (D th = 1, 70 × 10 -6 m 2 • s -1 ), ce qui entraîne une plus grande influence du terme de convection dans l'équation de transport de la fraction de masse de vapeur par rapport au terme de diffusion. De plus, puisque la fraction de masse de vapeur à l'interface dépend de la température, le taux d'évaporation est également inhomogène. Récemment, des observations similaires ont été faites pour une gouttelette d'évaporation unique dans un environnement turbulent homogène et isotrope ( [START_REF] Michael S Dodd | Analysis of droplet evaporation in isotropic turbulence through dropletresolved dns[END_REF]) et pour des gouttelettes d'évaporation dispersées dans un environnement de cisaillement turbulent homogène ( [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF]). Dans [START_REF] Scapin | Finite-size evaporating droplets in weakly-compressible homogeneous shear turbulence[END_REF], les auteurs L. itérations initiales de diffusion. Ensuite, après l'activation de l'évaporation, une chute soudaine de la température est observée avant de commencer à se stabiliser après deux temps de retournement des tourbillons. Ceci est dû au fait que la magnitude du taux d'évaporation atteint une valeur maximale au début de la simulation. En effet, la couche limite de masse n'est pas encore complètement formée, ce qui provoque une augmentation du gradient de fraction de masse de vapeur à l'interface. Par la suite, la vapeur générée à l'interface augmente le niveau de saturation, diminuant le taux d'évaporation et ralentissant la diminution de la température. A la fin de la simulation, on observe une pente négative dans l'évolution de la température du liquide, ce qui signifie que l'équilibre thermique n'a pas encore été atteint.

Un comportement opposé est observé pour la température du liquide dans le cas 2. La figure 4.27b montre une augmentation de la température du liquide causée par le gradient de température initial entre la température du liquide et celle du gaz. Comme dans le cas 1, les températures du liquide et de l'interface ont des valeurs différentes en raison des sous-itérations initiales. Ensuite, l'augmentation de la température ralentit lorsque le flux de chaleur du gaz vers le liquide diminue. Il en résulte une stabilisation de la température, représentée par un plateau après environ 15 temps de rotation des tourbillons. En comparant les deux images de la fig. Analyse du mélange turbulent scalaire Évolution de la variance scalaire Cette section est consacrée à l'étude du mélange scalaire en utilisant le second moment de la fraction massique normalisée de la vapeur. Ceci a été fait précédemment pour des écoulements monophasique ( [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]), diphasique dispersés ( [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 | Effects of the preferential segregation of droplets on evaporation and turbulent mixing[END_REF]) et denses ( [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]). Une décroissance exponentielle caractérise l'évolution temporelle de la variance scalaire dans les écoulements monophasés. Le terme de taux de dissipation dans l'équation de transport de la variance explique ce comportement. Dans les DNS d'écoulements diphasiques dispersés, des techniques lagrangiennes sont mises en oeuvre pour suivre l'évolution temporelle du diamètre, de la vitesse et de la position des gouttelettes. Ces gouttelettes créent des poches de vapeur lorsque le champ de vitesse turbulent les convecte. Cette méthode donne un bon aperçu des effets d'une source scalaire sur l'évolution temporelle de la variance scalaire. Cependant, elle a quelques limitations, par exemple, les interactions entre le scalaire et l'interface liquide/gaz ne sont pas résolues, et un modèle idéalisé pour la vaporisation du liquide est utilisé.

Dans ce travail, l'interface et l'écoulement interne des structures liquides sont résolus, de manière similaire à [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF]. De plus, la température, la fraction massique de vapeur et la dynamique des fluides sont couplées dans notre formalisme par le biais du taux d'évaporation. Par conséquent, la fraction de masse de vapeur n'est pas un scalaire passif, et le mélange dépend simultanément des couches limites massique et thermique. ) de la fraction massique de vapeur normalisée pour les deux cas étudiés. En général, les courbes de la fig. 4.31 ont un comportement similaire. La dispersion augmente lorsque la simulation commence jusqu'à ce qu'elle atteigne une valeur maximale après environ un temps de retournement des tourbillons. Ensuite, la dispersion commence à décroître jusqu'à zéro, proche de l'état de saturation. Ce comportement est très différent des résultats obtenus avec l'écoulement monophasique mais similaire aux résultats obtenus avec les simulations d'écoulement dispersé diphasique. Comme discuté par [START_REF] Reveillon | Spray vaporization in nonpremixed turbulent combustion modeling: a single droplet model[END_REF], l'augmentation initiale de la variance scalaire est causée par les termes sources apparaissant dans son équation de transport lorsqu'un liquide s'évaporant génère de la vapeur. La valeur maximale de la RMS est atteinte lorsque le taux de dissipation scalaire équilibre la génération de dispersion. Après avoir atteint le maximum, le mélange turbulent devient plus influent dans le système que la vaporisation, ce qui entraîne une décroissance de la variance scalaire. En outre, [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF] a observé un comportement similaire dans des simulations d'écoulement diphasique dense. Dans ce cas, l'augmentation initiale de la dispersion est due à la génération de vapeur à l'interface et à la diversité de l'épaisseur de la couche limite formée autour des structures liquides. En outre, la fraction volumique initiale du liquide sur l'évolution temporelle de la variance scalaire est étudiée. On observe qu'une diminution du volume de gaz disponible entraîne une augmentation de la valeur maximale de la variance scalaire et une décroissance plus nette, ce qui se traduit par un temps caractéristique plus court pour atteindre l'équilibre.

Dans les sections suivantes, une analyse plus détaillée de la génération et du mélange de scalaires est effectuée en calculant les PDF scalaires et les PDF conjointes.

Evolution des PDFs scalaires

Le comportement asymptotique du mélange scalaire a été largement étudié et est bien connu dans la littérature. Pour un environnement diphasique confiné avec un liquide qui se vaporise entouré d'air sec, la PDF scalaire commence comme une fonction de Dirac correspondant à Y * v = 0 lorsque l'évaporation n'a pas commencé. Puis, une fois l'état d'équilibre atteint, la PDF scalaire est une autre fonction de Dirac à Y * v = 1, où la génération de vapeur n'est plus possible en raison de l'absence de gradient de vapeur. Dans ce travail, nous sommes intéressés par l'étude de l'évolution temporelle des formes de la PDF scalaire. Dans ce but, la PDF de la fraction massique de vapeur normalisée pour quatre instants successifs pour le cas 1 est tracée dans la fig. 4.31. Chaque PDF correspond à une valeur spécifique de la fraction massique de vapeur moyenne : Y * v = 0, 18, 0, 36, 0, 54, 0, 72. La tendance générale décrite ci-dessus est facilement observable. Un comportement intéressant qui est observé pour les courbes présentées dans 4.32a est la présence de valeurs de Y * v > 1. Par conséquent, la fraction massique de vapeur n'est pas limitée à 1 lorsqu'elle est normalisée par sa valeur à l'état de saturation, même s'il s'agit de sa valeur finale. Ce comportement a au moins deux explications ; premièrement, à partir du tableau 4.3, on peut observer que la température à l'équilibre est inférieure à la température du liquide. Par conséquent, la fraction massique de vapeur à l'interface d'équilibre est plus petite que la fraction massique de vapeur initiale. Deuxièmement, il existe une large gamme de structures liquides avec différentes longueurs caractéristiques, ce qui entraîne une évolution temporelle différente des profils de température de chaque structure liquide. Cela signifie que certaines petites structures liquides n'ont pas le même comportement que le système global. Par exemple, une petite gouttelette convexe en contact avec un gaz à haute température va s'évaporer et augmenter en température simultanément, même si la température moyenne du système tend à diminuer. Ceci est également le résultat de la résolution de la température et du flux interne de la structure liquide. De plus, ces résultats démontrent le couplage fort entre les équations de conservation de l'énergie et des espèces et la variété des phénomènes capturés par notre formalisme compressible. Dans l'analyse du mélange scalaire passif dans les écoulements monophasiques, des valeurs de la fraction massique normalisée de vapeur supérieures à l'unité ne sont généralement pas observées. Ceci est dû au fait que la quantité totale de vapeur est habituellement présente à la condition initiale de la simulation. Dans la simulation réalisée par [START_REF] Duret | Dns analysis of turbulent mixing in two-phase flows[END_REF], la fraction de vapeur de l'interface est constante et imposée lors de l'initialisation de la simulation. Par conséquent, l'évolution temporelle de la température et sa relation avec la fraction massique de vapeur ne sont pas résolues. En conséquence, la fraction massique de vapeur à l'équilibre est égale à la fraction massique de vapeur imposée à l'interface.

Les premières PDF présentées dans 4.32a correspondent à un temps normalisé de t * 1 = 0, 60. À ce moment, le RMS est presque à son maximum (Y * v = 0, 299) qui est de 130% de sa valeur initiale (voir fig. 4.10a). La PDF présente un pic prononcé à Y * v = 0 correspondant à l'air sec, et une longue queue con-tenant toutes les valeurs possibles de la fraction massique de vapeur normalisée.

Ceci correspond aux premiers instants après le développement de la couche limite de masse, où le mélange est principalement gouverné par la diffusion des couches de vapeur autour de l'interface. Cette étape du mélange scalaire est caractérisée par de grands gradients de température et de fraction massique de vapeur à l'interface, expliquant le taux de changement observé dans l'évolution temporelle de la T * l , T * g , Y * v et ωΓ (voir figures 4. 27, 4.29, 4.30, et 4.28).

La PDF suivante correspond à t * 2 = 1, 95 et à une RMS de Y * v = 0, 256, ce qui signifie que la décroissance a déjà commencé ; ainsi, la contribution du taux de dissipation scalaire est plus importante que la contribution du terme source de vaporisation dans l'équation de transport de la variance. Dans cet état, les effets du champ de vitesse turbulent sur le mélange scalaire sont plus critiques. En conséquence, la vapeur est transportée dans le champ de gaz loin de l'interface, ce qui entraîne un déplacement vers la droite et une diminution du pic de la PDF. En outre, la queue de la PDF a été réduite, indiquant une réduction de la fraction de masse de vapeur à l'interface. De plus, la forme de la PDF présente plusieurs changements, décrivant divers scénarios, montrant une grande dispersion de la température de l'interface et une large gamme de longueurs caractéristiques des structures liquides et gazeuses.

Par la suite, à la troisième PDF, qui correspond à un temps de t * 3 = 4.3, un decalage du pic et un rétrécissement de la longeur de la PDF sont observés, montrant une homogénéisation de la fraction de masse de vapeur dans la phase gazeuse loin de l'interface, qui est également caractérisée par une diminution de la RMS (Y * v = 0.156). De plus, la probabilité de valeurs inférieures à Y * v = 0.25 est proche de zéro, ce qui signifie que la vapeur est répartie sur l'ensemble du volume gazeux. Une queue est également présente sur le côté droit de la PDF, correspondant aux plus grandes valeurs près de l'interface.

Enfin, le rétrécissement se poursuit dans la quatrième PDF (t * 4 = 8), conformément au comportement asymptotique décrit ci-dessus. À ce stade, le RMS a diminué à 43% de sa valeur initiale (Y * v = 0.099). Le pic de la PDF est plus grand, et sa forme présente moins d'irrégularités, ce qui indique un niveau plus élevé d'homogénéisation de la fraction massique de la vapeur. dernières valeurs de PDF qui changent à chaque instant. Cela pourrait être interprété comme le résultat de l'évolution de l'épaisseur de la couche limite de masse avec l'augmentation de la fraction de masse de vapeur moyenne. Il convient de mentionner que ces résultats décrivent un comportement général du système et que l'épaisseur de la couche limite de masse varie en fonction des conditions thermodynamiques locales dans chaque structure et chaque cellule proche de l'interface.

Les PDF de la fraction massique de vapeur pour le cas 2 aux mêmes instants que pour le cas 1 sont présentées dans la fig. 4

Conclusion générale

Ce travail a pour objectif principal de développer un formalisme faiblement compressible pour la simulation d'écoulements diphasiques turbulents avec changement de phase pour l'étude des processus d'évaporation et de mélange turbulent.

Le premier chapitre présente des exemples de procédés industriels où l'évaporation et le mélange turbulent jouent un rôle essentiel. Et comment une compréhension plus approfondie de ces phénomènes peut aider à optimiser leur conception et à réduire leur consommation d'énergie et leurs émissions polluantes. Ensuite, l'état de l'art des méthodes numériques dédiées au DNS de l'évaporation est étudié. Nous avons d'abord remarqué la difficulté accrue lorsque les transferts de chaleur et de masse sont considérés dans les simulations d'écoulements diphasiques, et la plupart des publications dédiées à ce type d'écoulement utilisent un formalisme incompressible. Nous avons obtenu des résultats importants qui ont contribué à une meilleure compréhension du processus d'évaporation, en particulier dans des structures liquides ou gazeuses isolées. Malheureusement, des phénomènes importants, tels que les variations de densité et de pression dues à la génération de vapeur à l'interface, ne peuvent être capturés.

En outre, les publications qui incluent les effets de compressibilité dans leur simulation ont utilisé une formulation à faible nombre de Mach des équations de Navier-Stokes. L'une des caractéristiques de cette stratégie est une pression thermodynamique spatialement constante, ce qui pourrait affecter les résultats lorsque plusieurs structures de gaz avec des longueurs caractéristiques différentes sont présentes dans le domaine. Ceci nous motive à développer un formalisme faiblement compressible utilisant une formulation complètement compressible des équations de Navier-Stokes pour simuler des écoulements diphasiques avec changement de phase. Ainsi, les processus d'évaporation et de mélange turbulent scalaire peuvent être étudiés simultanément dans un environnement fermé.

Le code interne ARCHER est adapté pour simuler des écoulements diphasiques incompressibles avec changement de phase pour atteindre notre objectif. À cette fin, la méthode de capture d'interface conservative en masse (CLSVOF) est modifiée pour prendre en compte la variation de masse liquide due à la vaporisation, en utilisant un terme puit dans l'équation de transport VOF. Le terme 'sink' est estimé à l'aide du taux d'évaporation et de la densité de surface locale, calculés à l'aide de la reconstruction PLIC. De plus, pour le traitement des discontinuités de vitesse à l'interface, trois méthodes proposées dans la littérature ( [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Lc Malan | A geometric vof method for interface resolved phase change and conservative thermal energy advection[END_REF]) sont testées et comparées dans une configuration cylindrique statique 2D avec un taux d'évaporation constant. Nous choisissons pour notre formalisme la méthode de projection basée sur deux vitesses et l'implémentation explicite du saut de vitesses proposé par [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. Ensuite, la thermodynamique du formalisme est validée avec les résultats obtenus avec la loi D 2 . De plus, la longueur de la zone de recirculation derrière un cylindre convexe avec et sans évaporation est comparée aux résultats de référence obtenus avec un cylindre solide. Dans le cas sans évaporation, un bon accord a été obtenu avec la référence, et un détachement et un élargissement de la zone de recirculation ont été observés lorsque l'évaporation est considérée. Ces résultats montrent les effets du flux de Stefan sur la dynamique du système. L'étape suivante consiste à étendre le formalisme à un formalisme faiblement compressible en utilisant la méthode basée sur la pression proposée par [START_REF] Duret | A pressure based method for vaporizing compressible two-phase flows with interface capturing approach[END_REF]. En outre, la méthode CLSVOF pour la représentation de l'interface et la méthode de projection a été modifiée pour tenir compte de l'évaporation du liquide et des effets de compressibilité, tels que la dilatation thermique. À cette fin, une équation de pression est dérivée et mise en oeuvre. Les propriétés de conservation de la masse totale du formalisme compressible sont validées en étudiant l'évolution temporelle de la masse de liquide et de gaz dans un cylindre statique 2D avec un taux d'évaporation constant. Ensuite, une configuration THI diphasique compressible 3D avec changement de phase est présentée pour montrer la capacité du formalisme à gérer les régimes d'atomisation.

Dans le chapitre suivant, notre formalisme faiblement compressible est mis en oeuvre pour analyser les processus d'évaporation et de mélange turbulent scalaire. Pour cette partie, de petits changements dans le formalisme ont été introduits pour considérer la masse molaire de la vapeur et la capacité ther-mique du mélange dans le calcul de la densité et de la température. En outre, deux simulations sont effectuées avec une configuration THI compressible à deux phases avec changement de phase, en faisant varier la température initiale du liquide pour avoir différents régimes d'évaporation. En estimant la température et la fraction massique de vapeur à l'état d'équilibre, les effets du régime d'évaporation sont observés. Dans le cas d'un taux d'évaporation important, la température d'équilibre est plus petite que la température initiale du liquide. Le comportement inverse est observé dans le cas d'un faible taux d'évaporation. Ceci est attribué à la quantité d'énergie consommée par le liquide pendant l'évaporation.

Le processus d'évaporation est étudié en comparant l'évolution temporelle de la température du liquide et du gaz, du taux d'évaporation et de la fraction massique de la vapeur. Tout d'abord, comme prévu, nous observons que la température du liquide et du gaz tend vers sa valeur d'équilibre, mais le comportement est significativement affecté par le taux d'évaporation. Ensuite, le processus de mélange turbulent scalaire est analysé via l'évolution temporelle de la variance et des formes PDF de la fraction de masse de vapeur. Ceci nous a permis d'observer les différents états du mélange turbulent. Depuis le stade initial, où le mélange est principalement contrôlé par la diffusion de la couche de vapeur près de l'interface et son évolution vers une forme similaire à une fonction de Dirac à l'équilibre. Ici, des résultats significativement différents ont été obtenus par rapport aux configurations monophasique ou diphasique dispersées. Nous concluons que l'effet de l'interface et de l'évaporation sur le mélange scalaire n'est pas négligeable. Ensuite, les densités de probabilité normalisées de la fraction de masse de vapeur sont présentées, en observant que les densités de probabilité évoluent vers une forme similaire, proche d'une distribution log-normale, dans les derniers instants de la simulation. Enfin, la relation entre le taux d'évaporation et la courbure de l'interface est analysée à l'aide de l'évolution temporelle des PDFs conjointes. Ici, un taux d'évaporation plus important est observé dans les courbures positives, ce qui correspond à la section concave de l'interface. pérature non-homogènes à l'interface, différentes structures de gaz avec une pression thermodynamique indépendante, une densité de gaz qui dépend de la température et de la pression locales, etc, le nombre d'applications possibles est vaste. Bien que, plusieurs perspectives de développement numérique peuvent être étudiées pour améliorer les formalismes incompressibles et compressibles. Cette section explorera les perspectives possibles de ce travail liées aux développements numériques et aux applications physiques.

Perspectives

Applications physiques

En commençant par le formalisme incompressible, une configuration de gouttelettes qui tombent avec changement de phase pourrait aider à accéder à la gamme d'applicabilité des modèles classiques trouvés dans la littérature. Par exemple, une continuation du travail effectué par [START_REF] Duret | Evaporating droplets: comparisons between dns and modelling[END_REF], où les effets du flux de Stefan, de la chaleur latente et de la régression de l'interface dans le transfert de chaleur et de masse d'une gouttelette sont étudiés, pourrait aboutir à des corrélations plus réalistes pour les nombres de Nusselt et de Sherwood utilisés dans la modélisation de l'évaporation. Une autre application possible de notre formalisme est l'analyse des particules capturées par les gouttelettes, généralement utilisée pour étudier la pollution atmosphérique. Des améliorations significatives pourraient être apportées dans ce domaine car les modèles ne prennent pas en compte les flux de Stefan, qui devraient influencer la capture des particules car ils repoussent les particules de l'interface lorsque l'évaporation se produit.

Pour le formalisme compressible, une étude plus détaillée menée sur la configuration THI diphasique pourrait être utilisée pour proposer de nouveaux modèles à l'échelle sous-maille pour les simulations LES et RANS. De plus, l'étude de l'évolution temporelle de la PDF de la fraction masique de vapeur a montré que la forme de la PDF dans les derniers moments de la simulation n'est pas bien représentée par les distributions gaussiennes ou β, contrairement aux résultats obtenus dans les simulations monophasiques ou diphasiques dispersées. Ceci est particulièrement intéressant pour les simulations de combustion turbulente, où l'approche de la PDF supposée est couramment utilisée pour caractériser le mélange et estimer le taux de réaction chimique.

Une autre application possible de notre formalisme compressible est l'étude du "boil off" dans les conteneurs de fluides cryogéniques. Ce phénomène est important pour le stockage et le transport de combustibles liquéfiés, tels que le gaz naturel et l'hydrogène. Le "boil off" consiste en une vaporisation spontanée lorsque le combustible liquide est transporté dans un réservoir. En con-séquence, la pression augmente et une soupape de sécurité libère une partie du combustible. Dans le cas de l'hydrogène, cette libération constitue un pourcentage important de la perte de carburant.

Perspectives de développement numérique

Il existe plusieurs améliorations possibles de notre formalisme, notamment la mise en oeuvre d'un calcul cohérent des flux de masse et de quantité de mouvement dans la simulation des écoulements diphasiques. Un candidat est le schéma numérique implémenté dans ARCHER par [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF]. Mais ce schéma doit être modifié pour tenir compte des discontinuités de vitesse générées par le changement de phase. En outre, un problème essentiel de la simulation présentée dans le chapitre 4 est le coût de calcul associé. Par conséquent, de nouvelles stratégies sont nécessaires pour réduire le nombre d'extrapolations des variables et de solveurs de Poisson requis pour les simulations. En outre, pour certaines configurations, une méthode de raffinement de maillage adaptatif (AMR) pourrait aider à capturer les couches limites tout en réduisant le nombre de points de maillage dans l'ensemble du domaine.

Actuellement, les simulations avec le formalisme compressible sont limitées à un environnement fermé car une condition limite d'écoulement appropriée n'est pas encore implémentée. Ceci est un défi en raison du traitement des ondes acoustiques. L'avantage de conditions limites appropriées est que la gamme d'applications de notre formalisme peut être étendue aux jets liquides ou aux films liquides.

De plus, de nouvelles stratégies numériques permettant de considérer simultanément l'évaporation et l'ébullition peuvent améliorer considérablement notre formalisme. Ceci a déjà été fait dans la littérature avec une méthode levelset dans [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the leidenfrost effect[END_REF].

1 1

 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 State-of-the-art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.2 Level-set method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.3 Front Tracking methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.4 Hybrid methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Experimentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.2 Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.3.3 Direct numerical simulation of phase change . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 Turbulent mixing in single-phase flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.4.2 Turbulent mixing in two-phase flows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Page 1 2 1. 1

 121 Figure Page 1 Examples of industrial processes involving evaporation and turbulent mixing: a) spray combustion in a compression ignition engine [1], b) diagram of a spray drying operation [2]. . . . . . . . . . . . . . . 2 1.1 Classification of numerical methods for for the interface representation in two-phase flows [3]. . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2 (a) Interface between two materials. (b) PLIC approximation segments of the interface shown in (a) [4]. . . . . . . . . . . . . . . . . . . 7 1.3 VOF reconstruction of the solution for advection in two dimensions. (a) The original interface. (b) The original SLIC reconstruction. (c) The Hirt and Nichols reconstruction. (d) PLIC reconstruction [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Level set method; (a) interface representation; (b) level set contours for a circular interface of radius 5 units [6]. . . . . . . . . . . . . 9 1.5 Interface representation with the Front tracking method [5]. . . . . . . 1.6 Direct numerical simulations of a vaporizing isolated droplet with different methods for the interface capturing and computation of the liquid velocity field (a) [7] (b) [8] and (c) [9]. . . . . . . . . . . . . . 1.7 Direct numerical simulations of a vaporizing droplets with a lowmach formalism of the Navier-Stokes equations in the gas phase (a) [10] (b) [11]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 . 6

 26 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Control volume dV with boundary dΓ containing a portion of the interface Γ [5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 One-dimensional case studied for the spatial discretization of the diffusion term of the energy conservation equation. . . . . . . . . . . 2.3 Normalized liquid mass time evolution for a constant evaporation rate for the small Stefan flow configuration. Solid line: Theory, • : Method A, : Method B, *: Method C. . . . . . . . . . . . . . . . . . . . 2.4 Normalized liquid mass time evolution for a constant evaporation rate for the large Stefan flow configuration for the mesh of 128 2 . Solid line: Theory, • : Method A, : Method B, *: Method C. . . . . . . 2.5 a) L ∞ , b) L 1 and c) L 2 norms of the liquid mass temporal evolution for a constant evaporation rate for the strong Stefan flow configuration. Solid line: 2nd order, Dashdotted: 1st order, • : Method A, : Method B, *: Method C. . . . . . . . . . . . . . . . . . . . . . . . . . One-fluid velocity fields for our 3 methods (t = 0.001 s). . . . . . . . . 2.7 Liquid velocity fields for our 3 methods (t = 0.001 s). . . . . . . . . . . 2.8 Gas velocity field and its extension in the liquid phase for method C (t = 0.001 s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 (a) Temperature field and (b) vapor mass fraction field (t = 0.08 s). . . 2.10 Normalized square diameter versus time (s) for a T g = 573 K. Solid: 64 3 , Dashed: 48 3 , Dotted: 32 3 . . . . . . . . . . . . . . . . . . . . . 2.11 Radial temperature profile in the gaseous phase (t = 0.08 s). Solid line: theoretical solution of D 2 law, • : 32 3 , : 48 3 : 64 3 . . . . . . . . 2.12 Radial vapor mass fraction profile in the gaseous phase (t = 0.08 s). Solid line: theoretical solution of D 2 law, • : 32 3 , : 48 3 : 64 3 . . . . . 2.13 (a) Temperature field and (b) vapor mass fraction field (t = 6 ms) for a initial velocity of u l = 1 m • s -1 . . . . . . . . . . . . . . . . . . . . . 2.14 Normal evaporation rate ( ωn) and its extension around the drop (t = 6 ms). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Velocity streamlines behind the water cylinder and visualisation of vu l = 0 iso-contour in the gas phase. Right: evaporating drop, left: non-evaporating drop. . . . . . . . . . . . . . . . . . . . . . . . . . 2.16 v -u l profile over a vertical line in the center of the domain for the water convected cylinder without evaporation (a) and with evaporation (b). The interface locations are represented by the vertical dash-dotted lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Velocity streamlines behind the n-decane cylinder and visualisation of vu l = 0 iso-contour. Right: evaporating drop, left: nonevaporating drop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.18 vu l profile over a vertical line in the center of the domain for the n-decane convected cylinder without evaporation (a) and with evaporation (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Summary of a euler time step for the compressible flow solver. . . . . 3.2 Gas mass temporal evolution. Solid line: Reference solution, •: simulation results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Liquid mass temporal evolution. Solid line: Reference solution, •: simulation results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Instantaneous temperature field with interface location (Level Set 0 isocontour). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Instantaneous vapor mass fraction field with interface location (Level Set 0 isocontour). . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Top: instantaneous divergence of the real velocity field with interface location (Level Set 0 isocontour). Bottom: Visualization of the velocity divergence inside the encapsulated gas structure highlighted in red in the top figure. . . . . . . . . . . . . . . . . . . . . . . . 4.1 Iterative algorithm to determine the equilibrium conditions of the simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Examples of industrial processes involving evaporation and turbulent mixing: a) spray combustion in a compression ignition engine [1], b) diagram of a spray drying operation [2].

Figure 1 . 2 :

 12 Figure 1.2: (a) Interface between two materials. (b) PLIC approximation segments of the interface shown in (a) [4].

  .3 three different methods for the reconstruction of the interface are shown.

Figure 1 . 3 :

 13 Figure 1.3: VOF reconstruction of the solution for advection in two dimensions. (a) The original interface. (b) The original SLIC reconstruction. (c) The Hirt and Nichols reconstruction. (d) PLIC reconstruction [5].

Figure 1 . 4 :

 14 Figure 1.4: Level set method; (a) interface representation; (b) level set contours for a circular interface of radius 5 units [6].

Figure 1 . 5 :

 15 Figure 1.5: Interface representation with the Front tracking method [5].

Figure 1 . 6 :

 16 Figure 1.6: Direct numerical simulations of a vaporizing isolated droplet with different methods for the interface capturing and computation of the liquid velocity field (a) [7] (b) [8] and (c) [9].

Figure 1 . 7 :

 17 Figure 1.7: Direct numerical simulations of a vaporizing droplets with a low-mach formalism of the Navier-Stokes equations in the gas phase (a) [10] (b) [11].

Fig. 1 .

 1 Fig. 1.7 shows two examples of simulations of an evaporating liquid in a closed and weakly compressible environment. Fig. 1.7a is a single droplet in a Homogeneous isotropic turbulence environment ([11]) and fig. 1.7b is a simulation of dispersed droplets evaporating in a weakly-compressible homogeneous shear turbulence ([10]).

Figure 2 . 1 :

 21 Figure 2.1: Control volume dV with boundary dΓ containing a portion of the interface Γ [5].

Figure 2 . 2 :

 22 Figure 2.2: One-dimensional case studied for the spatial discretization of the diffusion term of the energy conservation equation.

Figure 2 .

 2 Figure 2.3 shows the time evolution of the normalized liquid mass m * = m l m 0

Figure 2 . 3 :Figure 2 . 4 :

 2324 Figure 2.3: Normalized liquid mass time evolution for a constant evaporation rate for the small Stefan flow configuration. Solid line: Theory, • : Method A, : Method B, *: Method C.

Figure 2

 2 Figure 2.5: a) L ∞ , b) L 1 and c) L 2 norms of the liquid mass temporal evolution for a constant evaporation rate for the strong Stefan flow configuration. Solid line: 2nd order, Dashdotted: 1st order, • : Method A, : Method B, *: Method C.

  Figure 2.6: One-fluid velocity fields for our 3 methods (t = 0.001 s). L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022

Figure 2 . 8 :

 28 Figure 2.8: Gas velocity field and its extension in the liquid phase for method C (t = 0.001 s).

Figure 2

 2 Figure 2.9: (a) Temperature field and (b) vapor mass fraction field (t = 0.08 s).

2 Figure 2 . 10 :Figure 2 . 11 :

 2210211 Figure 2.10: Normalized square diameter versus time (s) for a T g = 573 K. Solid: 64 3 , Dashed: 48 3 , Dotted: 32 3 .

Figure 2 . 12 :

 212 Figure 2.12: Radial vapor mass fraction profile in the gaseous phase (t = 0.08 s). Solid line: theoretical solution of D 2 law, • : 32 3 , : 48 3 : 64 3 .

Figure 2 .

 2 Figure 2.13: (a) Temperature field and (b) vapor mass fraction field (t = 6 ms) for a initial velocity of u l = 1 m • s -1 .

Figure 2 . 14 :

 214 Figure 2.14: Normal evaporation rate ( ωn) and its extension around the drop (t = 6 ms).

Figure 2 . 15 :

 215 Figure 2.15: Velocity streamlines behind the water cylinder and visualisation of vu l = 0 iso-contour in the gas phase. Right: evaporating drop, left: nonevaporating drop.

Figure 2 .

 2 Figure 2.16: vu l profile over a vertical line in the center of the domain for the water convected cylinder without evaporation (a) and with evaporation (b). The interface locations are represented by the vertical dash-dotted lines.

Figure 2 .

 2 Figure 2.18: vu l profile over a vertical line in the center of the domain for the n-decane convected cylinder without evaporation (a) and with evaporation (b).

( 3 . 51 )

 351 Looking now at the gas velocity extension, as explained in the last chapter,L. Germes MartínezDNS of phase change June 28, 2022

Figure 3 . 1 :

 31 Figure 3.1: Summary of a euler time step for the compressible flow solver.

Figure 3 . 2 :

 32 Figure 3.2: Gas mass temporal evolution. Solid line: Reference solution, •: simulation results.

Figure 3 . 3 :

 33 Figure 3.3: Liquid mass temporal evolution. Solid line: Reference solution, •: simulation results.

Figure 3 . 4 :

 34 Figure 3.4: Instantaneous temperature field with interface location (Level Set 0 isocontour).

Figure 3 . 5 :

 35 Figure 3.5: Instantaneous vapor mass fraction field with interface location (Level Set 0 isocontour).

Figure 3 . 6 :

 36 Figure 3.6: Top: instantaneous divergence of the real velocity field with interface location (Level Set 0 isocontour). Bottom: Visualization of the velocity divergence inside the encapsulated gas structure highlighted in red in the top figure.

Figure 4 . 1 :

 41 Figure 4.1: Iterative algorithm to determine the equilibrium conditions of the simulations.

Figure 4 . 2 :

 42 Figure 4.2: Instantaneous temperature fields with interface visualization (Level Set 0 isocontour) at t 2 for the two cases studied.

Figure 4 .Figure 4 . 3 :Figure 4 . 4 :Figure 4 . 5 :

 4434445 Figure 4.4 shows the volume rendering of the vapor mass fraction at the same time as figures 4.2 and 4.3. Here, we can observe the interactions of the va-

Figure 4 . 6 :

 46 Figure 4.6: Evaporation rate temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 4 . 7 :

 47 Figure 4.7: Normalized gas temperature temporal evolution for case 1 (dashed line) and case 2 (solid line).

Figure 4 . 8 :

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

Figure 4 . 9 :

 49 Figure 4.9: Gas Lewis number temporal evolution for case 1 (dashed line) and case 2 (solid line).

Fig. 4 .

 4 Fig. 4.10 shows the temporal evolution of the RMS Y* v = (Y * v -Y * v ) 2

Figure 4 . 10 :

 410 Figure 4.10: Normalized vapor mass fraction RMS (a) and max RMS (b) for case 1 (dash line) and case 2 (solid line).

Figure 4 . 3 (

 43 Figure 4.11: (a) Instantaneous PDFs of Y * 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

  4.10a). The PDF has a pronounced peak at Y * v = 0 corresponding to dry air, and a long tail containing all the possible values of the normalized vapor mass fraction. This correspond to the first instants after the development of the mass boundary layer, where the mixing is mostly governed by diffusion of the vapor layers around the interface. This stage of the scalar mixing is characterized by large temperature and vapor mass fraction gradients at the interface, explaining the rate of change observed in the temporal evolution of the T * l , T * g , Y * v and ωΓ (see figures 4.5, 4.7, 4.8, and 4.6).

Figure 4 . 13 :

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

Fig. 4 .

 4 Fig. 4.13 displays the normalized PDF of Y * v (P (Y * v )Y * v with respect to Y *

Figure 4 . 14 : 3 (

 4143 Figure 4.14: Instantaneous PDFs of Y * v s at t * 1 (solid line), t * 2 (dashed line), t * 3

Figure 4 . 16 :

 416 Figure 4.16: Instantaneous interface visualization (zero level set iso-contour) colored by the magnitude of the evaporation rate for the same time instant than fig. 4.17.

Figure 4 . 17 :

 417 Figure 4.17: Instantaneous joint PDFs for ω in function of the normalized curvature for case 1 at t * 1 , t * 2 , t * 3 and t * 4 .

Figure 4 . 18 :

 418 Figure 4.18: Graphical representation of the curvature sign in a deformed liquid structure.

  For example,

Figure 4 . 19 :

 419 Figure 4.19: Instantaneous joint PDFs for ω in function of the curvature for case 2 at t * 1 , t * 2 , t * 3 and t * 4 .

Figure 4 .

 4 Figure 4.20: (a) Champ de température et (b) champ de fraction de masse de vapeur (t = 0.08 s).

Figure 4 .

 4 Figure 4.21: (a) Champ de température et (b) champ de fraction massique de vapeur (t = 6 ms) pour une vitesse initiale de u l = 1 m • s -1 .

  4.23).

Figure 4 . 22 :

 422 Figure 4.22: Champ de température instantanée avec localisation de l'interface (isocontour du niveau 0).

Figure 4 . 23 :

 423 Figure 4.23: Champ de fraction massique de vapeur instantanée avec localisation de l'interface (isocontour du niveau 0).

Figure 4 . 24 :Figure 4 . 25 :

 424425 Figure 4.24: Champs de température instantanés avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 pour les deux cas étudiés.

Figure 4 . 26 :Figure 4 . 27 :

 426427 Figure 4.26: Champ de fraction massique de vapeur instantanée avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 pour les deux cas étudiés.

Figure 4 . 28 :

 428 Figure 4.28: Évolution temporelle du taux d'évaporation pour le cas 1 (ligne pointillée) et le cas 2 (ligne continue).

Figure 4 . 29 :

 429 Figure 4.29: Évolution temporelle de la température normalisée du gaz pour le cas 1 (ligne pointillée) et le cas 2 (ligne solide).

Figure 4 . 30 :

 430 Figure 4.30: Évolution temporelle de la fraction de masse de vapeur normalisée pour le cas 1 (ligne pointillée) et le cas 2 (ligne continue).

La figure 4 .

 4 31 montre l'évolution temporelle du RMSY * v = (Y * v -Y * v ) 2 et sa valeur maximale théorique Y * v ,max = (Y * v (1 -Y * v

Figure 4 . 31 :

 431 Figure 4.31: Fraction de masse de vapeur normalisée RMS (a) et RMS max (b) pour le cas 1 (ligne pointillée) et le cas 2 (ligne solide).

Figure 4 . 3 (

 43 Figure 4.32: (a) PDFs instantanées de Y * v à t * 1 (ligne pleine), t * 2 (ligne pointillée), t * 3 (ligne pointillée) et t * 4 (ligne pointillée), et (b) Zoom des valeurs supérieures de Y * v aux mêmes moments pour le cas 1.

1 -→ t * 4 Figure 4 . 33 : 3 (

 144333 Figure 4.33: PDFs instantanées de Y * v à t * 1 (ligne solide), t * 2 (ligne pointillée), t * 3

  Dans ce travail, des résultats prometteurs ont été obtenus pour contribuer à une meilleure compréhension des processus d'évaporation et de mélange turbulent scalaire dans un régime d'atomisation dense. Considérant les phénomènes capturés par notre formalisme, par exemple, un taux d'évaporation et une tem-L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022

  Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2.2 Interface jump conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Interface capturing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

	2	Incompressible two-phase flow simulations with phase change
		27
	2.1	Introduction	27
	2.2	Governing equations	27
	2.2.1 2.3	Numerical Methods	33
	2.3.1

2.3.2 Projection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.3 Energy transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3.4 Vapor mass fraction transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.3.5 Summary of methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

  Champs de température instantanés avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 pour les deux cas étudiés. . . 4.25 Champs de fraction massique de vapeur instantanée avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 pour les deux cas étudiés. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 Champ de fraction massique de vapeur instantanée avec visualisation de l'interface (isocontour du level set égal à zéro) à t 2 pour les deux cas étudiés. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Évolution temporelle des températures normalisées du liquide (ligne pleine) et de l'interface (ligne pointillée) pour (a) le cas 1 et (b) le cas 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Évolution temporelle du taux d'évaporation pour le cas 1 (ligne pointillée) et le cas 2 (ligne continue). . . . . . . . . . . . . . . . . . . . 4.29 Évolution temporelle de la température normalisée du gaz pour le cas 1 (ligne pointillée) et le cas 2 (ligne solide). . . . . . . . . . . . . . . 4.30 Évolution temporelle de la fraction de masse de vapeur normalisée pour le cas 1 (ligne pointillée) et le cas 2 (ligne continue). . . . . . 4.31 Fraction de masse de vapeur normalisée RMS (a) et RMS max (b)

	xviii	LIST OF FIGURES
	4.24

. . . . . 4.2 Instantaneous temperature fields with interface visualization (Level Set 0 isocontour) at t 2 for the two cases studied. . . . . . . . . . . . . . pour le cas 1 (ligne pointillée) et le cas 2 (ligne solide). . . . . . . . . . 4.32 (a) PDFs instantanées de Y * v à t * 1 (ligne pleine), t * 2 (ligne pointillée), t * 3 (ligne pointillée) et t * 4 (ligne pointillée), et (b) Zoom des valeurs supérieures de Y * v aux mêmes moments pour le cas 1. . . . . . . . . . 4.33 PDFs instantanées de Y * v à t * 1 (ligne solide), t * 2 (ligne pointillée), t *

  INSA Rouen Normandie -CORIA June 28, 2022 the interface temperature; unrealistic values can be obtained. Which will results in a erroneous determination of the Dirichlet boundary condition for Y v and the evaporation rate ([START_REF] Gibou | A fourth order accurate discretization for the laplace and heat equations on arbitrary domains, with applications to the stefan problem[END_REF]).

Table 2 .

 2

		226	1.78 × 10 -5	0.046	1000	0.029	-	-
	Water	1000	1.137 × 10 -3	0.6	4180	0.018	2.3 × 10 6	373
	N-decane	750	5.65 × 10 -4	0.14	2207	0.142	3.25 × 10 5	447

3: Physical properties for the air, water and n-decane 2.4.1 Validation cases Interface regression due to a constant vaporisation rate

Table 2 .

 2 4: Comparison of the length of the recirculation zone for two different liquid.

	Liquid	Evaporation Re	L *	Ref. %Diff.
	Water	off	40 4.62 4.69	1.4
	Water	on	40 4.91 4.69	4.7
	N-decane	off	40 4.65 4.69 0.85
	N-decane	on	40 6.20 4.69	32

Table 4 .

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

	Fluid	ρ kg m 3	µ kg ms	λ W mK	C p	J kgK	M kg mol	h lv	J kg	T B (K)
	Air	25	1.88 × 10 -5	0.045	1060	0.029	-		-
	Vapor	-	-	-		2098	0.142	-		-
	N-decane	750	5.65 × 10 -4	0.14	2207	0.142	3.25 × 10 5	447
			Fluid		Sc	P r		Le
			Air	1.96 0.44 4.44
			N-decane	-	8.90	-

5

P a, ρ 0 = 750 kg • m -3 and γ l = 1.215. Additionally, the initial dimensionless numbers for the simulations are listed in table 4.2.

Table 4 .

 4 2: Dimensionless numbers at the initial state of the simulations.

Table 4 .

 4 3: Different variables at the equilibrium state for the two cases studied.

  De plus, la méthode est adaptée à une configuration THI diphasique avec L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022 changement de phase, similaire à [13, 14]. Dans ce travail, nous nous intéressons à l'inclusion des équations de conservation de l'énergie et des espèces dans le formalisme. Cela permet d'estimer simultanément la pression, la densité et le taux de vaporisation en fonction de la température et de la fraction massique de la vapeur. Par conséquent, un modèle plus réaliste des processus d'atomisation et d'évaporation est décrit.

  Développement du formalisme compressible avec transfert de chaleur et de masseComme nous l'avons mentionné dans le chapitre 2, le formalisme incompressible dédié aux simulations d'écoulements diphasiques avec changement de phase est limité aux environnements non confinés, c'est-à-dire qu'au moins une condition aux limite de type sortie libre doit être considérée. La divergence de vitesse créée à l'interface en raison du saut de vitesse doit être évacuée du domaine pour maintenir l'hypothèse d'incompressibilité. De plus, pour être conservatif au niveau de la masse dans une simulation en milieu fermé avec un liquide qui se vaporise, la densité du gaz doit varier avec le temps en fonction de la masse molaire de la vapeur. Cette dernière limitation a été traitée dans la littérature pour le formalisme incompressible en prédisant la densité du gaz au prochain pas de temps en utilisant la quantité de liquide vaporisé au pas de temps actuel ([START_REF] Calimez | Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs[END_REF]). Néanmoins, cette approche est incapable de prédire la variation spatiale de la pression ou de la densité. Par conséquent, un formalisme faiblement compressible pour les simulations d'écoulements diphasiques avec changement de phase dans un environnement fermé est considéré. De cette façon, les changements de densité et de pression induits par la vapeur générée à l'interface sont pris en compte, ce qui permet aux simulations de traiter de nombreuses structures liquides et gazeuses avec des conditions thermodynamiques indépendantes.

INSA Rouen Normandie -CORIA June 28, 2022

  Germes Martínez DNS of phase change June 28, 2022 ont montré les variations du taux d'évaporation dues aux petites déformations causées par le champ de vitesse turbulent. De plus, leurs résultats ont également montré une augmentation du taux d'évaporation due à l'interaction de la fraction de masse de la vapeur avec le champ de vitesse turbulent. Ces résultats soulignent la limitation des modèles classiques de vaporisation basés sur des gouttelettes sphériques avec un taux d'évaporation constant. La figure 4.26 montre le rendu volumique de la fraction de masse de vapeur en même temps que les figures 4.24 et 4.25. Ici, nous pouvons observer les interactions entre la fraction de masse de vapeur et le champ de vitesse turbulent dans le domaine. De même, l'accumulation de vapeur est observée avec une tendance à être complètement homogène lorsque le système atteint l'équilibre. Un autre phénomène intéressant est le sillage de vapeur derrière les plus petites structures liquides lorsqu'elles sont transportées à travers le domaine. Le même comportement est observé dans les simulations présentées dans la section 2.4.2.Il y a plusieurs résultats possibles pour la durée de vie de ces petites structures liquides, par exemple, elles peuvent rester convectées par le champ de vitesse jusqu'à ce qu'elles s'évaporent complètement ou s'effondrent avec une autre structure liquide. Ces scénarios possibles considèrent que le système n'est pas à l'équilibre. De plus, il est noté que la taille et la forme des sillages de vapeur dépendent de la taille de la structure, de sa vitesse moyenne, de la distance à une autre interface, etc. L'évolution temporelle de la température moyenne en volume du liquide (T l ) et de l'interface T Γ normalisée par la température d'équilibre pour les deux cas est présentée dans les figures 4.27. t * = t τ est le temps normalisé par le temps de retournement des tourbillons défini comme :

	τ =	k ξ
	où (•) est l'opérateur de moyenne de volume, et ξ est le taux moyen de dis-
	sipation d'énergie cinétique. Comme prévu, la température du liquide pour le
	cas 1 diminue à mesure qu'il s'approche de l'état d'équilibre (voir fig. 4.27a).
	De plus, la température du liquide suit la température à l'interface, diminuant
	malgré la température plus élevée du gaz à cause du saut thermique causé par

la chaleur latente. Ce résultat montre la précision de la GFM mise en oeuvre pour tenir compte de la condition de saut d'énergie et du couplage entre la température du liquide et celle du gaz. Au début de la simulation, les températures de l'interface et du liquide ont des amplitudes différentes en raison des sous-L. Germes Martínez INSA Rouen Normandie -CORIA June 28, 2022

  .33. Dans ce cas, on peut également trouver des valeurs de Y * v > 1 mais avec une probabilité plus faible que dans le cas 1. Ceci est dû au fait que la température d'équilibre est plus élevée que la température initiale du liquide (voir tableau 4.3). Ainsi, une raison probable pour avoir des valeurs de Y * v > 1 est que les plus petites structures liquides sont chauffées à des températures plus élevées que la température d'équilibre. Cela se produit surtout dans les premiers instants de la simulation. Les probabilités pour les valeurs de Y * v > 1 diminuent avec le temps jusqu'à devenir inexistantes dans la dernière PDF. Pour les mêmes valeurs de Y * v , on peut apprécier une influence du régime d'évaporation dans la forme des PDFs. Par exemple, au premier instant, le pic de la PDF a considérablement diminué, et le temps nécessaire pour arriver à la même Y * v (t * 1 = 7.66) est presque 13 fois plus grand que dans le cas 1. Pour le RMS (Y * v = 0.182), une augmentation similaire par rapport à la valeur initiale est observée (125%). La réduction de l'évaporation peut expliquer le temps pour arriver à la même Y * v , puisque moins de vapeur est générée à l'interface, et l'influence du flux de Stefan est plus faible. Par conséquent, le champ de vitesse de la turbulence a plus de temps pour agir sur le champ de la fraction de masse de la vapeur, ce qui explique la réduction du pic de la PDF. Le même phénomène se produit dans les autres PDFs. La deuxième PDF correspond à un temps de t * 2 = 10.7 et un Y * v = 0.137 représentant 94% de sa valeur initiale. À ce moment, le pic de la PDF continue de diminuer. Dans les troisième et quatrième PDFs (t * 3 = 17.51 et t * 4 = 41.65), la largeur a été considérablement réduite. Cela indique qu'un meilleur mélange turbulent est atteint par rapport au cas 1.

June 28, 2022

1.3 Phase change

L. Germes MartínezINSA Rouen Normandie -CORIA June 28, 2022

http://docs.mercurve.rdb.is/ L. Germes Martínez INSA Rouen Normandie -CORIA June

[START_REF] López | A volume of fluid method based on multidimensional advection and spline interface reconstruction[END_REF] 2022

Acknowledgements

moves toward its value at equilibrium 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.

In fig. 4.14b, it appears that the mean normalized vapor mass fraction in case 2 will not go beyond the value of 0.8. These results reveal the challenges of numerically reaching an equilibrium state. In this case, the low evaporation rate cause a relative long characteristic time needed to reach equilibrium. During this time, small numerical errors or dissipation inherent to the numerical scheme, which accumulate as the simulation progresses, could change the equilibrium temperature and vapor mass fraction of the system. Moreover, these changes are easier to observe in the results related to Y * v s due to its sensitivity to small changes in temperature (see eq. 2.22). Even though the equilibrium vapor mass fraction seems to not tend exactly to the one predicted with the iterative algorithm presented in section 4.4.1 for this specific configuration, we have a lot of information to describe the heat and mass transfer overall behavior in this case. Note that this work is one the first DNS to have access to such detailed information about simultaneous heating and vaporization at the interface in a turbulent configuration. One of the prospect of this work is to investigate or develop numerical methods to improve the energy and mass conservation.

The PDFs of the evaporation rate for the same times are presented in fig. 4.15. Here, the PDF evolves from a Dirac function shape at ω = 0, corresponding to the beginning of the simulation before the evaporation is activated, to another Dirac function shape when the system reaches the equilibrium. Also, an important dispersion of the evaporation rate is observed at the earlier moments. Additionally, a long tail to the right corresponding to the larger evaporation rate