Guillaume Sahut
email: guillaume.sahut@univ-grenoble-alpes.fr

Giovanni Ghigliotti

Guillaume Balarac

Manuel Bernard

Vincent Moureau

Philippe Marty

Numerical simulation of boiling on unstructured grids

Keywords: Boiling, Two-phase flows, Multiphase flows, Phase change, Bubble growth, Unstructured grids, Simplicial grids, Heat and mass transfer, Mass transfer rate, Level Set method, Ghost Fluid Method, Node-centered Finite Volume

come   L'archive ouverte pluridisciplinaire

Introduction

Boiling is a physical phenomenon which plays an important role in various industrial processes. It occurs especially in steam cycles of thermal power stations, nuclear power plants, two-phase loops used in cooling of electronic components and spray cooling processes. It is also used in the condenser heat exchangers of every refrigerating machine or heat pump. The study of boiling proves to be a major stake not only from a theoretical and academic point of view but also from a technical and industrial point of view, where any improvement of boiling heat transfer would have a major impact given the number of applications.

Boiling, and in general phase-changing flows, are a specific case of two-phase flows. From a computational point of view, two families of approaches exist to simulate twophase flows. The first approach consists in using Lagrangian markers to explicitly identify the position of the interface. Methods belonging to this family are referred to as fronttracking methods. The second consists in identifying the interface position through the use of an auxiliary field, that is advected within an Eulerian framework. Methods of this kind are referred to as front-capturing methods. The two most prominent examples of front-capturing methods are the Volume Of Fluid method (VOF) [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF], based on the advection of a volume fraction function, and the Level Set method (LS) [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], based on the advection of a function representing the signed distance function to the interface.

In addition to the method used to locate the interface, numerical methods to compute two-phase flows can be differentiated on the strategy deployed to handle discontinuities (both in physical properties of the fluids and boundary conditions) at the interface. A first approach is to alleviate the discontinuity by spreading the variation through a few cells across the interface. This eases the numerical solution of the system of PDEs, but may introduce unphysical spreading of the discontinuity. Methods of this kind are denoted smooth interface methods. Another approach consists in retaining the discontinuity at the discretization level. While this seems satisfactory from a modelling point of view, a certain number of numerical techniques need to be deployed in order to retain the discontinuity despite solving PDEs that naturally tend to smoothen the solution. Such methods are denoted sharp interface methods. One of the most prominent realizations of sharp interface methods is the Ghost Fluid Method (GFM) proposed by Fedkiw et al. [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF]. In the GFM, each physical field is divided in one field per phase, and each field of each phase is artificially extended beyond the interface in a continuous fashion. The extended values are called ghost values. They allow the computation of derivatives across the interface, since the discontinuity is artificially removed for the field of a given phase.

Direct Numerical Simulations (DNS) of boiling flows were pioneered by Welch [START_REF] Welch | Local simulation of two-phase flows including interface tracking with mass transfer[END_REF] in 1995. The author developed a two-dimensional, moving-mesh finite-volume method for a single, weakly deformable bubble. A simple interface model based on surface tension and surface energy is used. The mass transfer rate at the interface is computed using moving grid triangles. The computations are performed on interface-dependent triangular grids. The results showed basic capabilities to track interfaces with phase change. Son and Dhir [START_REF] Son | Numerical Simulation of Film Boiling Near Critical Pressures With a Level Set Method[END_REF] used the Continuum Surface Force model [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] to compute the mass transfer rate and the velocity discontinuity at the interface in the context of nucleate boiling to investigate bubble release pattern depending on the heat flux. The interface is captured by a Level Set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. Juric and Tryggvason [16] presented a numerical method to simulate liquid-vapor phase change in which the interface is followed by a Front-Tracking method [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF][START_REF] Juric | A front-tracking method for dendritic solidification[END_REF][START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]. Interfacial source terms for surface tension, mass transfer and latent heat are added as indicator functions of the interface. The authors used finite differences on cartesian grids to simulate film boiling and bubble departure. Based on the work of Welch [START_REF] Welch | Local simulation of two-phase flows including interface tracking with mass transfer[END_REF], Welch and Wilson [START_REF] Welch | A volume of fluid based method for fluid flows with phase change[END_REF] presented in 2000 two-dimensional numerical simulations of incompressible two-phase flows with phase change, using a VOF-based interface tracking method in conjunction with a mass transfer model and a model for surface tension. The authors derived a one-dimensional analytical test problem featuring a thin thermal layer propagating with the moving phase interface. The authors emphasized that this test problem isolates the ability of a method to accurately calculate the thermal layers responsible for driving the mass transfer in boiling flows. In 2007, Tanguy et al. [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] used a Level Set function to capture the interface coupled to the Ghost Fluid Method [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] to compute discontinuities at the interface, in order to simulate vaporizing flows (liquid-to-vapor phase change with multi-component species) in two-dimensional cartesian grids. The method is validated against various test-cases of vaporizing droplets yielding good agreement with analytical results. Also in 2007, Gibou et al. [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF] used the same coupling for two-dimensional film boiling simulation, again with good accuracy. These two references demonstrate the efficiency of the Level Set/Ghost Fluid Method coupling, and thus mark a milestone in two-phase flow simulations with phase change. In 2014, Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] proposed a comparison of different numerical methods suited to simulations of two-phase flows with phase change. The liquid-vapor interface is captured using the Level Set method. The accuracies of the Ghost Fluid Method (sharp interface) and the delta function method (smooth interface) are compared to compute the normal velocity jump condition. The authors showed that smoothing the velocity jump condition at the interface could lead to a misleading mass prediction, whereas the Ghost Fluid Method performed well. Moreover, the authors demonstrated that high-order extrapolation methods on the thermal field allowed performing accurate and robust simulations for a thermally controlled bubble growth. Simulations of the growth of static and rising bubbles are presented. The computations are performed on two-dimensional axisymmetric cartesian grids.

Based on the work of Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], our goal is to present and assess a numerical method for the simulation of multiphase flows with phase change on unstructured grids. The numerical developments have been implemented in the YALES2 code, a communitybased code developed by a network of several French labs within the Scientific Interest Group (GIS) "SUCCESS". YALES2 is a numerical code whose general purpose is the solving of fluid mechanics problems. It is designed for three-dimensional unstructured grids in order to simulate single-phase and two-phase flows in complex geometries. Our methodology has been designed for unstructured grids, with distributed memory parallelism in order to take advantage of distributed computational power. Our configuration of interest is a vapor bubble surrounded by a liquid phase. Both the vapor and liquid phases are of the same chemical composition. We focus on the growth of the bubble due to the liquid phase heated above the saturation temperature. In addition to the complexity of the physical phenomena encountered at the interface when phase change is not considered (action of surface tension and parameter discontinuity, density and viscosity, as well as pressure discontinuity between the two phases), boiling contributes to the pressure discontinuity and also implies a velocity discontinuity at the interface between the liquid and vapor phases. Classical numerical methods for the simulation of two-phase flows are not designed to take into account a discontinuous velocity field whose discontinuity is moreover coupled to the temperature field. Substantial efforts are then required to extend existing algorithms in order to simulate boiling flows, especially on unstructured grids. We present these developments along with some promising results concerning numerical simulations of boiling flows on unstructured grids. Section 2 introduces the physics of boiling, emphasizing the coupling between the Navier-Stokes equations and the heat equation by means of the velocity and pressure discontinuities at the interface computed as functions of the mass transfer rate. Section 3 presents our numerical method for the simulation of multiphase flows without phase change on unstructured grids in two and three dimensions. Section 4 extends the method to take into account phase change. Section 5 presents numerical results for the validation of the method in two and three dimensions, first with a fixed uniform and constant mass transfer rate, and second with a mass transfer rate computed as a function of the heat fluxes at the interface. Finally, Section 6 concludes this work and highlights some perspectives.

Governing equations and interface conditions

The liquid and vapor phases are assumed non-miscible and both phases are assumed incompressible. The subscripts 'liq' and 'vap' are used to name physical properties in the liquid and vapor phases, respectively. These properties are uniform in each phase but are usually different from one phase to the other (e.g. density, viscosity...). The Sharp Interface model is used in which the interface between the two phases is infinitely thin [START_REF] Montazeri | Sharp interface models for two-phase flows: Insights towards new approaches[END_REF]. Consequently, physical properties admit a discontinuity at the interface denoted by

[A] Γ = A liq -A vap , (1)
where A is a per-phase uniform physical field and discontinuous at the interface Γ, and A liq and A vap are taken in a given infinitesimal neighborhood of the interface.

Physical model

In each phase, mass conservation, or the incompressibility hypothesis, reads

∇ • u = 0, ( 2
)
where u is the fluid velocity in the given phase. While this hypothesis is very closely satisfied for the liquid phase, care must be taken when applying it to the vapor phase. Within the boiling context, we assume that the vapor temperature does not depart significantly from the saturation temperature, so that thermal dilatation can be neglected. Indeed, in the test cases presented in this paper, this is the case. Following the notations of Nguyen et al. [START_REF] Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF] and Gibou et al. [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF], we denote the interface normal velocity by W = Dn, where n is the interface normal vector. Mass conservation across the interface reads (see also [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF])

ṁ = ρ liq (u liq • n -D) = ρ vap (u vap • n -D) , (3)
where ṁ is the mass transfer rate (mass per unit of time and surface) across the interface accounting for phase change, ρ is the density and n is the interface normal vector pointing to the liquid phase. As a result, the magnitude of the interface normal velocity D is given by

D = u liq • n - ṁ ρ liq = u vap • n - ṁ ρ vap , (4)
leading to the discontinuity of the normal velocity at the interface, given by

[u • n] Γ = ṁ 1 ρ Γ . (5)
Moreover, in order to specify that the tangential components of the velocity are continuous across the interface, it is more convenient to rewrite Eq. ( 5) as [START_REF] Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF][START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]]

[u] Γ = ṁ 1 ρ Γ n. (6)
In each phase, momentum conservation leads to the Navier-Stokes equations

∂ρu ∂t = ∇ • (T -ρuu) + f , ( 7
)
where T is the symmetric stress tensor given by

T = -P I + 2µS, (8)
P is the pressure and µ is the dynamic viscosity of the considered phase, I is the identity tensor, S is the symmetric rate-of-strain tensor given by

S = 1 2 ∇u + (∇u) T , (9)
and f is the total body force applied to the fluid particle. Neglecting the gravitational forces and assuming a uniform surface tension coefficient [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF], the total body force reduces to

f = σκnδ Γ , ( 10
)
where σ is the surface tension coefficient, κ is the interface curvature and δ Γ is a multidimensional Dirac distribution localized at the interface [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]. Equations ( 7) and [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF] constitute the Whole-Domain Formulation of Navier-Stokes equations [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]. Refering to Fig. 1, we define point P moving at the interface velocity W = Dn, i.e. the relative velocity of the fluid with respect to P is u -Dn. Rewriting Eq. ( 7) in the reference frame moving with point P leads to

∂ ∂t {ρ (u -Dn)} = ∇ • {T -ρ (u -Dn) (u -Dn)} + σκnδ Γ , (11)
where

S = 1 2 ∇ (u -Dn) + (∇ (u -Dn)) T . ( 12
)
The pressure discontinuity at the interface is obtained by integrating Eq. ( 11) over the volume V of thickness δ surrounding point P and taking the limit of the integral when δ tends to zero. Indeed, since V is a thin band centered on the interface, this limit gives 5 the expression of the discontinuity of each volumetric term appearing in Eq. [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF]. This demonstration has been given in [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF] without phase change. We now extend it to take phase change into account. In order for the volume V to follow the motion of point P , we use the Reynolds transport theorem to integrate the left-hand side of Eq. ( 11) over

V , i.e. V ∂ ∂t {ρ (u -Dn)} dV = d dt V ρ (u -Dn) dV - S (Dn • n S ) ρ (u -Dn) dS. ( 13
)
In the limit δ → 0 + , the first term of the right-hand side of Eq. ( 13) vanishes. The complete integral of Eq. ( 11) over V then reads (14) Moreover, the integrals over S in Eq. ( 14) reduce to the sum of the integrals over S 1 and S 3 . Since S 1 is located in the liquid phase and S 3 , in the vapor phase, and since n S1 = n and n S3 = -n, Eq. ( 14) is rewritten

0 = -[ρ (u -Dn) (u • n -2D)] Γ -[P ] Γ n + 2 [µS] Γ n + σκn, ( 15
)
where δ is taken small enough for κ to be assumed uniform on the portion of the interface included in V . Taking the normal component of Eq. ( 15) and simplifying the first term of the right-hand side by means of Eq. (3) leads to

[P ] Γ = σκ + 2 [µn • Sn] Γ -ρ (u • n -D) 2 Γ , ( 16
)
where S is initially given by Eq. ( 12) and can be reduced here to Eq. ( 9) since n is uniform on S 1 and S 3 . Finally, using Eqs (3) and ( 5), the last term of the right-hand side of Eq. ( 16) is rewritten

ρ (u • n -D) 2 Γ = ṁ2 1 ρ Γ , (17)
leading to the pressure discontinuity at the interface,

[P ] Γ = σκ + µn • ∇u + (∇u) T n Γ -ṁ2 1 ρ . (18)
In the Whole-Domain Formulation, the contribution of surface tension forces to the pressure discontinuity at the interface is already taken into account by means of the total body force f in Eqs. ( 7) and [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF]. Moreover, as stated in [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF], independently of the formulation, the contribution of viscous effects is already taken into account when computing the divergence of the stress tensor T in Eq. [START_REF] Dalmon | Fluids-membrane interaction with a full eulerian approach based on the level set method[END_REF]. Only the contribution of the mass transfer rate needs to be added separately. The main difficulty of this formulation is the need of an accurate discretization of subgrid quantities, often leading to the smearing of the interface [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF].

Another formulation of the Navier-Stokes equations for two-phase flows is the Jump-Condition Formulation [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF] in which surface tension forces are not directly included in the equations. In this formulation, Eq. ( 7) is rewritten

∂u ∂t + ∇ • (uu) = - ∇P ρ + 1 ρ ∇ • µ ∇u + (∇u) T , (19)
where the density is moreover extracted from the derivatives, emphasizing that it is constant and uniform in each phase. The contributions of surface tension forces and mass transfer rate to the pressure discontinuity at the interface given by Eq. ( 18) are added as external ingredients. In this work, we use the Jump-Condition Formulation for compatibility with the Sharp Interface model [START_REF] Montazeri | Sharp interface models for two-phase flows: Insights towards new approaches[END_REF]. Energy conservation can be formulated by means of different primitive variables as the internal energy or the enthalpy. The internal energy is favored to describe thermodynamics phenomena in closed systems, and temperature (or enthaly) is more commonly used for description of quasi-isobaric phenomena (usually open systems) [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. In this work, we focus on simulations of boiling in open systems, we then express energy conservation with the temperature variable. In each phase, energy conservation reads

∂T ∂t + u • ∇T = 1 ρc p ∇ • (λ∇T ) , ( 20
)
where T is the temperature, c p is the heat capacity at constant pressure and λ is the thermal conductivity of the considered phase. Without phase change, energy conservation implies that the heat flux at the interface -λ∇T | Γ is entirely transmitted from one phase to the other, i.e. one has [-λ∇T • n] Γ = 0. When phase change is considered, a portion of the heat flux, equal to [-λ∇T • n] Γ , is absorbed by the interface neighborhood in order to provide the latent heat L v required to break liquid molecular bonds and transform some quantity ṁ of liquid into vapor per interface surface and time units. The remaining part of the heat flux is transmitted from one phase to the other. This phenomenon leads to the following power balance,

[-λ∇T • n] Γ = ṁL v , (21)
accounting for the latent heat of vaporization as an energy sink at the interface. Equation ( 21) acts as a definition of the mass transfer rate ṁ. Moreover, we make the assumption that the temperature is continuous at the interface and that the interface temperature T Γ is equal to the uniform and constant saturation temperature T sat of the fluid, i.e.

T Γ = T sat . (22)

Interface capturing method

The liquid-vapor interface is captured using the Level Set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. The signed distance function φ to the interface Γ is defined at time t for x belonging to the computational domain Ω as

φ (t, x) =      -min xΓ∈Γ(t) x -x Γ , if x ∈ vapor phase, 0, if x ∈ Γ(t), min xΓ∈Γ(t) x -x Γ , if x ∈ liquid phase, (23)
where the set {x ∈ Ω | φ (t, x) = 0} captures the interface at time t. In order to accurately capture the interface throughout the simulation, the signed distance function φ is advected by the interface velocity u Γ using the standard advection equation

∂φ ∂t + u Γ • ∇φ = 0, (24)
where u Γ is equivalently given by [START_REF] Nguyen | A boundary condition capturing method for incompressible flame discontinuities[END_REF][START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]]

u Γ = u liq - ṁ ρ liq n (25)
and

u Γ = u vap - ṁ ρ vap n. (26)
Equations ( 25) and [START_REF] Sahut | Simulation numérique de l'ébullition sur maillages non structurés[END_REF] denote that the transport of the interface can be formulated in both phases. The interface normal vector n is given by

n = ∇φ ∇φ (27)
and the interface curvature κ, by

κ = -∇ • n, (28)
which is rewritten in the form given by Goldman [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF],

κ = ∇φ T H (φ) ∇φ -∇φ 2 Tr (H (φ)) ∇φ 3 , (29)
where H (φ) is the hessian matrix of φ and Tr is the trace operator. If one expresses the interface velocity by means of Eq. ( 26), then Eq. ( 24) is rewritten

∂φ ∂t + u vap • ∇φ = ṁ ρ vap , (30)
assuming that ∇φ = 1. In order to solve Eq. ( 30) on grid nodes, the interface velocity u Γ , while physically defined only at the subgrid interface, is computed on grid nodes.

Numerical method to simulate multiphase flows without phase change on unstructured grids

This work is based on YALES2, a library of finite-volume-based solvers developed for realistic turbulent two-phase flow simulations with low Mach numbers [START_REF] Moureau | Design of a massively parallel cfd code for complex geometries[END_REF]. The YALES2 solver has been validated for various applications such as combustion [START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flame-pdf modeling[END_REF][START_REF] Guédot | Design of implicit high-order filters on unstructured grids for the identification of large-scale features in large-eddy simulation and application to a swirl burner[END_REF], bio-mechanics [START_REF] Chnafa | Image-based large-eddy simulation in a realistic left heart[END_REF], hydro-electricity [START_REF] Wilhelm | Analysis of head losses in a turbine draft tube by means of 3D unsteady simulations[END_REF], or wind energy [START_REF] Bénard | Large-eddy simulation of wind turbines wakes including geometrical effects[END_REF]. It is able to handle efficiently unstructured meshes with several billions of elements, thus enabling the Large Eddy Simulation (LES) of laboratory and semi-industrial configurations.

YALES2 uses a level set method to capture the interface motion [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF]. The step of time advection of the level set function does not preserve the profile of the signed distance function to the interface. After advection of the level set function, the function φ is then reshaped, or reinitialized, to the signed distance function to the advected interface using the geometrical approach proposed by Janodet et al. [START_REF] Janodet | An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows[END_REF]. The interface curvature is computed using Eq. ( 29) in which high-order differential operators provided by Bernard et al. [START_REF] Bernard | A framework to perform highorder deconvolution for finite-volume method on simplicial meshes[END_REF] are used.

Another important ingredient of the method is the pressure-velocity coupling. The incompressible Navier-Stokes equation ( 19) is solved in both liquid and vapor phases using the projection method based on fractional time steps developed by Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and improved by Kim and Moin [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF]. The methodology is first presented without phase change. The extension to phase change will be detailed in the next section. Knowing u n and P n-1 2 , the first step of the projection method used to define u n+1 , the velocity field at next iteration, is the computation of a velocity predictor u * given by

u * -u n ∆t = -∇ • (u n u n ) + 1 ρ n ∇ • µ n ∇u n + (∇u n ) T , (31)
which is a priori not divergence-free. The second step is the resolution of a Poisson equation to compute the pressure P n+ 1 2 , given by

∇ • ∇P n+ 1 2 ρ n = ∇ • u * ∆t , (32)
where P n+ 1 2 is the only unknown. The last step is the correction of the velocity predictor u * by means of the updated pressure P n+ 1 2 to obtain the updated velocity u n+1 , i.e.

u n+1 -u * ∆t = - ∇P n+ 1 2 ρ n , (33)
ensuring ∇ • u n+1 = 0 by Eq. ( 32). In Eqs. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF][START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF], the density and viscosity, despite being uniform and constant in each phase, are time-dependent on one given grid node since their values are altered by the passage of the interface. The left-hand side of Eq. ( 32) is composed of two differential operators, the gradient operator applied to P and the divergence operator applied to ∇P/ρ. The computation of a derivative in one phase should be performed using only the physical quantities defined in the said phase, but the numerical discretization of the differential operator requires using also grid nodes located in the other phase. Then the jump of the quantity at the interface has to be substracted in the discretization. Consequently, one has to compute [P ] Γ and [∇P/ρ] Γ .

Without phase change, and neglecting the viscous stress tensor jump, the pressure jump condition at the interface is given by the Laplace pressure, i.e.

[P ] Γ = σκ.

The jump of the quantity ∇P/ρ can be formally defined as

∇P ρ Γ = 1 ρ liq ∇P liq Γ - 1 ρ vap ∇P vap Γ , (35)
where P liq Γ and P vap Γ are taken infinitely close to the interface, respectively in the liquid and vapor phases. Without phase change, Eq. ( 35) can be simplified using the continuity of the velocity at the interface. Indeed, by linearity of the jump operator, Eq. ( 19) implies

∂u ∂t Γ + [∇ • (uu)] Γ = - ∇P ρ Γ + 1 ρ ∇ • µ ∇u + (∇u) T Γ , (36)
where the continuity of the velocity and its derivatives at the interface implies that the two jumps of the left-hand side are null. Moreover, as stated in [START_REF] Dalmon | Fluids-membrane interaction with a full eulerian approach based on the level set method[END_REF], the two jumps of the right-hand side are directed toward the interface normal vector and given by

n • ∇P ρ Γ = µn • ∇ 2 u ρ Γ . (37)
Dalmon et al. [START_REF] Dalmon | Fluids-membrane interaction with a full eulerian approach based on the level set method[END_REF] also emphasize that, similarly to Eq. ( 18) for the pressure jump, the contribution of viscous effects to the jump of the pressure gradient is already taken into account when computing the divergence of the stress tensor T in Eq. [START_REF] Dalmon | Fluids-membrane interaction with a full eulerian approach based on the level set method[END_REF]. Indeed, the jump of the pressure gradient can be decomposed as

n • ∇P ρ Γ = n • ∇P ρ predictor, Eq. (31) Γ + n • ∇P ρ Poisson, Eq. (32) Γ , (38)
where the exponents of the right-hand side indicate the step of the projection method in which the corresponding jump component is imposed. Consequently, provided that the velocity predictor u * is computed by Eq. ( 31), i.e. with µ inside the divergence, then the remaining part of the jump of the pressure gradient that has to be enforced in Eq. (32), i.e. the second term of the right-hand side of Eq. ( 38), simplifies to [START_REF] Dalmon | Fluids-membrane interaction with a full eulerian approach based on the level set method[END_REF] ∇P ρ

Poisson, Eq. ( 32)

Γ = 0. ( 39
)
To detail the numerical implementation, we refer to Fig. 2 representing the control volume V p associated to grid node x p crossed by the interface Γ. As a preliminary, we give the spatial discretization of the divergence of the velocity used in the finite volume framework of YALES2 to solve Eqs. ( 31) and (32) at node x p belonging to the liquid phase. One has To each pair of nodes are associated two facets of the control volume, each linking the barycenter of a neighboring element to the middle of the pair. The sum of the fluxes of a physical field at the two facets is equal to the flux of the same physical field at the "face" linking the two barycenters of the neighboring elements to the pair, represented in dense red dashes. If this face is denoted S i and if n i is its outward normal vector to the control volume of p, then the vector Ap,q i represented in red is equal to S i n i . For more clarity, S i and n i are not represented. The interface Γ is represented in blue.

∇ • u| p = 1 V p N i=1 u p + u qi 2 • A p,qi , (40)
• p • q 1 • q 2 • q 3 • q 4 • q 5 + + + + + • • • • • A p,q2 A p,q3 A p,q4 A p,q5 A p,q1 Γ V p
where N is the number of neighboring nodes of node x p , A p,qi is the oriented face of the control volume V p as illustrated in Fig. 2, and where the phase and temporal indices are omitted. We also give the spatial discretization of the laplacian of the pressure used to solve Eq. ( 32) at node x p , one has

∇ • ∇P n+ 1 2 ρ n p = 1 V p N i=1 1 ρ n p,qi ∇ p,qi P n+ 1 2 • A p,qi , (41)
where

ρ n p,qi = ρ liq if x qi ∈ liquid phase, ρ Γ := θρ liq + (1 -θ) ρ vap otherwise, (42)
θ ∈ [0; 1] is the distance between x p and the interface along the pair (x p , x qi ), relatively to the length of the pair, and ∇ p,qi P n+ 1 2 is the pressure gradient along the pair (x p , x qi ) defined below. We now describe the method used to impose the jumps at the interface given by Eqs. ( 34) and (39) in Eq. ( 41). The only node pair crossed by the interface is the pair (x p , x q1 ) where node x p is located in the liquid phase and node x q1 in the vapor phase. Following the spirit of a sharp interface method, the jumps must be imposed at the intersection point between this node pair and the interface. In order to do so, we use the Ghost Fluid Method (GFM) [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], and two pressure fields (one per phase) are considered. This is only an intermediate step necessary for the derivation. Once the jumps are applied at the interface, only one pressure field remains in the whole domain, with liquid pressure values in the liquid phase and vapor pressure values in the vapor phase (see Eq. ( 51)). A liquid pressure field, P liq , is defined in the liquid phase and in a vapor region close to the interface, in order to be able to define differential operators in the liquid neighborhood of the interface (e.g., at node x p ). Similarly, a vapor pressure field, P vap , is defined in the vapor phase and in a liquid region close to the interface. Since we discretize the pressure laplacian at node x p which is located in the liquid phase, we use the liquid pressure field in Eq. ( 41), even at node x q1 . Consequently, the pressure jump at node x q1 must be computed and applied in Eq. (41) to remove the (liquid) pressure discontinuity on the same node in the gradient operator. Equations ( 35) and (39) lead to

[∇P ] Γ = [ρ] Γ ρ liq ∇P liq Γ . (43)
As a result, the first-order Taylor expansion of the pressure jump at grid node x qi from the pressure jump at the interface reads

[P ] qi = [P ] Γ + (1 -θ) [ρ] Γ ρ liq ∆x p,qi • ∇P liq Γ . ( 44
)
Then the pressure gradient ∇ p,qi P liq used in Eq. ( 41) is given along any node pair (x p , x qi ) located in the liquid phase by

∇ p,qi P liq = P liq qi -P liq p ∆x p,qi 2 ∆x p,qi , (45)
and along any node pair crossed by the interface, as for the pair (x p , x q1 ), by

∇ p
where the 'liq' superscript is added to emphasize that Eq. ( 41) is solved using only pressure values compatible with the liquid phase (either physical values in the liquid phase, or vapor values in the vapor phase to which the pressure jump is substracted), and where the temporal dependence is removed for clarity. We make the assumption that the pressure gradient ∇P liq Γ is directed along node pair (x p , x q1 ) and is thus equal to ∇ p,q1 P liq . Using Eq. ( 46), Eq. ( 48) is then rewritten

ρ Γ ρ liq P liq q1 -P liq p ∆x p,
u * p + u * q1 ∆t • A p,q1 - 1 ρ Γ [P ] Γ ∆x p,q1 2 ∆x p,q1 • A p,q1 . (51)
The right-hand side of Eq. ( 51) is the contribution (multiplied by the control volume V p ) of node pair (x p , x q1 ) to the computation of the pressure laplacian at node x p discretized in Eq. (41) where the pressure jump is computed and applied at the interface. For node pairs not crossed by the interface, the second term of the right-hand side of Eq. (51) vanishes. Equation ( 32) is then solved using the Deflated Pre-Conjugate Gradient (DPCG) solver [START_REF] Malandain | Simulation massivement parallèle des écoulements turbulents à faible nombre de Mach[END_REF][START_REF] Malandain | Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines[END_REF]. The same methodology is used to apply the pressure jump at the interface for the computation of the pressure gradient in the correction step of the velocity predictor given by Eq. ( 33).

Numerical implementation to take into account phase change

The difference between two-phase flow simulations with and without phase change is the discontinuity of the velocity field at the interface in the case of phase change given by Eq. ( 6). Moreover, phase change requires the resolution of the heat equation in order to compute the mass transfer rate ṁ at the interface.

Velocity-pressure coupling

With phase change, the projection method given by Eqs. ( 31), ( 32) and ( 33) is still applied in each phase. Nevertheless, since the velocity is discontinuous at the interface, special care has to be taken when differentiating across the interface. In order to overcome this difficulty, the Ghost Fluid Method (GFM) [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] is again used. The velocity field is divided into one field per phase. These velocity fields are defined in domains that overlap at least a few numerical cells around the interface. Both velocity fields are artificially extended (or extrapolated) beyond the interface with the method detailed below (the method does not depend on the GFM). Consequently, we define a liquid velocity field u liq and a vapor velocity field u vap . We define both velocity fields u liq and u vap in the whole domain. Physical values of the liquid (vapor) velocity are contained in the liquid (vapor) part of u liq (u vap ), and the vapor (liquid) part of u liq (u vap ) is used to store the liquid (vapor) ghost values.

The liquid velocity field is extrapolated across the interface in the direction of the interface normal vector n Γ solving the constant extrapolation equation given by Aslam [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF],

∂u liq ∂τ + ∇u liq • n Γ = 0, ( 52
)
where τ is a pseudo-time emphasizing that this is purely a numerical step. Equation ( 52) is solved only on the vapor nodes of a narrow band around the interface. The same extrapolation is computed for u vap on the liquid nodes of the narrow band. Once Eq. ( 52) is solved for both liquid and vapor phases, the velocity divergence for both phases can be computed close to the interface in Eqs. ( 31) and [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF]. Moreover, the superposition at one node of the fields of the two phases enables the computation of the jump at this node.

In the projection method, the pressure jump at the interface is given by Eq. ( 18). Due to the discontinuity of the velocity at the interface, Eq. (39) does not hold anymore, and one has

[P ] qi = [P ] Γ + (1 -θ)∆x p,qi • ρ vap 1 ρ ∇P Γ + [ρ] Γ ρ liq ∇P liq Γ . (53)
Following the same methodology as in the previous section, Eq. (51) becomes

1 ρ Γ P vap q1 -P liq p ∆x p,q1 2 ∆x p,q1 • A p,q1 = 1 2 u * ,liq p + u * ,liq q1 ∆t • A p,q1 - 1 ρ Γ [P ] Γ ∆x p,q1 2 ∆x p,q1 • A p,q1 -(1 -θ) ρ vap ρ Γ ∆x p,q1 • [u * ] Γ -u n+1 Γ ∆t ∆x p,q1 2 ∆x p,q1 • A p,q1 , (54)
where Eq. ( 33) is used to replace

[∇P/ρ] Γ by [u * ] Γ -u n+1
Γ /∆t. The right-hand side of Eq. ( 54) is the contribution (multiplied by the control volume V p ) of node pair (x p , x q1 ) to the computation of the pressure laplacian at node x p discretized in Eq. (41) where the pressure and velocity jumps are computed and applied at the interface in the case of phase change. For node pairs not crossed by the interface, the second and third terms of the right-hand side of Eq. (54) vanish. Moreover, in Eq. ( 54), the jump of the velocity predictor, due to the jump of u n , is computed as

[u * ] Γ = [u * ] p + [u * ] q1 2 , (55)
using the physical and ghost values discussed above, whereas the jump of the future velocity u n+1 Γ is computed using Eq. ( 6) which actually accounts for the phase change phenomenon by means of the mass transfer rate. The same methodology is used to apply the pressure and velocity jumps at the interface for the computation of the pressure gradient in the correction step of the liquid velocity predictor given by Eq. [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]. Similarly, the projecton method is applied in the vapor phase to update the vapor velocity field.

Heat equation and mass transfer rate

The mass transfer rate ṁ is computed at time n + 1. Since the mass transfer rate has physical meaning only at the interface, special care is taken for its computation. The mass transfer rate defined by Eq. ( 21) can be expanded to

ṁΓ = -λ liq ∇T liq | Γ + λ vap ∇T vap | Γ L v • n Γ . (56)
In order to compute ṁΓ precisely at the interface, it is necessary to evaluate accurately the temperature gradient of each phase at the interface position. The main difficulty is then to avoid differentiating across the interface. Indeed, the computation at the interface position or at nodes of the first layer, as x p in Fig. 2, would request informations from both sides of the interface, which is not possible. In this work, gradients are extended from nodes of the second layer, x q3 and x q4 in Fig. 2, to the interface location using Taylor series expansions. These expansions are based on high-order differential operators introduced in [START_REF] Bernard | A framework to perform highorder deconvolution for finite-volume method on simplicial meshes[END_REF]. The computation of the liquid temperature gradient at the interface location is detailed hereinafter. The procedure is the same for the vapor phase. First of all, second-order accurate gradient G O2 and first-order accurate Hessian H O1 are computed at nodes of the second layer, x q3 and x q4 . These operators are then extrapolated at nodes of the first layer, x p , with second-and first-order accuracy, respectively,

   G O2 p,qj (T liq ) = G O2 qj (T liq ) + H O1 qj (T liq ) • x p -x qj , H O2 p,qj (T liq ) = H O1 qj (T liq ) . (57)
Differential operators at node x p are then obtained by averaging Taylor series expansions from nodes x q3 and x q4 . Once differential operators are known at nodes of the first layer, another expansion is performed from node x p to the interface location, to determine the liquid temperature gradient at the interface,

∇T liq | p,Γ = G O2 Γ,p (T liq ) = G O2 p (T liq ) + H O1 p (T liq ) • (x Γ -x p ) . ( 58
)
The same methodology applied in the vapor phase yields ∇T vap | q1,Γ . The interface normal vector is computed at the same interface location by

n p,q1;Γ = (1 -θ) n p + θn q1 . (59)
The mass transfer rate is then computed precisely at the interface location crossing the node pair (x p , x q1 ) by

ṁp,q1 = -λ liq ∇T liq | p,Γ + λ vap ∇T vap | q1,Γ L v • n p,q1;Γ . (60)
Finally, since values of ṁ are needed on nodes to solve the level set advection equation and Navier-Stokes equations, the value of ṁ stored at node x p is computed as the simple average of the values ṁp,qj computed on node pairs x p , x qj crossing the interface. The heat equation is solved using the following semi-implicit scheme [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]

ρ i c p,i ∆t T n+1 i -∇ • λ i ∇T n+1 i = ρ i c p,i ∆t T n i -ρ i c p,i u n i • ∇T n i , (61)
where the subscript i denotes the phase. As in [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], we make the assumption that boiling occurs at a uniform and constant saturation temperature. The immersed Dirichlet boundary condition given by Eq. ( 22) is then used at the interface when solving Eq. ( 61). In Eq. ( 61), the advection term ρ i c p,i u n i • ∇T n i is included in the right-hand side as a source term. The term semi -implicit refers to the fact that the advection term is still computed explicitly. In order to compute the diffusion term at time n + 1 and at node x p with the condition at the interface given by Eq. ( 22), one needs to compute a ghost liquid temperature value T liq,G q1 on node x q1 . Consequently, as for the velocity, we define two temperature fields T liq and T vap in order to ease the computation of ghost values by the GFM. The discretization of the laplacian of liquid temperature is given at node x p by

∇ • (λ liq ∇T liq )| p = 1 V p λ liq T liq,G q1 -T liq p ∆x p,q1 2 ∆x p,q1 •A p,q1 + 1 V p 5 j=2 λ liq T liq qj -T liq p ∆x p,qj 2 ∆x p,qj •A p,qj , (62)
where the ghost value T liq,G q1 is computed by linear extrapolation between the node x p and the interface whose temperature is set to T sat , i.e.

T liq,G q1 = T liq p + T sat -T liq p θ . (63)
The explicit advection term is computed using the standard second-order accurate finitevolume gradient operator for liquid nodes whose neighboring nodes are all located in the liquid phase. For liquid nodes which have a neighbor in the vapor phase, such as x p , the strategy previously presented, based on high-order differential operators from [START_REF] Bernard | A framework to perform highorder deconvolution for finite-volume method on simplicial meshes[END_REF], for the computation of the temperature gradients used to define the mass transfer rate in Eq. (60), is reused to compute ∇T | liq,n p

. Finally, in order to compute the temporal derivative, the high-order differential operators from [START_REF] Bernard | A framework to perform highorder deconvolution for finite-volume method on simplicial meshes[END_REF] are used to extrapolate the temperature field across the interface. This provides nodes that are swept over by the interface between times n and n + 1 with valid physical temperature values both at times n and n + 1 [START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF].

The complete linear system used to solve Eq. ( 61

θ ∆x p,q1 2 • A p,q1 T sat , (64)
which is solved using the BiCGStab(2) linear solver [START_REF] Sleijpen | BiCGStab( ) for linear equations involving unsymmetric matrices with complex spectrum[END_REF]. While the relative distance θ appears in both the left-and right-hand sides of Eq. ( 64) in denominators, our simulations did not require the use of a threshold on its value for stability. The same method is used to update the vapor temperature field.

Remarks on the interface description

Once u vap and ṁ are known, Eq. ( 30) is used to advect the signed distance function φ to the interface, in which the mass transfer rate is considered as a source term. As in [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], the convention stating that the interface is advected by the velocity of the enclosed phase is used. Since we focus on bubble growth, the liquid-vapor interface is then advected by the vapor velocity field. Figure 3 presents the actual full algorithm used in our solver.

Moreover, we emphasize here that u vap and ( ṁ/ρ vap ) n are two velocity vectors of different but dependent physical origins. The vapor velocity is due to the flow dynamics and the mass transfer rate is due to phase change. The vapor velocity field acts on the mass transfer rate by means of the vapor temperature field, and the mass transfer rate intervenes in the projection method used to update the vapor velocity field. As a result, these two velocity fields are coupled. This emphasizes the need of an overall accurate numerical method to simulate boiling flows since a small numerical error in the computation of one quantity can rapidly propagate to others, and eventually lead to an abnormal deformation of the interface.

Numerical results

Our numerical method is assessed on two benchmarks already used in Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] in the case of two-dimensional axisymmetric cartesian grids. These test-cases deal with bubble growth and give analytical expressions of the bubble radius at any time of the simulation. The first one consists in the growth of a bubble with a fixed mass transfer rate (the heat equation is not solved). The second one uses the evolution of the temperature field to define the mass transfer rate computed at the interface using Eq. ( 21), and therefore determine the time-dependent bubble radius.

Growth of a static bubble with a constant and uniform mass transfer rate

Our numerical method is evaluated on the case of a 2D static growing water bubble with a fixed mass transfer rate from [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. The computational domain is a square of side length L = 8 × 10 -3 m. The initial bubble radius is R 0 = 1 × 10 -3 m and the imposed

Compute φ n+1

Eq. [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF] Reinitialize φ n+1 → φ n+1 [START_REF] Janodet | An unstructured conservative level-set algorithm coupled with dynamic mesh adaptation for the computation of liquid-gas flows[END_REF] and compute n n+1 and κ n+1 [START_REF] Bernard | A framework to perform highorder deconvolution for finite-volume method on simplicial meshes[END_REF] Eqs. ( 27) and ( 29 Eq. ( 52)

Level Set

Temperature † Mass transfer rate

Projection method

Initialization of φ 0 , T 0 liq , T 0 vap , u 0 liq = u 0 vap = 0 and computation of ṁ0, † Γ Figure 3: Algorithm used in our solver. The computations of the narrow band around the interface and of the relative distance θ are omitted. They are actually performed at initialization, and twice in the temporal loop: once after the advection of the Level Set function and once after its reinitialization. Moreover, the computation of the time step and the exit condition of the temporal loop are also omitted. mass transfer rate is ṁ = 1 × 10 -1 kg m -2 s -1 . The theoretical bubble radius R th (t) is given at any time t ≥ 0 by

ρ liq 1 × 10 3 kg m -3 ρ vap 1 kg m -3 σ 7 × 10 -2 N m -1 µ liq 1 × 10 -3 kg m -1 s -1 µ vap 1.78 × 10 -5 kg m -1 s -1
R th (t) = R 0 + ṁ ρ vap t. ( 65
)
The simulations are performed until final time t 1 = 1 × 10 -2 s needed for the bubble radius to double. The other physical parameters of interest are listed in Table 1. The computations have been performed on the four following grid cell sizes : 4 × 10 -4 m, 2 × 10 -4 m, 1 × 10 -4 m, 5 × 10 -5 m, where the grid cell size is approximated on unstructured grids by a commercial meshing software. The coarsest grid used is shown in Fig. 4.

In order to evaluate the accuracy of the method, relative errors on the bubble radius are

ξ 2 (R) = (x1,x2)∈Λ (R (x Γ ) -R 1 ) 2 (x1,x2)∈Λ R 2 1 (66) and ξ ∞ (R) = 1 R 1 max (x1,x2)∈Λ |R (x Γ ) -R 1 | . ( 67
)
In these equations, Λ is the set of node pairs crossed by the interface, x Γ is the interface location interpolated on node pair (x 1 , x 2 ) ∈ Λ by

x Γ = (1 -θ x1,x2 ) x 1 + θ x1,x2 x 2 , (68)
with θ x1,x2 the relative distance involved in Eq. ( 42), and R (x Γ ) is given by

R (x Γ ) = x Γ , (69)
since the bubble is centered on the origin of the (e x , e y ) frame. Figure 5a shows the interface location at initial and final times as well as the liquid and vapor velocity fields at final time, on the finest unstructured grid considered. One can see the excellent agreement between the computed and theoretical interfaces at final time. As expected, the vapor is almost static while the liquid is ejected radially from the interface. Figure 5b 20

shows the relative errors, ξ 2 (R) and ξ ∞ (R), at final time. The method has a convergence rate close to one for the final bubble radius with respect to grid refinement for both relative errors, ξ 2 (R) and ξ ∞ (R). On the finest grid, the relative error ξ ∞ (R) is below 1%. This test-case has also been performed on three-dimensional unstructured grids. Figures 6a and 6b show the interface at initial and final times in three dimensions on the finest grid considered. Figure 6c shows the corresponding relative errors, ξ 2 (R) and ξ ∞ (R), with grid refinement. In three dimensions, the method also has a convergence order close to one for the final bubble radius with respect to grid refinement for both relative errors, ξ 2 (R) and ξ ∞ (R). We emphasize here the excellent agreement between the theoretical and numerical results quantified through the stringent error based on the infinity norm, ξ ∞ (R).

The proposed numerical method is then able to simulate phase change with high accuracy by means of a fixed mass transfer rate on two-and three-dimensional unstructured grids. Indeed, the method has demonstrated its ability to accurately model the interface motion and capture the interface location when the mass transfer rate is uniform and constant. The ability to accurately compute the interface location on unstructured grids opens the path to numerical simulations of liquid-vapor phase change on complex geometries.

Growth of a static bubble in a superheated liquid

The next step is the validation of the coupling between the Navier-Stokes equation and the heat equation by means of the mass transfer rate computed using Eq. ( 56) and the methodology thereafter. This test-case consists in a vapor bubble surrounded by a superheated liquid [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. The temperature field is initialized in the vapor phase at saturation temperature, and is initialized in the liquid phase by means of a radially symmetric profile whose analytical formulation has been derived in three dimensions by Scriven [START_REF] Scriven | On the dynamics of phase growth[END_REF]. Recently we have derived an analytical formulation in two dimensions for the same problem [START_REF] Sahut | Simulation numérique de l'ébullition sur maillages non structurés[END_REF]. To the best of our knowledge, such two-dimensional derivation did not exist in the literature. The mathematical demonstrations of both three-and two-dimensional analytical formulations are detailed in [START_REF] Sahut | Simulation numérique de l'ébullition sur maillages non structurés[END_REF] and recalled below. The temperature gradient is non-zero in the liquid side of the interface, leading to a non-zero mass transfer rate at the interface. Due to the initially stationary bubble and to the symmetry of the problem, the liquid-vapor interface motion is then only due to phase change. The gravitational forces being neglected, the interface is expected to remain spherical in three dimensions and circular in two dimensions for the whole temporal evolution. The duration of the simulation is equal to the theoretical duration needed for the bubble radius to double.

Two-dimensional case

In order to avoid the influence of anisotropic control volumes, we first present the results obtained on two-dimensional grids built only with equilateral triangles, leading to regular hexagonal control volumes. The computational domain is a rhombus of side length L = 2.4 × 10 -2 m of which a portion is shown in Fig. 7. The grids used are listed in Table 2. The initial bubble radius is R 0 = 1 × 10 -3 m and the simulation is performed for the physical time prescribed by the analytical solution for doubling the radius. The

ξ ∞ (R) ξ 2 (R) Slope 1 c
Figure 6: Three-dimensional extension of the simulations illustrated in Fig. 5 to tetrahedral grids. On a, for the finest grid (∆x = 5 × 10 -5 m), the initial interface is shown in blue, the theoretical and computed interfaces are shown at final time in green and yellow, respectively. On b, on the same grid, the complete computed interface is shown at final time. On c, the relative errors on the bubble radius, ξ 2 (R) and ξ∞(R), with grid refinement are shown at final time.

theoretical liquid temperature field is given at any time t by [START_REF] Sahut | Simulation numérique de l'ébullition sur maillages non structurés[END_REF] where s is the reduced variable defined as

T (s) -T ∞ T sat -T ∞ = exp β 2 -s 2 , (70)
s = r 2 √ α liq t , (71)
r = x 2 + y 2 1 2
is the radial variable, α liq = λ liq / (ρ liq c p,liq ) is the thermal diffusivity of the liquid phase, T ∞ is the temperature in the liquid phase when r → ∞, and β is the growth rate defined as

β = 1 - ρ vap ρ liq -1 2 . ( 72
)
Equation ( 70) is used to initialize the liquid temperature field. The theoretical bubble radius R th (t) is given, for some initial t 0 > 0 which will be computed later, at any time

t ≥ t 0 by R th (t) = 2β α liq t. (73)
The physical properties of interest corresponding to a water vapor bubble in liquid water are listed in Table 3 where the heat capacity at constant pressure of the liquid phase has been decreased by two orders of magnitude with respect to the value used in Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] (see discussion in Section 5.2.2). Figure 8 shows the computed interface at final time and the final temperature, velocity and pressure fields for the finest grid. One can notice the high accuracy on the bubble shape at final time : the bubble is perfectly circular and its radius is in excellent agreement with the theoretical one. Figure 9 shows the relative errors ξ 2 (R) and ξ ∞ (R).

ρ liq 958 kg m -3 ρ vap 0.59 kg m -3 σ 5.9 × 10 -2 N m -1 µ liq 2.82 × 10 -4 kg m -1 s -1 µ vap 1.23 × 10 -6 kg m -1 s -1 λ liq 0.6 W m -1 K -1 λ vap 0.026 W m -1 K -1 c p,liq 42.16 J kg -1 K -1 c p,vap 2034 J kg -1 K -1 L v 2.257 × 10 6 J kg -1 T sat 373 K

Three-dimensional case

Our numerical tests are now extended to three dimensions using tetrahedral unstructured grids. Note that this test-case is similar to the test-cases performed by Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] on two-dimensional axisymmetric cartesian grids. One noticeable difficulty arising with such grids is the anisotropy of control volumes. The computational domain is a cube of side length L = 1.2 × 10 -2 m. The initial bubble radius is R 0 = 1 × 10 -3 m and the simulation is performed for the physical time prescribed by the analytical solution for doubling the radius. The liquid temperature field is initialized by [START_REF] Scriven | On the dynamics of phase growth[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Sahut | Simulation numérique de l'ébullition sur maillages non structurés[END_REF]

T (s) -T ∞ T sat -T ∞ = 2β 3 Ja exp β 2 + 2 β 2 +∞ s x -2 exp -x 2 -2 β 3 x -1 dx, ( 74
)
where Ja is the Jakob number defined as

Ja = ρ liq c p,liq (T ∞ -T sat ) ρ vap L v , (75)
β is implicitly defined as

Ja = 2β 3 exp β 2 + 2 β 2 +∞ β x -2 exp -x 2 -2 β 3 x -1 dx, (76)
and = 1 -ρ vap /ρ liq . The theoretical bubble radius is given by Eq. (73). In Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], the simulations are performed for Jakob numbers ranging from 3 to 10. The authors have simulated the test-case using different numerical methods. We chose to compare our method to the variant giving the smallest final relative errors on the bubble radius for the highest tested Jakob numbers, since increasing the Jakob number decreases the width of the thermal boundary layer at the interface, thus making the mass transfer rate more challenging to compute on a given grid [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. We performed our simulations with Jakob numbers of 3, 5, 7 and 9. Table 4 reproduces the errors given for these Jakob numbers in Table 6 of [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. We performed our simulations on the three three-dimensional unstructured tetrahedral grids listed in Table 5. In order to speed up computations, the characteristic cell sizes of these grids are multiplied by a factor 1.6 with respect to the values used in [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] and listed in Table 4. Moreover, the liquid heat capacity at constant pressure is decreased by two orders of magnitude with respect to the value used in [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] (see Table 3), with the goal of decreasing the computing time. Indeed, Eq. (75) shows that if c p,liq decreases, then T ∞ increases (all other parameters being unchanged). As a result, the liquid temperature field, given at time t by Eq. ( 74), has a steeper initial slope close to the interface. Consequently, the mass transfer rate ṁ and so the interface 25 Relative errors (-)

ξ ∞ (R) ξ 2 (R)
Figure 9: Relative errors on the final bubble radius, ξ 2 (R) and ξ∞(R), for the two-dimensional simulations on unstructured grids illustrated in Fig. 8.

∆x (m) Ja=3 Ja=5 Ja=7 Ja=9 6 of Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. The cell sizes in the first column have been computed from the number of cells and the size of the domain given by the authors.

velocity due to phase change (right-hand side of Eq. ( 30)) are then increased, leading to faster simulations. Table 6 shows ξ ∞ (R) on the three grids and for the four Jakob numbers considered. The relative error ξ ∞ (R) decreases with grid refinement for all Jakob numbers except for Ja = 3 on the finest grid. One can indeed observe an error increase for the smallest Jakob number on the finest grid. This behavior can be related to the results of Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] listed in Table 4 in which the error on the bubble radius reaches a minimum value on the finest grid and does not clearly converge with grid refinement below this minimum value. Also, ξ ∞ (R) increases with the Jakob number.

As stated above and in Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF], increasing the Jakob number reduces, by Eq. ( 74), the width of the thermal boundary layer at the interface. Consequently, the mass transfer rate and the interface velocity due to phase change are computed with less accuracy since a lower number of grid nodes are located in the thermal boundary layer at every iteration. Then, the bubble radius, measured at final time from the signed distance function advected by the interface velocity due to phase change, is also less accurate. The convergence of the bubble radius with grid refinement and decrease of the Jakob number has been obtained by Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] on two-dimensional axisymmetric cartesian grids, as summarized in Table 4. In the present work, this convergence is extended to ∆x (m) Number of cells N (-)

1.5 × 10 -4 554 800 7.5 × 10 -5

1 933 621 3.75 × 10 -5 9 186 699

Table 5: Three-dimensional tetrahedral grids. The cell sizes correspond to the ones listed in Table 4 multiplied by 1.6. The refined area around the interface is shown in Fig. 10.

∆x (m) Ja=3 Ja=5 Ja=7 Ja=9

1.5 × 10 -4 6.8% 20% 34% 47% 7.5 × 10 -5 0.8% 1.7% 7.9% 16% 3.75 × 10 -5 1.7% 0.6% 1.0% 1.9% three-dimensional unstructured grids with comparable accuracy: in our simulations, the bubble radius at final time is generally slightly more accurate for low Jakob numbers (3 and 5) slightly less accurate for high Jakob numbers (7 and 9) than the one obtained by Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF]. Since our grid cells are larger than the ones used in Tanguy et al. [START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF] by a factor 1.6, and since our grids are three-dimensional, which implies that the sizes of the tetrahedra mentioned in Table 5 are only approximated by the meshing software, further comparison would not be relevant. Table 6 is then considered as the validation of our solver for the growth of static bubbles.

For a Jakob number of 5, Eq. (76) gives β = 5.304, Eq. (73) gives the initial time t 0 = 5.98 × 10 -4 s and the final time t 1 = 2.39 × 10 -3 s and Eq. (75) gives T ∞ = 538 K. For illustrative purposes, Fig. 10 shows a portion of the coarsest grid used to validate this test-case in three dimensions. Figure 11 shows the initial interface and the final computed and theoretical interfaces for the finest grid, as well as the relative errors on the final bubble radius, ξ 2 (R) and ξ ∞ (R), for the three grids used. The final bubble radius converges with grid refinement for both relative errors, ξ 2 (R) and ξ ∞ (R). We emphasize here the excellent agreement between the theoretical and numerical results for the finest grid, whereas, for ξ ∞ (R), one single inaccurate nodal value can lead to a high relative error. This is a very stringent criterion constituting a severe test for the numerical accuracy. On the finest grid, the relative error ξ ∞ (R) is equal to 0.6%, as reported in Table 6. Figure 12 shows the radial temperature at initialization and final time. As expected, the liquid temperature is very close to T ∞ in all the liquid phase, except in a region close to the interface, the thermal boundary layer, where a steep gradient responsible for the interface movement is observed. In this area, the liquid temperature decreases from T ∞ = 538 K in the liquid phase to T sat = 373 K at the interface, showing that the immersed Dirichlet boundary condition given by Eq. ( 22) is satisfied throughout the simulation. Indeed, the use of an implicit algorithm to impose saturation temperature at the subgrid position of the interface does not allow for a direct verification of the respect of such constraint. Figure 13 shows the liquid temperature and velocity fields at final time on the finest grid. One can see the small thickness of the boundary layer, compared to the bubble diameter, which requires a fine enough grid to be captured. As expected, the liquid is ejected from the interface and the liquid velocity field is thus aligned with the interface normal vector. The velocity jump at the interface due to phase change is responsible for the liquid motion which can exit the domain thanks to the outlet boundary conditions used on all six faces of the cubic domain.

Our solver is able to maintain the spherical shape of the bubble throughout the simulation as expected from the theoretical solution. Since this test-case is very severe relatively to the accuracy of the interface capturing or tracking method used, especially on unstructured grids, we believe that the present results have a strong interest for further simulations of two-phase flows with phase change on unstructured grids.

Conclusion and perspectives

In this work, a numerical method for simulations of two-phase flows with phase change due to heat transfer (boiling) on unstructured grids has been developed within the YALES2 code and assessed in two and three dimensions. The method has first been validated on the test-case of a static growing bubble with a fixed uniform and constant mass transfer rate, in two-and three-dimensional unstructured grids. Compared to the analytical solution, the bubble radius at final time is first-order accurate with grid refinement in the L ∞ sense. The ability to accurately compute the interface location on unstructured grids opens the path to numerical simulations of liquid-vapor phase change 8 to tetrahedral grids. On a, for the finest grid (∆x = 3.75 × 10 -5 m), the initial interface is shown in blue, the theoretical and computed interfaces are shown at final time in green and yellow, respectively. On b, on the same grid, the complete computed interface is shown at final time. On c, the relative errors on the bubble radius, ξ 2 (R) and ξ∞(R), with grid refinement are shown at final time. Local errors in the interface position, as shown in red on b, strongly affect the value of ξ∞(R) whereas their effect on ξ 2 (R) is more limited, as shown in c. on complex geometries. The method has then been extended to take into account the computation of the mass transfer rate at the interface from the thermal fluxes on both sides of the interface. To this purpose, the heat equation is solved in both phases and two temperature fields are used to take into account the discontinuity of the heat flux at the interface. An immersed Dirichlet boundary condition is imposed at the interface in order to ensure that boiling always occurs at saturation temperature. The mass transfer rate is computed using a new framework to extrapolate the liquid temperature gradient from the liquid phase to the interface, and the vapor temperature gradient from the vapor phase to the interface, in order to compute the mass transfer rate precisely at the interface. The overall implementation has been validated against the test-case of a bubble at rest growing in a superheated liquid, in two and three dimensions. Excellent agreement with the theoretical bubble radius has been obtained at final time on twoand three-dimensional unstructured grids. Moreover, the bubble radius converges with grid refinement in all cases. In three dimensions, the relative L ∞ -norm of the error on the bubble radius is below 2% on the finest grid considered. These results demonstrate the ability of our numerical method to accurately simulate two-phase flows with phase change where the mass transfer rate is computed from the thermal fluxes at the interface, thus taking a step towards realistic boiling simulations in industrial context. To the best of our knowledge, this work has provided the first numerical simulations of bubble growth in three dimensions with a bubble radius converging with grid refinement in the L ∞ -norm sense. The L ∞ convergence of the radius implies that the position of the bubble surface converges to the exact solution at every surface point. It is then considered as an important ingredient for predictive numerical simulations.

The natural extension of this work is the simulation of nucleate boiling, a mode of heat transfer widely used in industrial applications, occurring when a liquid is in contact with a solid whose temperature is above the liquid boiling point, leading to the formation of vapor bubbles on the solid surface. Numerical simulations of nucleate boiling require a methodology to take the motion of the contact line (where the solid, liquid and vapor phases meet) into account, as well as the contact angle existing between the solid surface and the liquid-vapor interface. The numerical simulation of nucleate boiling on unstructured grids is currently under development as an extension of this work. The inclusion of the contact line and contact angle in our method is expected to enable predictive three-dimensional direct numerical simulations of nucleate boiling at the bubble scale in the coming years.

Figure 1 :

 1 Figure 1: The volume V is a thin band of width δ and length δ centered on the interface Γ moving at velocity W = Dn, and delimited by the surface S = S 1 ∪ • • • ∪ S 4 . The point P is at rest with respect to the interface.

ρP n S dS + 2 S

 2 (u -Dn) {(u -2Dn) • n S } dS -S µSn S dS + Γ σκn dγ.

Figure 2 :

 2 Figure 2: Control volume Vp associated to node p. The neighboring nodes of p are numbered counterclockwise. The barycenter of each element and the middle of each pair of nodes are respectively represented by the symbols + and •. The control volume associated to node p is represented in simple black dashes.To each pair of nodes are associated two facets of the control volume, each linking the barycenter of a neighboring element to the middle of the pair. The sum of the fluxes of a physical field at the two facets is equal to the flux of the same physical field at the "face" linking the two barycenters of the neighboring elements to the pair, represented in dense red dashes. If this face is denoted S i and if n i is its outward normal vector to the control volume of p, then the vector Ap,q i represented in red is equal to S i n i . For more clarity, S i and n i are not represented. The interface Γ is represented in blue.

) 2 imposing 2 to obtain u n+1 liq and u n+1 vap Eq. ( 33 )

 2233 Compute heat equation advection terms for T n+1 liq and T n+1 vap Eq. (57) Solve linear systems for T n+1 liq and T n+1 vap with implicit diffusion imposing TΓ = Tsat Eq. (64) Perform constant extrapolation of u * liq and u * vap across the interface along n n+1 Eq. (52) Solve Poisson equation to compute P n+ 1Eqs. (32), (54), (18), (55) and (6)Correct u * liq and u * vap with ∇P n+ 1Perform constant extrapolation of u n+1 liq and u n+1 vap across the interface along n n+1

Figure 4 :

 4 Figure 4: Coarsest grid used in two dimensions (∆x = 4 × 10 -4 m) for the simulation of the growth of a static bubble with a constant and uniform mass transfer rate. The small circle represents the initial interface of radius R 0 = 1 × 10 -3 m, whereas the bigger circle represents the final theoretical interface (the radius is twice the initial one).

Figure 5 :

 5 Figure 5: Two-dimensional simulations on triangular grids. On a, the interface location is shown at initial and final times on the finest grid considered (∆x = 5 × 10 -5 m). The initial interface is represented in blue, the computed interface at final time, in black, and the theoretical interface at final time, in red. The liquid and vapor velocity fields are also shown at final time. On b, the relative errors, ξ 2 (R) and ξ∞(R), with grid refinement are shown at final time.

Figure 7 :

 7 Figure 7: Coarsest two-dimensional grid (∆x = 9.375 × 10 -5 m) composed only of equilateral triangles used is our simulations. The small circle represents the initial interface of radius R 0 = 1 × 10 -3 m, whereas the bigger circle represents the final theoretical interface (the radius is twice the initial one).

Figure 8 :

 8 Figure 8: Results in two dimensions on the finest triangular grid used (∆x = 2.344 × 10 -5 m). On a, the initial interface is shown in blue, the final computed interface, in black, and the final theoretical interface in red. The final liquid and vapor temperature fields are shown on b, the final velocity fields are shown on c and the final pressure field is shown on d.

Figure 10 :

 10 Figure10: Refined three-dimensional cubic grid used in our simulations. The small circle represents the initial interface location of radius R 0 = 1 × 10 -3 m, the bigger circle represents the final theoretical interface location (the radius is twice the initial one). The grid is refined only in an area englobing a narrow band around the interface in order to speed up computations. This grid has a characteristic cell size around the interface of ∆x = 1.5 × 10 -4 m and counts 554 800 grid cells. This is the coarsest grid used in our 3D tests.

Figure 11 :

 11 Figure11: Three-dimensional extension of the simulations illustrated in Fig.8to tetrahedral grids. On a, for the finest grid (∆x = 3.75 × 10 -5 m), the initial interface is shown in blue, the theoretical and computed interfaces are shown at final time in green and yellow, respectively. On b, on the same grid, the complete computed interface is shown at final time. On c, the relative errors on the bubble radius, ξ 2 (R) and ξ∞(R), with grid refinement are shown at final time. Local errors in the interface position, as shown in red on b, strongly affect the value of ξ∞(R) whereas their effect on ξ 2 (R) is more limited, as shown in c.

Figure 12 :Figure 13 :

 1213 Figure 12: Radial profile of the three-dimensional temperature field shown at a initialization (R 0 = 1 × 10 -3 m) and b final time (R 1 = 2 × 10 -3 m) for ∆x = 3.75 × 10 -5 m.

  ) at node x p is given by

	  ρ liq c p,liq ∆t	+	1 V p	λ liq	∆x p,q1 θ ∆x p,q1	2 • A p,q1 +	1 V p	5 j=2	λ liq	∆x p,q1 ∆x p,qj	2 • A p,qj	  T liq,n+1 p
										-	1 V p	5 j=2	λ liq	∆x p,qj ∆x p,qj	2 • A p,qj T liq,n+1 qj
	=	ρ liq c p,liq ∆t	T liq,n p	-ρ liq c p,liq u liq p • ∇T | liq,n p	+	1 V p	λ liq	∆x p,q1

Table 1 :

 1 Physical properties used in Section 5.1.

Table 2 :

 2 Two-dimensional triangular grids.

Table 3 :

 3 Physical properties used in Section 5.2. Only c p,liq differs (by two orders of magnitude) from the value used in Tanguy et al.[START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF].

Table 4 :

 4 Relative errors on the bubble radius at final time for different Jakob numbers on three twodimensional axi-symmetric cartesian grids extracted from Table

Table 6 :

 6 Relative error on the bubble radius at final time, ξ∞(R), for different unstructured grids and Jakob numbers in three dimensions.

Acknowledgments

This work has been supported by the LabEx Tec21 (Investissements d'Avenir-Grant Agreement No. ANR-11-LABX-0030). The authors would like to thank Patrick Bégou for support on algorithms implementation, Geoffroy Vaudor and Romain Janodet for their help on the use of the geometrical approach for the reinitialization of the Level Set function, Bruno Voisin for his help on some mathematical derivations, and Savinien Pertant for various fruitful discussions on the numerical method and its implementation.