# A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows

Leandro Germes Martinez, Benjamin Duret, Julien Reveillon, Francois-Xavier Demoulin

► To cite this version:

Leandro Germes Martinez, Benjamin Duret, Julien Reveillon, Francois-Xavier Demoulin. A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows. International Conference of Liquid Atomization and Sprays Systems (ICLASS 2021), Aug 2021, Edimbourg, United Kingdom. pp.5847, 10.2218/iclass.2021.5847. hal-03372960

HAL Id: hal-03372960

https://hal.science/hal-03372960

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows

Leandro Germes Martinez, Benjamin Duret, Julien Reveillon, Francois-Xavier Demoulin

► To cite this version:

Leandro Germes Martinez, Benjamin Duret, Julien Reveillon, Francois-Xavier Demoulin. A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows. International Conference of Liquid Atomization and Sprays Systems (ICLASS 2021), Aug 2021, Edimbourg, United Kingdom. pp.5847, 10.2218/iclass.2021.5847. hal-03372960

HAL Id: hal-03372960

https://hal.science/hal-03372960v1

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows

Leandro Germes Martínez*, Benjamin Duret, Julien Reveillon, François-Xavier Demoulin

CORIA-UMR 6614 - Normandie Université, CNRS-Université et INSA de Rouen, Campus 675, avenue de l'Université, Saint Etienne du Rouvray 76800, France

*Corresponding author email: germesml@coria.fr

## Abstract

This work intent to provide numerical tools, allowing simulations of two-phase flow incorporating compressible effects and a proper treatment of the jump conditions found at the interface due to large density ratios and phase change. For achieving this task, first, the incompressible level-set method for vaporizing two-phase flows proposed by Tanguy et al. [12] is revisited and adapted to a conservative interface representation: The Coupled Level-set/Volume of Fluid method. In this context, results are presented for simulations of the evaporation of an isolated static droplet. Then, the formalism is extended into a compressible one using as framework the pressure-based method proposed by Duret et al. [4]. For validation, a 3D simulation of a static evaporating droplet is studied for showing the proper mass balance between the two phases. Finally, a 3D two-phase compressible Homogenous Isotropic Turbulence (HIT) configuration is presented to demonstrate the potential of this method in presence of breakup, coalescence and evaporation processes.

## Keywords

CLSVOF; Turbulent flows; Phase change; Ghost Fluid Method; Compressibility

## Introduction

Liquid-gas flows are part of large number of engineering systems, e.g., internal combustion engines, heat exchangers, multiphase reactors, etc. Since heat and mass transfers are mainly driven by the process of phase change, its understanding is an important part in the design and optimization of this type of systems. The majority of the existing simulations devoted to the aforementioned systems use an incompressible numerical formalism. Such assumption is not able to consider the temporal and spatial density gradients inside gas structures such as bubbles enclosed in an evaporating liquid. Consequently, outflow boundary conditions are required to allow the velocity field to be divergence-free in each phase.

When dealing with the simulation of two-phase flows with phase change, numerous difficulties appears due to the diversity of physical processes that takes place. Being the most challenging the discontinuity created at the interface in the velocity field. Most numerical methods dedicated to the convective transport of the scalars in DNS of two-phase flows are originally developed for handling divergence-free velocity fields. Hence, they cannot use this velocity field directly. A solution has been proposed by [6] and then modified by [9]. The idea of this method is to compute a divergence-free liquid-velocity field for advecting of the liquid volume fraction, by subtracting the velocity jump at the interface from the original velocity field. Another method has been proposed by [8]. Similar to [9], its purpose is to construct a divergence-free velocity field to advect the liquid volume fraction. In this case, [8] takes the original velocity field and extend the velocity of the liquid phase into the gas phase using the Aslam's extension method [1]. A third method is the one proposed by [12], here, the ghost fluid method is used for extending both the liquid and gas velocity fields, and this way obtaining a continuous liquid velocity field. The main objective of this paper is to present numerical strategies for simulating two-phase flows which take into account the surface tension, compressible effects, phase change and
# A weakly-compressible DNS formalism dedicated to evaporating turbulent liquid-gas flows

Leandro Germes Martínez*, Benjamin Duret, Julien Reveillon, François-Xavier Demoulin

CORIA-UMR 6614 - Normandie Université, CNRS-Université et INSA de Rouen, Campus 675, avenue de l'Université, Saint Etienne du Rouvray 76800, France

*Corresponding author email: germesml@coria.fr

## Abstract

This work intent to provide numerical tools, allowing simulations of two-phase flow incorporating compressible effects and a proper treatment of the jump conditions found at the interface due to large density ratios and phase change. For achieving this task, first, the incompressible level-set method for vaporizing two-phase flows proposed by Tanguy et al. [12] is revisited and adapted to a conservative interface representation: The Coupled Level-set/Volume of Fluid method. In this context, results are presented for simulations of the evaporation of an isolated static droplet. Then, the formalism is extended into a compressible one using as framework the pressure-based method proposed by Duret et al. [4]. For validation, a 3D simulation of a static evaporating droplet is studied for showing the proper mass balance between the two phases. Finally, a 3D two-phase compressible Homogenous Isotropic Turbulence (HIT) configuration is presented to demonstrate the potential of this method in presence of breakup, coalescence and evaporation processes.

## Keywords

CLSVOF; Turbulent flows; Phase change; Ghost Fluid Method; Compressibility

## Introduction

Liquid-gas flows are part of large number of engineering systems, e.g., internal combustion engines, heat exchangers, multiphase reactors, etc. Since heat and mass transfers are mainly driven by the process of phase change, its understanding is an important part in the design and optimization of this type of systems. The majority of the existing simulations devoted to the aforementioned systems use an incompressible numerical formalism. Such assumption is not able to consider the temporal and spatial density gradients inside gas structures such as bubbles enclosed in an evaporating liquid. Consequently, outflow boundary conditions are required to allow the velocity field to be divergence-free in each phase.

When dealing with the simulation of two-phase flows with phase change, numerous difficulties appears due to the diversity of physical processes that takes place. Being the most challenging the discontinuity created at the interface in the velocity field. Most numerical methods dedicated to the convective transport of the scalars in DNS of two-phase flows are originally developed for handling divergence-free velocity fields. Hence, they cannot use this velocity field directly. A solution has been proposed by [6] and then modified by [9]. The idea of this method is to compute a divergence-free liquid-velocity field for advecting of the liquid volume fraction, by subtracting the velocity jump at the interface from the original velocity field. Another method has been proposed by [8]. Similar to [9], its purpose is to construct a divergence-free velocity field to advect the liquid volume fraction. In this case, [8] takes the original velocity field and extend the velocity of the liquid phase into the gas phase using the Aslam's extension method [1]. A third method is the one proposed by [12], here, the ghost fluid method is used for extending both the liquid and gas velocity fields, and this way obtaining a continuous liquid velocity field. The main objective of this paper is to present numerical strategies for simulating two-phase flows which take into account the surface tension, compressible effects, phase change and
large density ratio in a turbulent flow. We propose to achieve this task by extending the projection method proposed by [12] into a compressible formalism. Also by combining it with an accurate and conservative interface capturing method (CLSVOF method), allowing to simulate atomization processes (high Weber and Reynolds numbers).

In the following part of this work, first, the general formalism for a compressible two-phase flow, the interface jump conditions and the numerical strategies are described. Then, in order to validate the thermodynamics of our incompressible formalism, numerical results are compared with the D² law. Finally, a 3D HIT configuration is presented to show the capacity of the method to handle compressible turbulent two-phase flows with phase change.

## Governing Equations

In this work, a system of two immiscible fluids is studied; composed by a mono-component liquid and an ideal mixture of gases (e.g. air and vaporized liquid). Mass, energy and momentum are transported between the two phases through an interface while supposing thermodynamic equilibrium.

The continuity and momentum equations can be written as:

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) &= 0 \\ \frac{\partial \rho \vec{u}}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u} \otimes \vec{u}) &= -\vec{\nabla} P + \vec{\nabla} \cdot (2\mu\bar{\epsilon} - \frac{2}{3}\mu(\vec{\nabla} \cdot \vec{u})\bar{I}) + \rho \vec{f} \end{aligned} \right. \quad (1) $$

where $\bar{\epsilon}$ is the strain rate tensor, defined as: $\bar{\epsilon} = \frac{1}{2}(\vec{\nabla}\vec{u} + \vec{\nabla}\vec{u}^T)$

When both fluids are considered to be compressible, the system of equations 1 must be closed by two equations of state. For the liquid, the Tait's equation is considered. Concerning the gas phase, the equation of state for ideal gases is used. The evaporation depends mostly of the temperature and the vapor concentration at the interface. Hence, the need of solving the following conservation equations. The energy conservation equation can be expressed in term of temperature:

$$ \rho c_p \left( \frac{\partial T}{\partial t} + \vec{u} \cdot \vec{\nabla} T \right) = \vec{\nabla} \cdot \left( \lambda \vec{\nabla} T \right) + \alpha_T T \left( \frac{\partial P}{\partial t} + \vec{u} \cdot \vec{\nabla} P \right) + \bar{\tau} \vec{\nabla} \vec{u} + \dot{Q} \quad (2) $$

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

$$ \frac{\partial \rho_g Y_v}{\partial t} + \vec{\nabla} \cdot (\rho_g \vec{u} Y_v) = \vec{\nabla} \cdot (\rho_g D_m \vec{\nabla} Y_v) \quad (3) $$

## Interface jump conditions

Considering that we are working with a one-fluid formalism, appropriate jump conditions must be defined for the discontinuous variables at the interface ($\rho$, $\mu$ and $\lambda$). When mass transfer is considered across the interface, additional jumps conditions are required to fully define our system. The mass conservation across the interface can be expressed as [3]:

$$ \rho_l (\vec{u}_l^\Gamma - \vec{u}^\Gamma) \cdot \vec{n} = \rho_g (\vec{u}_g^\Gamma - \vec{u}^\Gamma) \cdot \vec{n} = \dot{\omega} \quad (4) $$

where $\dot{\omega}$ represent the mass flux across the interface due to evaporation (i.e. the vaporization rate). $\vec{u}^\Gamma$ is the velocity of the interface; which can be described, from the liquid side, as the sum of the liquid velocity at the interface $(\vec{u}_l^\Gamma)$ and the interface regression speed $(\vec{s}_d)$:

$\vec{u}^\Gamma = \vec{u}_l^\Gamma - \|\vec{s}_d\|\vec{n}$. $\vec{s}_d$ is collinear to the normal of the interface and defined as: $\vec{s}_d = -\frac{\dot{\omega}}{\rho_l}\vec{n}$

Additionally, in the presence of phase change, an extra energy jump condition must be considered to respect the energy conservation equation:

$$ \lambda_l \vec{\nabla} T_l^\Gamma \cdot \vec{n} - \lambda_g \vec{\nabla} T_g^\Gamma \cdot \vec{n} = -h_{lv} \dot{\omega} \quad (5) $$
large density ratio in a turbulent flow. We propose to achieve this task by extending the projection method proposed by [12] into a compressible formalism. Also by combining it with an accurate and conservative interface capturing method (CLSVOF method), allowing to simulate atomization processes (high Weber and Reynolds numbers).

In the following part of this work, first, the general formalism for a compressible two-phase flow, the interface jump conditions and the numerical strategies are described. Then, in order to validate the thermodynamics of our incompressible formalism, numerical results are compared with the D² law. Finally, a 3D HIT configuration is presented to show the capacity of the method to handle compressible turbulent two-phase flows with phase change.

## Governing Equations

In this work, a system of two immiscible fluids is studied; composed by a mono-component liquid and an ideal mixture of gases (e.g. air and vaporized liquid). Mass, energy and momentum are transported between the two phases through an interface while supposing thermodynamic equilibrium.

The continuity and momentum equations can be written as:

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u}) &= 0 \\ \frac{\partial \rho \vec{u}}{\partial t} + \vec{\nabla} \cdot (\rho \vec{u} \otimes \vec{u}) &= -\vec{\nabla} P + \vec{\nabla} \cdot \left( 2\mu\bar{\epsilon} - \frac{2}{3}\mu (\vec{\nabla} \cdot \vec{u}) \bar{I} \right) + \rho \vec{f} \end{aligned} \right. \quad (1) $$

where $\bar{\epsilon}$ is the strain rate tensor, defined as: $\bar{\epsilon} = \frac{1}{2} (\vec{\nabla}\vec{u} + \vec{\nabla}\vec{u}^T)$

When both fluids are considered to be compressible, the system of equations 1 must be closed by two equations of state. For the liquid, the Tait's equation is considered. Concerning the gas phase, the equation of state for ideal gases is used. The evaporation depends mostly of the temperature and the vapor concentration at the interface. Hence, the need of solving the following conservation equations. The energy conservation equation can be expressed in term of temperature:

$$ \rho c_p \left( \frac{\partial T}{\partial t} + \vec{u} \cdot \vec{\nabla} T \right) = \vec{\nabla} \cdot \left( \lambda \vec{\nabla} T \right) + \alpha_T T \left( \frac{\partial P}{\partial t} + \vec{u} \cdot \vec{\nabla} P \right) + \bar{\tau} \vec{\nabla} \vec{u} + \dot{Q} \quad (2) $$

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

$$ \frac{\partial \rho_g Y_v}{\partial t} + \vec{\nabla} \cdot (\rho_g \vec{u} Y_v) = \vec{\nabla} \cdot (\rho_g D_m \vec{\nabla} Y_v) \quad (3) $$

## Interface jump conditions

Considering that we are working with a one-fluid formalism, appropriate jump conditions must be defined for the discontinuous variables at the interface ($\rho$, $\mu$ and $\lambda$). When mass transfer is considered across the interface, additional jumps conditions are required to fully define our system. The mass conservation across the interface can be expressed as [3]:

$$ \rho_l (\vec{u}_l^\Gamma - \vec{u}^\Gamma) \cdot \vec{n} = \rho_g (\vec{u}_g^\Gamma - \vec{u}^\Gamma) \cdot \vec{n} = \dot{\omega} \quad (4) $$

where $\dot{\omega}$ represent the mass flux across the interface due to evaporation (i.e. the vaporization rate). $\vec{u}^\Gamma$ is the velocity of the interface; which can be described, from the liquid side, as the sum of the liquid velocity at the interface ($\vec{u}_l^\Gamma$) and the interface regression speed ($\vec{s}_d$):

$\vec{u}^\Gamma = \vec{u}_l^\Gamma - \|\vec{s}_d\|\vec{n}$. $\vec{s}_d$ is collinear to the normal of the interface and defined as: $\vec{s}_d = -\frac{\dot{\omega}}{\rho_l}\vec{n}$

Additionally, in the presence of phase change, an extra energy jump condition must be considered to respect the energy conservation equation:

$$ \lambda_l \vec{\nabla} T_l^\Gamma \cdot \vec{n} - \lambda_g \vec{\nabla} T_g^\Gamma \cdot \vec{n} = -h_{lv} \dot{\omega} \quad (5) $$
Finally, to respect the species conservation for a monocomponent liquid, the jump condition at the interface for $Y_v$ is:

$$ \rho_g D_m \vec{\nabla} Y_v \cdot \vec{n}|_g^\Gamma = \dot{\omega}(1 - Y_{vs}) \quad (6) $$

where $Y_{vs}$ is the vapor mass fraction at the interface in the gaseous phase. In order to obtain it, thermodynamic equilibrium is assumed at the interface; which requires a direct relation between the energy and $Y_v$. To estimate $Y_{vs}$, the Clausius-Clapeyron relation is used.

The interface jump conditions of the different variables reveal the strong coupling between energy and the velocity jump. Handling this coupling is challenging and require robust and accurate numerical methods dedicated to DNS of two-phase flow with phase change.

## Numerical Methods

### Interface capturing method

The Coupled Level-set/Volume of Fluid (CLSVOF) method implemented in the ARCHER code [7] is extended to handle multiphase flows with phase change. Here, two additional terms have been added to the liquid volume fraction equation; a term due to the compressibility and other for the evaporation of the liquid:

$$ \frac{\partial \alpha_l}{\partial t} + \vec{\nabla} \cdot (\alpha_l \vec{u}_l) = - \frac{\alpha_l D \rho_l}{\rho_l Dt} - \frac{\dot{m}}{\rho_l} \quad (7) $$

where $\dot{m}$ represent the quantity of liquid volume evaporated in each cell containing an interface, and it is estimated by:

$$ \dot{m} = \rho_l \Sigma \| \vec{s}_d \| \quad (8) $$

where $\Sigma = \frac{S}{V_{cell}}$ is the interfacial surface density of a mixed cell, $S$ is the surface, calculated based in the PLIC reconstruction. The material derivative of the density is computed as:

$$ \frac{D \rho_l}{Dt} = \frac{1}{c_l^2} \frac{DP}{Dt} - \frac{\alpha_{Tl}}{\rho_l c_{pl}} \vec{\nabla} \cdot (\lambda \vec{\nabla} T) \quad (9) $$

The pressure and temperature term of eq. 9 are solved after the computation of the Poisson equation for the pressure and during the evolution of the temperature, respectively. Both of them will be explained in the next sections.

To couple the VoF method and the LS method it is necessary to solve both of their transport equations simultaneously. The discretization of both equations is made following the "coupled" second-order conservative operator split advection scheme defined in [11].

The signed distance function is used to obtain the geometrical information of the interface, e.g.,

the vector normal to the interface $(\vec{n} = \frac{\vec{\nabla}\phi}{|\vec{\nabla}\phi|})$ and the curvature $(\kappa = -\vec{\nabla} \cdot \vec{n})$ needed for the computation of the surface tension force. A redistancing algorithm is applied every iteration step following the same principle of previous work ([12], [5]).

### Flow solver

A projection method is used for the temporal resolution of the momentum equation. An Eulerian staggered grid, where the scalars variables ($P, T, \alpha_l, ...$) are defined on the cell center, and the velocity is computed on the faces of the cell.

In this work, we are going to use the method proposed by [12] for a reacting interface. The Ghost Fluid Method (GFM) is used to populate the ghost cells on each side of the interface. And, the velocity extensions allow to have an liquid/gas velocity fields which are continuous at the interface. For the compressible part of this work, the projection method of [12] had
Finally, to respect the species conservation for a monocomponent liquid, the jump condition at the interface for $Y_v$ is:

$$ \rho_g D_m \vec{\nabla} Y_v \cdot \vec{n}|_g^\Gamma = \dot{\omega}(1 - Y_{vs}) \quad (6) $$

where $Y_{vs}$ is the vapor mass fraction at the interface in the gaseous phase. In order to obtain it, thermodynamic equilibrium is assumed at the interface; which requires a direct relation between the energy and $Y_v$. To estimate $Y_{vs}$, the Clausius-Clapeyron relation is used.

The interface jump conditions of the different variables reveal the strong coupling between energy and the velocity jump. Handling this coupling is challenging and require robust and accurate numerical methods dedicated to DNS of two-phase flow with phase change.

## Numerical Methods

### Interface capturing method

The Coupled Level-set/Volume of Fluid (CLSVOF) method implemented in the ARCHER code [7] is extended to handle multiphase flows with phase change. Here, two additional terms have been added to the liquid volume fraction equation; a term due to the compressibility and other for the evaporation of the liquid:

$$ \frac{\partial \alpha_l}{\partial t} + \vec{\nabla} \cdot (\alpha_l \vec{u}_l) = - \frac{\alpha_l D \rho_l}{\rho_l Dt} - \frac{\dot{m}}{\rho_l} \quad (7) $$

where $\dot{m}$ represent the quantity of liquid volume evaporated in each cell containing an interface, and it is estimated by:

$$ \dot{m} = \rho_l \Sigma \| \vec{s}_d \| \quad (8) $$

where $\Sigma = \frac{S}{V_{cell}}$ is the interfacial surface density of a mixed cell, $S$ is the surface, calculated based in the PLIC reconstruction. The material derivative of the density is computed as:

$$ \frac{D \rho_l}{Dt} = \frac{1}{c_l^2} \frac{DP}{Dt} - \frac{\alpha_{Tl}}{\rho_l c_{pl}} \vec{\nabla} \cdot (\lambda \vec{\nabla} T) \quad (9) $$

The pressure and temperature term of eq. 9 are solved after the computation of the Poisson equation for the pressure and during the evolution of the temperature, respectively. Both of them will be explained in the next sections.

To couple the VoF method and the LS method it is necessary to solve both of their transport equations simultaneously. The discretization of both equations is made following the "coupled" second-order conservative operator split advection scheme defined in [11].

The signed distance function is used to obtain the geometrical information of the interface, e.g.,

the vector normal to the interface $(\vec{n} = \frac{\vec{\nabla}\phi}{|\vec{\nabla}\phi|})$ and the curvature $(\kappa = -\vec{\nabla} \cdot \vec{n})$ needed for the computation of the surface tension force. A redistancing algorithm is applied every iteration step following the same principle of previous work ([12], [5]).

### Flow solver

A projection method is used for the temporal resolution of the momentum equation. An Eulerian staggered grid, where the scalars variables ($P, T, \alpha_l, ...$) are defined on the cell center, and the velocity is computed on the faces of the cell.

In this work, we are going to use the method proposed by [12] for a reacting interface. The Ghost Fluid Method (GFM) is used to populate the ghost cells on each side of the interface. And, the velocity extensions allow to have an liquid/gas velocity fields which are continuous at the interface. For the compressible part of this work, the projection method of [12] had
been modified to take into account the compressible effects, i. e., density gradients due to an increase of pressure in an enclosed environment containing an evaporating liquid and spatial variation of the temperature field.

First, an intermediate velocity for both liquid and gas field $(\vec{u}_l^*, \vec{u}_g^*)$ is computed by solving the momentum equation (eq.1) without the pressure:

$$ \left\{ \begin{aligned} \vec{u}_l^* &= \vec{u}_l^m - \Delta t \left( \left( \vec{u}_l^m \cdot \vec{\nabla} \right) \vec{u}_l^m - \vec{\nabla} \cdot \left( 2\mu\bar{\epsilon}_l - \frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}_l^m \bar{I} \right) \right) \\ \vec{u}_g^* &= \vec{u}_g^m - \Delta t \left( \left( \vec{u}_g^m \cdot \vec{\nabla} \right) \vec{u}_g^m - \vec{\nabla} \cdot \left( 2\mu\bar{\epsilon}_g - \frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}_g^m \bar{I} \right) \right) \end{aligned} \right. \quad (10) $$

the term $\frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}^m \bar{I}$ had been added to adapt the projection method to a compressible formalism.

Then, the second term of the RHS for Poisson equation is obtained with:

$$ f = \begin{cases} \vec{\nabla} \cdot \vec{u}_l^* & \text{if } \phi > 0, \\ \vec{\nabla} \cdot \vec{u}_g^* & \text{if } \phi < 0, \end{cases} \quad (11) $$

The Poisson equation is computed with following strategy:

$$ -\vec{\nabla} \left( \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} \right) = \vec{\nabla} \cdot \vec{u}^{n+1} - f \quad (12) $$

where $\vec{\nabla} \cdot \vec{u}^{n+1}$ is computed with:

$$ \vec{\nabla} \cdot \vec{u}^{n+1} = - \left( \frac{\alpha_l}{\rho_l C_l^2} + \frac{\alpha_g}{\rho_g C_g^2} \right) \frac{DP}{Dt} - \left( \frac{\alpha_g \alpha_{Tg}}{\rho_g c_{pg}} + \frac{\alpha_l \alpha_{Tl}}{\rho_l c_{pl}} \right) \vec{\nabla} \cdot (\lambda \vec{\nabla} T) \quad (13) $$

The convective term of the pressure material derivative $(\vec{u} \cdot \vec{\nabla} P)$ is discretized with a 5th order WENO scheme. The Poisson solver used to solve this Helmholtz type equation 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:

$$ \left\{ \begin{aligned} \vec{u}_l^{n+1} &= \vec{u}_l^* - \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} & \text{if } \phi > 0, \\ \vec{u}_g^{n+1} &= \vec{u}_g^* - \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} & \text{if } \phi < 0, \end{aligned} \right. \quad (14) $$

The next steps are dedicated to the computation of the liquid and gas velocity extensions $(\vec{u}_l^{\text{ghost}}, \vec{u}_g^{\text{ghost}})$. For this, a procedure similar to [12] is used.

## Heat and mass transfer description

To follow the evolution of the temperature field, in the energy equation (eq. 2), the term responsible for the energy generation by friction $(\bar{\tau} \vec{\nabla} \vec{u})$ and the term of volumetric energy generation $(\dot{Q})$ are neglected in this work:

$$ \rho c_p \left( \frac{\partial T}{\partial t} + \vec{u} \cdot \vec{\nabla} T \right) = \vec{\nabla} \cdot (\lambda \vec{\nabla} T) + \alpha_T T \left( \frac{\partial P}{\partial t} + \vec{u} \cdot \vec{\nabla} P \right) \quad (15) $$

The eq. 15 is solved in two steps. First, an intermediary temperature $(T^*)$ is computed containing the convective and diffusive terms of the equations. Then, after the resolution of the Poisson equation for the pressure, the remaining term is added to $T^*$.

Strong temperature gradients can occur at the interface, which are a source of error and can lead to numerical heating. To avoid this, the Aslam's extension method [1] is used to extrapolate
been modified to take into account the compressible effects, i. e., density gradients due to an increase of pressure in an enclosed environment containing an evaporating liquid and spatial variation of the temperature field.

First, an intermediate velocity for both liquid and gas field ($\vec{u}_l^*, \vec{u}_g^*$) is computed by solving the momentum equation (eq.1) without the pressure:

$$ \left\{ \begin{aligned} \vec{u}_l^* &= \vec{u}_l^m - \Delta t \left( \left( \vec{u}_l^m \cdot \vec{\nabla} \right) \vec{u}_l^m - \vec{\nabla} \cdot \left( 2\mu\bar{\epsilon}_l - \frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}_l^m \bar{I} \right) \right) \\ \vec{u}_g^* &= \vec{u}_g^m - \Delta t \left( \left( \vec{u}_g^m \cdot \vec{\nabla} \right) \vec{u}_g^m - \vec{\nabla} \cdot \left( 2\mu\bar{\epsilon}_g - \frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}_g^m \bar{I} \right) \right) \end{aligned} \right. \quad (10) $$

the term $\frac{2}{3}\mu\vec{\nabla} \cdot \vec{u}^m \bar{I}$ had been added to adapt the projection method to a compressible formalism.

Then, the second term of the RHS for Poisson equation is obtained with:

$$ f = \begin{cases} \vec{\nabla} \cdot \vec{u}_l^* & \text{if } \phi > 0, \\ \vec{\nabla} \cdot \vec{u}_g^* & \text{if } \phi < 0, \end{cases} \quad (11) $$

The Poisson equation is computed with following strategy:

$$ -\vec{\nabla} \left( \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} \right) = \vec{\nabla} \cdot \vec{u}^{n+1} - f \quad (12) $$

where $\vec{\nabla} \cdot \vec{u}^{n+1}$ is computed with:

$$ \vec{\nabla} \cdot \vec{u}^{n+1} = - \left( \frac{\alpha_l}{\rho_l C_l^2} + \frac{\alpha_g}{\rho_g C_g^2} \right) \frac{DP}{Dt} - \left( \frac{\alpha_g \alpha_{Tg}}{\rho_g c_{pg}} + \frac{\alpha_l \alpha_{Tl}}{\rho_l c_{pl}} \right) \vec{\nabla} \cdot (\lambda \vec{\nabla} T) \quad (13) $$

The convective term of the pressure material derivative ($\vec{u} \cdot \vec{\nabla} P$) is discretized with a 5th order WENO scheme. The Poisson solver used to solve this Helmholtz type equation 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:

$$ \left\{ \begin{aligned} \vec{u}_l^{n+1} &= \vec{u}_l^* - \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} & \text{if } \phi > 0, \\ \vec{u}_g^{n+1} &= \vec{u}_g^* - \Delta t \frac{\vec{\nabla} P^{n+1}}{\rho^{n+1}} & \text{if } \phi < 0, \end{aligned} \right. \quad (14) $$

The next steps are dedicated to the computation of the liquid and gas velocity extensions ($\vec{u}_l^{\text{ghost}}$, $\vec{u}_g^{\text{ghost}}$). For this, a procedure similar to [12] is used.

## Heat and mass transfer description

To follow the evolution of the temperature field, in the energy equation (eq. 2), the term responsible for the energy generation by friction ($\bar{\tau} \vec{\nabla} \vec{u}$) and the term of volumetric energy generation ($\dot{Q}$) are neglected in this work:

$$ \rho c_p \left( \frac{\partial T}{\partial t} + \vec{u} \cdot \vec{\nabla} T \right) = \vec{\nabla} \cdot (\lambda \vec{\nabla} T) + \alpha_T T \left( \frac{\partial P}{\partial t} + \vec{u} \cdot \vec{\nabla} P \right) \quad (15) $$

The eq. 15 is solved in two steps. First, an intermediary temperature ($T^*$) is computed containing the convective and diffusive terms of the equations. Then, after the resolution of the Poisson equation for the pressure, the remaining term is added to $T^*$.

Strong temperature gradients can occur at the interface, which are a source of error and can lead to numerical heating. To avoid this, the Aslam's extension method [1] is used to extrapolate
the temperature in each phase. Then the extended values are used in the convection term close to the interface.

A 2nd order central difference scheme is used for the discretization of the diffusion term. When mass transfer is considered from the liquid towards the gas (evaporation), the energy jump at the interface (eq. 5) must be taken into account in eq. 15 through a GFM.

For the transport of the vapor mass fraction ($Y_v$), eq. 3 is solved only in the gas phase with $Y_{vs}$ as a Dirichlet boundary condition at the interface. Eq. 3 is reduced to:

$$ \frac{\partial Y_v}{\partial t} + (\vec{u} \cdot \vec{\nabla}) Y_v = \frac{\vec{\nabla} \cdot (\rho_g D_m \vec{\nabla} Y_v)}{\rho_g} \quad (16) $$

As for eq. 15, the diffusion term is discretized using a 2nd order central difference scheme. To account for the Dirichlet boundary condition at the interface, an 1st order Aslam-Chiu extension method is implemented, similar to the work of [2]. The interface temperature $T_Γ$ needed for the Clausius-Clapeyron relation is estimated with a linear interpolation of the liquid temperature extension.

## Results and discussion

### $D^2$ law

In order to validate the thermodynamics of the incompressible formalism, simulation results will be compared with the well-known $D^2$ law. The classical result is that the square of the droplet diameter decreases linearly with time; this behavior is described as the $D^2$ law [10]. To arrive to the normalized square radius equation, the continuity equation in spherical coordinates must be solved assuming incompressibility and that the liquid mass flux from the interface to the surrounding air ($\dot{m}' = 4\pi r^2 \rho \vec{u}$) is constant.

Among the conditions necessary for the $D^2$ law to be applicable, the isolated droplet must be in an infinite continuous medium. Translating this condition into a numerical domain is too numerically expensive. It would imply that the drop surface is sufficiently far from the domain boundary, in order to not disturb the spatial profile of the velocity, vapor mass fraction or temperature. To consider the effect of the boundary, we propose to integrate the species and energy conservation equations using the values of each variable at the limit of the numerical domain, yielding to the following equations:

$$ Y_v(r) = 1 + (Y_{vlim} - 1) \exp \left[ - \frac{\dot{m}'}{4\pi\rho D_m} \left( \frac{1}{r} - \frac{1}{r_{lim}} \right) \right] \quad (17) $$

$$ \ln \left( \frac{T(r) - \frac{h_{lv}}{c_p} - \frac{\varphi}{\dot{m}'}}{T_{lim} - \frac{h_{lv}}{c_p} - \frac{\varphi}{\dot{m}'}} \right) = - \frac{\dot{m}'}{4\pi\rho D_g^t} \left( \frac{1}{r} - \frac{1}{r_{lim}} \right) \quad (18) $$

$\dot{m}'$ is computed by evaluating eq. 17 at the interface ($Y_v(R_D) = Y_{vs}$). And $\varphi$ can be interpreted as the heat flow entering the liquid and it is obtained by evaluating eq. 18 at the interface ($T(R_D) = T_Γ$). This way, the length of the numerical domain is considered in our theoretical calculation, making the results more comparable. A water droplet with initial radius of $R_D = 150\mu m$ is considered. The initial droplet temperature is $T_l = 353K$ and the initial gas temperatures are $T_g = 573K$ and $T_g = 373K$. The domain is an cube with a side length of $8R_D$. Simulations are carried out with three mesh sizes: 32³, 48³ and 64³.

On fig. 1, the temperature and the vapor mass fraction fields are presented at $t = 0.08s$ for the finest grid (64³). As expected, a spherical symmetric property around the drop is observed for both fields. The temporal evolution of the normalized diameter of the drop is showed in fig. 2. In the first moments of the simulations, a strong regression speed is observed. This is because the liquid is not initialized at the equilibrium temperature. Then the curve start to
the temperature in each phase. Then the extended values are used in the convection term close to the interface.

A 2nd order central difference scheme is used for the discretization of the diffusion term. When mass transfer is considered from the liquid towards the gas (evaporation), the energy jump at the interface (eq. 5) must be taken into account in eq. 15 through a GFM.

For the transport of the vapor mass fraction ($Y_v$), eq. 3 is solved only in the gas phase with $Y_{vs}$ as a Dirichlet boundary condition at the interface. Eq. 3 is reduced to:

$$ \frac{\partial Y_v}{\partial t} + (\vec{u} \cdot \vec{\nabla}) Y_v = \frac{\vec{\nabla} \cdot (\rho_g D_m \vec{\nabla} Y_v)}{\rho_g} \quad (16) $$

As for eq. 15, the diffusion term is discretized using a 2nd order central difference scheme. To account for the Dirichlet boundary condition at the interface, an 1st order Aslam-Chiu extension method is implemented, similar to the work of [2]. The interface temperature $T_Γ$ needed for the Clausius-Clapeyron relation is estimated with a linear interpolation of the liquid temperature extension.

## Results and discussion

### $D^2$ law

In order to validate the thermodynamics of the incompressible formalism, simulation results will be compared with the well-known $D^2$ law. The classical result is that the square of the droplet diameter decreases linearly with time; this behavior is described as the $D^2$ law [10]. To arrive to the normalized square radius equation, the continuity equation in spherical coordinates must be solved assuming incompressibility and that the liquid mass flux from the interface to the surrounding air ($\dot{m}' = 4\pi r^2 \rho \vec{u}$) is constant.

Among the conditions necessary for the $D^2$ law to be applicable, the isolated droplet must be in an infinite continuous medium. Translating this condition into a numerical domain is too numerically expensive. It would imply that the drop surface is sufficiently far from the domain boundary, in order to not disturb the spatial profile of the velocity, vapor mass fraction or temperature. To consider the effect of the boundary, we propose to integrate the species and energy conservation equations using the values of each variable at the limit of the numerical domain, yielding to the following equations:

$$ Y_v(r) = 1 + (Y_{vlim} - 1) \exp \left[ - \frac{\dot{m}'}{4\pi\rho D_m} \left( \frac{1}{r} - \frac{1}{r_{lim}} \right) \right] \quad (17) $$

$$ \ln \left( \frac{T(r) - \frac{h_{lv}}{c_p} - \frac{\varphi}{\dot{m}'}}{T_{lim} - \frac{h_{lv}}{c_p} - \frac{\varphi}{\dot{m}'}} \right) = - \frac{\dot{m}'}{4\pi\rho D_g^t} \left( \frac{1}{r} - \frac{1}{r_{lim}} \right) \quad (18) $$

$\dot{m}'$ is computed by evaluating eq. 17 at the interface ($Y_v(R_D) = Y_{vs}$). And $\varphi$ can be interpreted as the heat flow entering the liquid and it is obtained by evaluating eq. 18 at the interface ($T(R_D) = T_Γ$). This way, the length of the numerical domain is considered in our theoretical calculation, making the results more comparable. A water droplet with initial radius of $R_D = 150\mu m$ is considered. The initial droplet temperature is $T_l = 353K$ and the initial gas temperatures are $T_g = 573K$ and $T_g = 373K$. The domain is an cube with a side length of $8R_D$. Simulations are carried out with three mesh sizes: 32³, 48³ and 64³.

On fig. 1, the temperature and the vapor mass fraction fields are presented at $t = 0.08 s$ for the finest grid (64³). As expected, a spherical symmetric property around the drop is observed for both fields. The temporal evolution of the normalized diameter of the drop is showed in fig. 2. In the first moments of the simulations, a strong regression speed is observed. This is because the liquid is not initialized at the equilibrium temperature. Then the curve start to
stabilize, resulting in a linear trend as the temperature at the interface approach an equilibrium temperature; as established in the $D^2$ law. This results show the ability of our method to handle strong temperature gradients at the interface efficiently. Fig. 3 show the temperature and vapor mass fraction profile versus a normalized radius ($r^* = r/R_D$) in the gas phase obtained by solving eq. 18 and 17 and the profiles from the three simulations (32³, 48³, 64³). A good agreement with the theory is found for all the simulations.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Temperature and vapor mass fraction fields ($t = 0.08$ s).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normalized square diameter versus time (s). $T_g = 573$ K: (Solid) 64³, (Dashed) 32³; $T_g = 373$ K: (Dashdot) 64³, (Dotted) 32³.</figcaption></figure>

## HIT configuration with phase change

Similarly to previous works, the idea is to study the influence of the evaporation phenomenon inside a compressible two-phase flows HIT (Homogeneous Isotropic Turbulence). In [5], 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 increase of the pressure due to the evaporation process was not considered. An important improvement has been proposed by [4] with a compressible formalism, where a phase change source term is directly introduced in the continuity and pressure equations; allowing a good conservation of the total mass while introducing the velocity jump at the interface. However, the influence of the temperature in the vapor mass fraction at the interface has not been implemented to the general formalism. Contrary to [4] the velocity jump term is introduced implicitly in the ghost points of the gas velocity field. The coupling between the energy, the continuity and the momentum equations is added by assuming thermodynamical equilibrium at the interface.

For this HIT compressible configuration, the forcing used in previous work [4] is applied. The thermodynamic properties are those of n-decane. The initial gas density is $\rho_g = 25 \text{ kg.m}^{-3}$.
stabilize, resulting in a linear trend as the temperature at the interface approach an equilibrium temperature; as established in the $D^2$ law. This results show the ability of our method to handle strong temperature gradients at the interface efficiently. Fig. 3 show the temperature and vapor mass fraction profile versus a normalized radius ($r^* = r/R_D$) in the gas phase obtained by solving eq. 18 and 17 and the profiles from the three simulations (32³, 48³, 64³). A good agreement with the theory is found for all the simulations.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Temperature and vapor mass fraction fields ($t = 0.08$ s).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normalized square diameter versus time (s). $T_g = 573$ K: (Solid) 64³, (Dashed) 32³; $T_g = 373$ K: (Dashdot) 64³, (Dotted) 32³.</figcaption></figure>

## HIT configuration with phase change

Similarly to previous works, the idea is to study the influence of the evaporation phenomenon inside a compressible two-phase flows HIT (Homogeneous Isotropic Turbulence). In [5], 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 increase of the pressure due to the evaporation process was not considered. An important improvement has been proposed by [4] with a compressible formalism, where a phase change source term is directly introduced in the continuity and pressure equations; allowing a good conservation of the total mass while introducing the velocity jump at the interface. However, the influence of the temperature in the vapor mass fraction at the interface has not been implemented to the general formalism. Contrary to [4] the velocity jump term is introduced implicitly in the ghost points of the gas velocity field. The coupling between the energy, the continuity and the momentum equations is added by assuming thermodynamical equilibrium at the interface.

For this HIT compressible configuration, the forcing used in previous work [4] is applied. The thermodynamic properties are those of n-decane. The initial gas density is $\rho_g = 25 \text{ kg.m}^{-3}$.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Temperature and vapor mass fraction profile in the gaseous phase (t = 0.08 s). Solid line: theoretical solution of D² law, + : 32³, ▷ : 48³ x : 64³.</figcaption></figure>

The pressure is $P = 4.11 \times 10^6 Pa$ and for the Tait equation, the parameters used are $B = 10^9 Pa$, $P_0 = 10^5 Pa$, $\rho_0 = 750 kg.m^{-3}$ and the adiabatic index for the liquid is $\gamma_l = 1.215$. The initial liquid volume fraction is $\Phi = 10\%$. The mesh size is $256^3$ and the targeted mean kinetic energy is $\bar{k} = 3.6 m^2.s^{-2}$. Fig. 4 shows a instantaneous temperature and vapor mass fraction fields and the iso-contour of the level-set function which represent 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 a increase in pressure due to the creation of vapor and a decrease in temperature because of the energy consumption induce by the phase change. Additionally, the influence of the convection due to the of the turbulent velocity field in the temperature field is observed. Similar phenomenons are found in the vapor mass fraction field.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Instantaneous temperature and vapor mass fraction fields and interface location (Level Set 0 isocontour)</figcaption></figure>

## Conclusions

The development and applications of a pressure based method for turbulent weakly-compressible two-phase flows with phase change are presented. The CLSVOF interface capturing method allows to combine the mass conservation properties of the VoF method and the sharp interface capturing of the LS method. For considering the spatial density gradients and the evaporation process the standard VoF equation has been modified. Also, the estimation of a phase change sink term using the local interface surface combine with a redistribution algorithm allow us to take into account the liquid evaporation; while assuring mass conservation.
A challenging task when simulating evaporating two-phase flows is the treatment of the addi-
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Temperature and vapor mass fraction profile in the gaseous phase (t = 0.08 s). Solid line: theoretical solution of D² law, + : 32³, ▷ : 48³ x : 64³.</figcaption></figure>

The pressure is $P = 4.11 \times 10^6 Pa$ and for the Tait equation, the parameters used are $B = 10^9 Pa$, $P_0 = 10^5 Pa$, $\rho_0 = 750 kg.m^{-3}$ and the adiabatic index for the liquid is $\gamma_l = 1.215$. The initial liquid volume fraction is $\Phi = 10\%$. The mesh size is $256^3$ and the targeted mean kinetic energy is $\bar{k} = 3.6 m^2.s^{-2}$. Fig. 4 shows a instantaneous temperature and vapor mass fraction fields and the iso-contour of the level-set function which represent 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 a increase in pressure due to the creation of vapor and a decrease in temperature because of the energy consumption induce by the phase change. Additionally, the influence of the convection due to the of the turbulent velocity field in the temperature field is observed. Similar phenomenons are found in the vapor mass fraction field.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Instantaneous temperature and vapor mass fraction fields and interface location (Level Set 0 isocontour)</figcaption></figure>

## Conclusions

The development and applications of a pressure based method for turbulent weakly-compressible two-phase flows with phase change are presented. The CLSVOF interface capturing method allows to combine the mass conservation properties of the VoF method and the sharp interface capturing of the LS method. For considering the spatial density gradients and the evaporation process the standard VoF equation has been modified. Also, the estimation of a phase change sink term using the local interface surface combine with a redistribution algorithm allow us to take into account the liquid evaporation; while assuring mass conservation.
A challenging task when simulating evaporating two-phase flows is the treatment of the addi-
tional jump conditions found at the interface. This is taken care by the use of accurate extensions of the discontinuous variables and the ghost fluid method. Implementation of the previous mentioned methods helps to minimize numerical diffusion across the interface and avoid artificial heating.

The thermodynamics and species evolution of the method is validated by comparing simulation results with the $D^2$ law. Finally, the implementation of the evaporation process in a 3D HIT compressible configuration present the ability of our method to handle strong interface curvature, thermal dilatation, collisions and breakups. Future work will be dedicated to the study of the interaction of the Stefan flow with the turbulence characteristics.

## Acknowledgements

This work was granted access to the HPC resources of IDRIS, TGCC and CINES under the allocation A0092B10101 made by GENCI (Grand Equipement National de Calcul Intensif). The CRIANN (Centre Regional Informatique et d'Applications Numeriques de Normandie) is also gratefully acknowledged (project No. 2017004) for their CPU resources.

## Nomenclature

<table><tr><td>γ<sub>l</sub></td><td>Adiabatic index</td></tr><tr><td>c<sub>p</sub></td><td>Heat capacity [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>λ</td><td>Thermal conductivity [w.m<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>α<sub>T</sub></td><td>Thermal dilatation coefficient [K<sup>-1</sup>]</td></tr><tr><td>τ̄</td><td>Viscous strain tensor [kg.m<sup>-1</sup>.s<sup>-2</sup>]</td></tr><tr><td>Y<sub>v</sub></td><td>Vapor mass fraction</td></tr><tr><td>D<sub>m</sub></td><td>Mass diffusivity [m<sup>2</sup>.s<sup>-1</sup>]</td></tr><tr><td>Y<sub>vs</sub></td><td>Vapor mass fraction at the interface</td></tr><tr><td>ṁ</td><td>Volumetric evaporation rate [kg.m<sup>-3</sup>.s<sup>-1</sup>]</td></tr><tr><td>α<sub>l</sub></td><td>Liquid volume fraction</td></tr><tr><td>Σ</td><td>Surface density [m<sup>-1</sup>]</td></tr><tr><td>ω̇</td><td>Evaporation rate [Kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>h<sub>lv</sub></td><td>Latent heat [J.Kg<sup>-1</sup>]</td></tr><tr><td>φ</td><td>Level-set function [m]</td></tr></table>

## References

[1] ASLAM, T. D. A partial differential equation approach to multidimensional extrapolation. *Journal of Computational Physics* 193, 1 (2004), 349–355.

[2] BOUALI, Z., DURET, B., DEMOULIN, F.-X., AND MURA, A. Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows. *International Journal of Multiphase Flow* 85 (2016), 326–335.

[3] CALIMEZ, X. Simulation d'une échelle par une méthode VOF d'écoulements diphasiques réactifs. PhD thesis, Châtenay-Malabry, Ecole centrale de Paris, 1998.

[4] DURET, B., CANU, R., REVEILLON, J., AND DEMOULIN, F. A pressure based method for vaporizing compressible two-phase flows with interface capturing approach. *International Journal of Multiphase Flow* 108 (2018), 42–50.

[5] DURET, B., LURET, G., REVEILLON, J., MÉNARD, T., BERLEMONT, A., AND DEMOULIN, F.-X. Dns analysis of turbulent mixing in two-phase flows. *International Journal of Multiphase Flow* 40 (2012), 93–105.

[6] MALAN, L. Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change. PhD thesis, University of Cape Town, 2018.

[7] MÉNARD, T., TANGUY, S., AND BERLEMONT, A. Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet. *International Journal of Multiphase Flow* 33, 5 (May 2007), 510–524.

[8] PALMORE JR, J., AND DESJARDINS, O. A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows. *Journal of Computational Physics* 399 (2019), 108954.

[9] SCAPIN, N., COSTA, P., AND BRANDT, L. A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. *Journal of Computational Physics* 407 (2020), 109251.

[10] SIRIGNANO, W. A., AND EDWARDS, C. F. Fluid dynamics and transport of droplets and sprays. *J. Fluids Eng.* 122, 1 (2000), 189–190.

[11] SUSSMAN, M., AND PUCKETT, E. G. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. *Journal of computational physics* 162, 2 (2000), 301–337.

[12] TANGUY, S., MÉNARD, T., AND BERLEMONT, A. A level set method for vaporizing two-phase flows. *Journal of Computational Physics* 221, 2 (2007), 837–853.
tional jump conditions found at the interface. This is taken care by the use of accurate extensions of the discontinuous variables and the ghost fluid method. Implementation of the previous mentioned methods helps to minimize numerical diffusion across the interface and avoid artificial heating.

The thermodynamics and species evolution of the method is validated by comparing simulation results with the $D^2$ law. Finally, the implementation of the evaporation process in a 3D HIT compressible configuration present the ability of our method to handle strong interface curvature, thermal dilatation, collisions and breakups. Future work will be dedicated to the study of the interaction of the Stefan flow with the turbulence characteristics.

## Acknowledgements

This work was granted access to the HPC resources of IDRIS, TGCC and CINES under the allocation A0092B10101 made by GENCI (Grand Equipement National de Calcul Intensif). The CRIANN (Centre Regional Informatique et d'Applications Numeriques de Normandie) is also gratefully acknowledged (project No. 2017004) for their CPU resources.

## Nomenclature

<table><tr><td>γ<sub>l</sub></td><td>Adiabatic index</td></tr><tr><td>c<sub>p</sub></td><td>Heat capacity [J.kg<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>λ</td><td>Thermal conductivity [w.m<sup>-1</sup>.K<sup>-1</sup>]</td></tr><tr><td>α<sub>T</sub></td><td>Thermal dilatation coefficient [K<sup>-1</sup>]</td></tr><tr><td>τ̄</td><td>Viscous strain tensor [kg.m<sup>-1</sup>.s<sup>-2</sup>]</td></tr><tr><td>Y<sub>v</sub></td><td>Vapor mass fraction</td></tr><tr><td>D<sub>m</sub></td><td>Mass diffusivity [m<sup>2</sup>.s<sup>-1</sup>]</td></tr><tr><td>Y<sub>vs</sub></td><td>Vapor mass fraction at the interface</td></tr><tr><td>ṁ</td><td>Volumetric evaporation rate [kg.m<sup>-3</sup>.s<sup>-1</sup>]</td></tr><tr><td>α<sub>l</sub></td><td>Liquid volume fraction</td></tr><tr><td>Σ</td><td>Surface density [m<sup>-1</sup>]</td></tr><tr><td>ω̇</td><td>Evaporation rate [Kg.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>h<sub>lv</sub></td><td>Latent heat [J.Kg<sup>-1</sup>]</td></tr><tr><td>φ</td><td>Level-set function [m]</td></tr></table>

## References

[1] ASLAM, T. D. A partial differential equation approach to multidimensional extrapolation. *Journal of Computational Physics* 193, 1 (2004), 349–355.

[2] BOUALI, Z., DURET, B., DEMOULIN, F.-X., AND MURA, A. Dns analysis of small-scale turbulence-scalar interactions in evaporating two-phase flows. *International Journal of Multiphase Flow* 85 (2016), 326–335.

[3] CALIMEZ, X. Simulation a petite échelle par une méthode VOF d'écoulements diphasiques réactifs. PhD thesis, Châtenay-Malabry, Ecole centrale de Paris, 1998.

[4] DURET, B., CANU, R., REVEILLON, J., AND DEMOULIN, F. A pressure based method for vaporizing compressible two-phase flows with interface capturing approach. *International Journal of Multiphase Flow* 108 (2018), 42–50.

[5] DURET, B., LURET, G., REVEILLON, J., MÉNARD, T., BERLEMONT, A., AND DEMOULIN, F.-X. Dns analysis of turbulent mixing in two-phase flows. *International Journal of Multiphase Flow* 40 (2012), 93–105.

[6] MALAN, L. Direct numerical simulation of free-surface and interfacial flow using the VOF method: cavitating bubble clouds and phase change. PhD thesis, University of Cape Town, 2018.

[7] MÉNARD, T., TANGUY, S., AND BERLEMONT, A. Coupling level set/VOF/ghost fluid methods: Validation and application to 3d simulation of the primary break-up of a liquid jet. *International Journal of Multiphase Flow* 33, 5 (May 2007), 510–524.

[8] PALMORE JR, J., AND DESJARDINS, O. A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows. *Journal of Computational Physics* 399 (2019), 108954.

[9] SCAPIN, N., COSTA, P., AND BRANDT, L. A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. *Journal of Computational Physics* 407 (2020), 109251.

[10] SIRIGNANO, W. A., AND EDWARDS, C. F. Fluid dynamics and transport of droplets and sprays. *J. Fluids Eng.* 122, 1 (2000), 189–190.

[11] SUSSMAN, M., AND PUCKETT, E. G. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. *Journal of computational physics* 162, 2 (2000), 301–337.

[12] TANGUY, S., MÉNARD, T., AND BERLEMONT, A. A level set method for vaporizing two-phase flows. *Journal of Computational Physics* 221, 2 (2007), 837–853.