Multicomponent Droplet Evaporation in a Geometric
Volume-Of-Fluid Framework

Edoardo Cipriano<sup>a,*</sup>, Abd Essamade Saufi<sup>a,1</sup>, Alessio Frassoldati<sup>a</sup>, Tiziano Faravelli<sup>a</sup>, Stéphane Popinet<sup>b</sup>, Alberto Cuoci<sup>a</sup>

<sup>a</sup>CRECK Modeling Lab, Department of Chemistry, Materials, and Chemical Engineering
"G. Natta", Politecnico di Milano, Piazza Leonardo da Vinci, 32, Milano, 20133, Italy

<sup>b</sup>Institut Jean Le Rond d'Alembert, CNRS UMR 7190, Sorbonne
Université, Paris, 75005, France

Abstract

This work proposes an innovative model for multicomponent phase change
in interface-resolved simulations. The two-phase system is described by a
geometric Volume-Of-Fluid (VOF) approach, and considers multiple compo-
nents in non-isothermal environments, relaxing the hypothesis of pure liquid
droplets usually studied in the literature. The model includes the Stefan flow
and implements the following solutions for the complications that arise when
studying liquid mixtures: i) a coupled approach for solving the interface jump
conditions; ii) a proper strategy to obtain a liquid velocity for the advection
of the volume fraction field, also applicable to static droplets with strong
density ratio; iii) and a geometric approach to discretize the scalar fields
transport equations. This model was implemented in the open-source code
Basilisk, and it was tested on a number of benchmark phase change prob-
lems, such as the fixed flux evaporation, the Stefan problem, Epstein Plesset,
and the Scriven problem. These test cases demonstrate the convergence of
the numerical methods to the analytical solutions. More complex configu-
rations, such as multicomponent isothermal and non-isothermal droplets are
compared using numerical benchmark solutions obtained assuming spherical
symmetry. The code, as well as the simulation setups, are released on the

<footnote>*Corresponding author.</footnote>

<footnote>Email address: edoardo.cipriano@polimi.it (Edoardo Cipriano)</footnote>

<footnote n="1">¹Present address: First Light Fusion Ltd., Unit 9/10, Oxford Pioneer Park, Mead Road, Oxford, OX5 1QU, Oxfordshire, United Kingdom (UK)</footnote>
Basilisk website, making it the first model and open-source implementation of multicomponent phase change in a VOF framework.

Keywords: VOF, droplet, phase change, evaporation, multicomponent

# 1. Introduction

Multicomponent liquid droplets are widespread in nature and have been studied at many levels of detail in a number of engineering applications. Most of the literature models rely on the approximation of spherical symmetry [1, 2], neglecting the effects of surface tension and droplet deformation. This simplification is accurate for droplets in microgravity conditions, where the effects of droplet deformation and buoyancy on the heat and mass transfer phenomena are negligible [3]. Therefore, they are useful for the analysis of complex chemical phenomena, such as combustion processes [4], but due to their limitations, multidimensional models describing phase change in gas-liquid systems are required.

The Geometric Volume-Of-Fluid (VOF) approach [5] is one of the most popular models available in the literature to describe multiphase systems. It allows the interface to be transported conservatively by solving a continuity equation and using detailed discretization methods that preserve the interface sharpness. The VOF model has been used for a wide variety of simulations involving gas-liquid systems, for example in free surface flows, breaking waves, atomizations, liquid jets, and droplet coalescence and breakup [6, 7, 8].

More recently, the VOF approach has also been considered for the simulation of two-phase systems with phase change. Phase change poses several numerical challenges that mainly arise from the presence of an expansion term localized at the gas-liquid interface. Such expansion leads to an additional convective flux, usually referred to as Stefan flow, which introduces a strong and localized discontinuity in the velocity field. Therefore, the transport of the interface and of the scalar fields must be carefully resolved in order to avoid interface smearing due to the expansion term.

The first attempt to include phase change in a Volume-Of-Fluid model was made by Welch and Wilson [9], who adapted the Front Tracking model for simulating boiling flows proposed by Juric and Tryggvason [10] and Son and Dhir [11] to a VOF framework to study the film boiling configuration [12, 13]. A similar approach was then adopted by Hardt and Wondra [14] and by Kunkelmann [15], who used a specific technique for the smearing of the source
terms given by the phase change, removing the volume expansion effect from
the interfacial cells. Consequently, the phase change expansion term does not
interfere with the interface advection process [15]. Their boiling model was
used for benchmark boiling evaporation problems, like the Stefan problem
and the film boiling [14], and for the simulation of nucleated bubbles [15].
Schlottke and Weigand [16] proposed an iterative approach to impose the
velocity jump conditions in the governing equations. This approach was
used for the simulation of vaporizing droplets with large deformations at high
Reynolds and Weber numbers. More recent works by Scapin et al. [17] and
Malan et al. [18] proposed similar methods to obtain a conservative advection
velocity from the one-field velocity by solving an additional Poisson equation.
These models focus on thermally-induced [18] or chemical species-induced
phase change [17], which simplifies the calculation of the vaporization rate,
and they were tested on benchmark evaporation problems like Stefan and
Scriven problems, and the *d-squared* law. Other literature models, such
as those proposed by Palmore and Desjardins [19], and Zhao et al. [20],
solve the coupled heat and mass transfer focusing on pure liquid droplets
and neglecting the presence of a generic number of multiple chemical species
both in the gas and in the liquid phases. The extension of these phase-change
models to multiple-species systems is not straight-forward, and it comes with
several challenges:

1. Multiple chemical species add non-linearity to the interface jump condition. Interface temperature and mass fractions, in gas and in liquid phase, are not known a-priori, and they must be obtained coupling energy and mass balances in a non-linear system of equations.

2. The scalar fields must be transported in both the gas and the liquid phase, considering the presence of a discontinuity in the velocity field and enforcing the boundary conditions at the gas-liquid interface.

3. A one-field approach for the species transport cannot be used, since the evaporation of a light species must cause the accumulation of a heavy species, and conservation of mass in the two different phases must be respected (Figure 1).

The numerical model proposed in this work addresses these challenges using a multicomponent phase change formulation, which is valid for an arbitrary number of chemical species in both the gas and liquid phases. The method proposed to obtain an extended liquid velocity for the transport of the volume fraction and the scalar fields, is also applicable to static droplets,
<figure><img src="image_1.png" /><figcaption>Figure 1: Species mass fraction profiles in a spherical liquid droplet. The green line is the mass fraction profile of a light component in liquid phase, the red line is the mass fraction of the same light component in gas phase, while the blue line is the mass fraction profile of the heavy chemical species in liquid phase.</figcaption></figure>

which is the most challenging phase change configuration. For static droplets, the velocity field is dominated by the Stefan flow, which must be handled with care to avoid non-physical interface transport. The mathematical model is described in Section 2, while Section 3 describes the numerical discretization of the transport equations. The resulting model is implemented in the Basilisk framework [21] and it is validated using benchmark boiling and droplet evaporation problems (Section 4).

## 2. Mathematical Formulation

Figure 2 depicts the control volume over which the transport equations are written. This volume includes two immiscible phases ($V_l$ and $V_g$) separated by a zero-thickness interface $\Gamma$. The different phases inside the same control volume $V$ are identified by introducing a scalar marker function defined in the whole domain ($V = V_l \cup V_g$) as an Heaviside function:

$$H(\mathbf{x}, t) = \begin{cases} 1, & \text{if } \mathbf{x} \in V_l \\ 0, & \text{if } \mathbf{x} \in V_g \end{cases} \qquad (1)$$

Exploiting the definition of the marker function, a generic physical property $\phi$ can be defined over the entire control volume:

$$\phi(\mathbf{x}, t) = \phi_l H(\mathbf{x}, t) + \phi_g (1 - H(\mathbf{x}, t)) \qquad (2)$$
<figure><img src="image_2.png" /><figcaption>Figure 2: Generic domain representation. Adapted from [22].</figcaption></figure>

where $\phi_l$ and $\phi_g$ are the phase properties (i.e. density, viscosity, thermal conductivity, specific heat, and diffusivity of the chemical species), for the liquid and for the gas phase, respectively. This work assumes the properties of the gas and the liquid phase to be constant during the entire simulation. Therefore, the generic property $\phi$ only changes according to the marker function $H(\mathbf{x}, t)$.

## 2.1. Interface Advection with Phase Change

Computationally, the marker function is approximated using the color function $c$ [23], defined as the volume fraction of the reference (liquid) phase in the control volume $V$:

$$c = \frac{1}{V} \int_V H(\mathbf{x}, t) dV \qquad (3)$$

which assumes value 1 in the liquid phase and 0 in the gas phase. According to the VOF approach, the advection of the interface can be performed by solving a conservation equation for the volume fraction. This equation can be derived from a mass balance on the liquid control volume $V_l$:

$$\frac{\partial}{\partial t} \int_{V_l} \rho_l dV_l + \oint_{S_l} \rho_l \mathbf{u}_l \cdot \mathbf{n} dS_l = - \oint_{S_\Gamma} \rho_l (\mathbf{u}_l - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma dS_\Gamma \quad (4)$$

where the term on the RHS is the liquid mass flux across the interface due to the phase change, and it can be expressed as a function of the vaporization rate, using the Rankine-Hugoniot relation [23]:

$$\rho_l (\mathbf{u}_l - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma = \rho_g (\mathbf{u}_g - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma = \dot{m} \qquad (5)$$
which is a mass balance across a moving gas-liquid interface. The last term on the RHS, $\dot{m}$, is the total vaporization rate per unit of surface. Introducing Equation 5 in 4, and assuming that the liquid density depends on c only, we obtain the transport equation for the volume fraction of the liquid phase:

$$ \frac{\partial}{\partial t} \int_V c dV + \oint_S c \mathbf{u}_l \cdot \mathbf{n} dS = - \int_V \frac{\dot{m}}{\rho_l} \delta_{\Gamma} dV \quad (6) $$

The term on the RHS is the source (or sink) of liquid volume according to the phase change phenomena. This equation can be solved given the velocity field of the liquid phase $\mathbf{u}_l$ and the total evaporation rate per unit of surface $\dot{m}$. The term $\delta_{\Gamma}$ applies the source only at the interface, and it is defined as the interface surface in the control volume divided by the volume itself, with resulting units of area over volume.

## 2.2. Incompressible Navier-Stokes Equations with Phase Change

The velocity and pressure fields are obtained solving the Navier-Stokes equations for incompressible flows with phase change. The incompressibility of the flow field is justified for both phases in the control volume due to the low Mach number ($Ma \ll 0.3$) of the systems of interest for this work. The momentum equation is therefore written as

$$ \frac{\partial}{\partial t} \int_V \rho \mathbf{u} dV + \oint_S \rho \mathbf{u} (\mathbf{u} \cdot \mathbf{n}) dS = \oint_S \boldsymbol{\tau} \cdot \mathbf{n} dS + \int_V \rho g dV + \int_V \mathbf{f}_{\sigma} dV \quad (7) $$

following the so-called one-field formulation, which means that a single equation is integrated over the entire control volume V using physical properties that vary abruptly at the gas-liquid interface (refer to [23] for a detailed derivation). The stress tensor is defined as $\boldsymbol{\tau} = \mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T) - p\mathbf{I}$, where $\mathbf{I}$ is the unit tensor, while the material properties $\rho$ and $\mu$ are the density and dynamic viscosity respectively, defined as in Equation 2.

The one-field continuity equation with phase change can be derived by writing the continuity equation (Equation 4) separately for the gas phase and for the liquid phase, and summing up the two expressions [9]. Considering incompressible flows and exploiting the interface mass balance (Equation 5), the resulting one-field continuity equation is obtained:

$$ \oint_S \mathbf{u} \cdot \mathbf{n} dS = \int_V \dot{m} \left( \frac{1}{\rho_g} - \frac{1}{\rho_l} \right) \delta_{\Gamma} dV \quad (8) $$
where the RHS of the equation is a source term localized at the gas-liquid interface which, physically, derives from the fact that the conversion from liquid to vapor during phase change (or vice versa) is accompanied by a significant change of volume. This process leads to an expansion term that makes the one-field velocity $\mathbf{u}$ non divergence-free at the interface.

### 2.3. Species Equations

Chemical species are transported by solving the conservation equation in terms of mass fractions written using a two-field formulation. Using this approach, two different governing equations are solved for the same species, one for the liquid phase and one for the gas phase. This approach was chosen because, when solving multiple chemical species, the mass fraction of a single species varies inversely in the two phases. This effect cannot be described using a one-field formulation with source term, which can only increase or decrease the mass fraction in both phases. Moreover, the separation of the fields limits the numerical diffusion across the interface. The species equation for a single species $i$ in the phase $k$ is derived by integrating the mass of the chemical species over the control volume $V_k$, with $k = l, g$:

$$ \frac{\partial}{\partial t} \int_{V_k} \rho_k \omega_{i,k} dV_k + \oint_{S_k} \rho_k \omega_{i,k} \mathbf{u}_k \cdot \mathbf{n} dS_k + \oint_{S_\Gamma} \rho_k \omega_{i,k} (\mathbf{u}_k - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma dS_\Gamma = \\ \oint_{S_k} \rho_k \mathcal{D}_{i,k} \nabla \omega_{i,k} \cdot \mathbf{n} dS_k + \oint_{S_\Gamma} \rho_k \mathcal{D}_{i,k} \left. \frac{\partial \omega_{i,k}}{\partial \mathbf{n}_\Gamma} \right|_k \cdot \mathbf{n}_\Gamma dS_\Gamma \quad (9) $$

where $\omega_{i,k}$ is the mass fraction of the $i$-th species in the $k$-th phase, while $\mathcal{D}_{i,k}$ is the diffusion coefficient. Analogously to the volume fraction transport equation, the integral over the interface surface is replaced by the interface mass balance, written for a single chemical species:

$$ \dot{m}_i = \rho_k \omega_{i,k} (\mathbf{u}_k - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma - \rho_k \mathcal{D}_{i,k} \left. \frac{\partial \omega_{i,k}}{\partial \mathbf{n}_\Gamma} \right|_k \cdot \mathbf{n}_\Gamma \quad (10) $$

Introducing Equation 10 in 9, the phase change contribution is collected into the term $\dot{m}_i$, which is the vaporization rate per unit of surface, for the single chemical species $i$:

$$ \frac{\partial}{\partial t} \int_V c_k \rho_k \omega_{i,k} + \oint_S c_k \rho_k \omega_{i,k} \mathbf{u}_k \cdot \mathbf{n} dS = \oint_S c_k \rho_k \mathcal{D}_{i,k} \nabla \omega_{i,k} \cdot \mathbf{n} dS - \int_V \dot{m}_i \delta_\Gamma dV \quad (11) $$
where the volume fraction $c_k$ is used to integrate Equation 9 over the whole control volume. This equation is solved for the gas and for the liquid phase separately and, eventually, the one-field mass fraction is reconstructed for visualization purposes using an arithmetic average based on the color function:

$$ \omega_i = \omega_{i,l}c + \omega_{i,g}(1-c) \quad (12) $$

## 2.4. Energy Equation

The energy transport equation is solved for temperature using a two-field formulation, following the same approach used for the transport of the chemical species mass fractions. The temperature equation can be derived from an energy balance on the control volume $V_k$, being $k = l, g$, with the assumption of Newtonian fluids and negligible viscous effects and pressure work term [24]:

$$ \frac{\partial}{\partial t} \int_{V_k} \rho_k c_{p,k} T_k dV_k + \oint_{S_k} \rho_k c_{p,k} T_k \mathbf{u}_k \cdot \mathbf{n} dS_k + \oint_{S_\Gamma} \rho_k c_{p,k} T_k (\mathbf{u}_k - \mathbf{u}_\Gamma) \cdot \mathbf{n}_\Gamma dS_\Gamma = \\ \oint_{S_k} k_k \nabla T_k \cdot \mathbf{n} dS_k + \oint_{S_\Gamma} k_k \frac{\partial T_k}{\partial \mathbf{n}_\Gamma} \Big|_k \cdot \mathbf{n}_\Gamma dS_\Gamma \quad (13) $$

where $T_k$ is the temperature of the k-th phase, $c_{p,k}$ is the constant pressure heat capacity, and $k_k$ is the thermal conductivity. Following the same approach explained in [20], introducing Equation 5 into 13, and integrating over the whole control volume, the two-field energy balance can be re-written as:

$$ \frac{\partial}{\partial t} \int_V c_k \rho_k c_{p,k} T_k dV + \oint_S c_k \rho_k c_{p,k} T_k \mathbf{u}_k \cdot \mathbf{n} dS = \\ \oint_S c_k k_k \nabla T_k \cdot \mathbf{n} dS + \int_V (\dot{q}_{\Gamma,k} - \dot{m} c_{p,k} T_k) \delta_\Gamma dV \quad (14) $$

where $q_{\Gamma,k}$ is the heat conduction across the interface, which is directly computed from the interface gradients on both sides, once the interface temperature value is obtained from the jump condition.

## 2.5. Multicomponent Interface Jump Condition

The vaporization rates appearing in Equation 6, 8, and 11, as well as the interface temperature in Equation 14 must be computed from the interface jump condition. For this purpose, the mass and energy balances are applied
<figure><img src="image_3.png" /><figcaption>Figure 3: Mass and energy fluxes at the gas-liquid interface. $\dot{q}_L$ and $\dot{q}_G$ are the conduction heat fluxes, while $\dot{j}_L$ and $\dot{j}_G$ are the species diffusion fluxes.</figcaption></figure>

over a control volume corresponding to a portion of the gas-liquid interface $\Gamma$ shown in Figure 3. The mass balance at the interface is written for each chemical species and it imposes the conservation of the convective and diffusive fluxes. Instead, the energy balance at the interface equates the difference in conductive fluxes to the energy required for the evaporation process [3]. Therefore, considering that accumulation of mass and energy is not allowed in a zero-thickness interface, and assuming that thermodynamic equilibrium conditions occur at the interface, the following system of equations is obtained:

$$ \left\{ \begin{aligned} \dot{m}_i &= \dot{m}\hat{\omega}_{i,l} - \rho_l \mathcal{D}_{i,l} \left. \frac{\partial \omega_{i,l}}{\partial \mathbf{n}_\Gamma} \right|_l = \dot{m}\hat{\omega}_{i,g} - \rho_g \mathcal{D}_{i,g} \left. \frac{\partial \omega_{i,g}}{\partial \mathbf{n}_\Gamma} \right|_g \\ \sum_{i=1}^{NS} \dot{m}_i \Delta h_{ev,i} &= k_l \left. \frac{\partial T_l}{\partial \mathbf{n}_\Gamma} \right|_l + k_g \left. \frac{\partial T_g}{\partial \mathbf{n}_\Gamma} \right|_g \\ \hat{x}_{i,g} &= k_{eq,i}(\hat{T})\hat{x}_{i,l} \end{aligned} \right. \quad (15) $$

where the hat over the symbols is used to refer to interface quantities: $\hat{x}_{i,g}$ and $\hat{x}_{i,l}$ are the mole fractions of the chemical species $i$, $\hat{\omega}_{i,g}$ and $\hat{\omega}_{i,l}$ are the mass fractions, and $\hat{T}$ is the interface temperature. The phase-equilibrium is computed using Raoult's Law, which defines the thermodynamic equilibrium constant as the ratio between the interface mole fractions in the gas and in the liquid phase, respectively:

$$ k_{eq,i}(\hat{T}) = \frac{\hat{x}_{i,g}}{\hat{x}_{i,l}} = \frac{P_{vap,i}(\hat{T})}{P} \quad (16) $$

where $P$ is the thermodynamic pressure of the system and $P_{vap,i}(\hat{T})$ is the vapor pressure of the generic species $i$, which can be computed using different
vapor-liquid equilibrium relations. Here, Antoine's Law was used:

$$P_{vap,i}(\hat{T}) = \exp \left[ A - \frac{B}{T+C} \right] \qquad (17)$$

where A, B, and C are empirical parameters that depend on the chemical species under investigation. In this work, the thermodynamic equilibrium constant is obtained under the assumptions of an ideal liquid mixture, an ideal gas, and a negligible Poynting correction [25]. Note that this formulation can easily be extended for non-ideal mixtures and more complex Equation of States.

The non-linear system of equations resulting from the interface jump condition is solved to determine the evaporation rate for every chemical species $\dot{m}_i$, the interface mass fractions ($\hat{\omega}_{i,l}$ and $\hat{\omega}_{i,g}$), and the interface temperature $\hat{T}$. The solution of this system of equations provides all the information required for the characterization of the interface and for the calculation of the source terms in the governing equations.

### 3. Numerical Discretization

The governing equations described in Section 2 are discretized on an adaptive Cartesian grid. The domain is divided in a number of square cells (cubic in 3D), and the governing equations are numerically approximated using the Finite Volume Method and a collocated arrangement. Two-dimensional configurations are discussed here, where $\Delta$ is the length of a generic cell, the subscript f is the generic edge (or face) of the cell, while the superscript n indicates the simulation time step. The time discretization used here follows the same formalism adopted in [26, 27, 28], which assumes that the volume fraction and the relative tracers (mass fractions and temperature) are defined at the time level $n - 1/2$, i.e. are lagging the velocity and pressure fields by half time step. The extension of the proposed numerical schemes to three dimensions is straight-forward.

#### 3.1. Geometric Volume-Of-Fluid Advection

Using the VOF approach, the color function described in Section 2 is computed in every cell of the domain, obtaining pure liquid cells ($c = 1$), pure gas cells ($c = 0$), and interfacial cells with intermediate values of volume fraction. The advection of the interface is solved conservatively using the
advection equation for the volume fraction of the liquid phase (Equation 6). In this work, the directional-split Geometric VOF scheme was adopted [29, 30]. This method was developed for Cartesian grids and it is conservative if a divergence-free velocity field is provided. This approach consists of two steps: the reconstruction step, where the interface is geometrically reconstructed in each interfacial cell according to the Piecewise Linear Interface Construction (PLIC) [29], and the flux calculation step, where the liquid volume fluxes across the cell faces are geometrically computed along each direction, and the volume fraction in each cell is updated at the next time level:

$$ \frac{c' - c^{n-1/2}}{\Delta t} = \frac{1}{\Delta} (F(\mathbf{u}_{\Gamma}^n)_{f+} - F(\mathbf{u}_{\Gamma}^n)_{f-}) + \int_V c_c \frac{\partial \mathbf{u}^n}{\partial x} dV \quad (18) $$

$$ \frac{c^{n+1/2} - c'}{\Delta t} = \frac{1}{\Delta} (G(\mathbf{u}_{\Gamma}^n)_{f+} - G(\mathbf{u}_{\Gamma}^n)_{f-}) + \int_V c_c \frac{\partial \mathbf{u}^n}{\partial y} dV \quad (19) $$

The terms F and G in Equation 18, 19 are the geometric fluxes across the faces, while the integrals on the RHS are dilation terms that correct errors in the divergence stemming from the dimensional splitting. Details about the specific implementation can be found in [22, 23, 29].

To avoid interface smearing, the velocity used for the interface advection $\mathbf{u}_{\Gamma}$ must be continuous across the interface. For this reason, we cannot approximate $\mathbf{u}_{\Gamma}$ with the one-field velocity $\mathbf{u}$. Instead, the interface velocity is obtained from Equation 5:

$$ \mathbf{u}_{\Gamma} = \mathbf{u}_l - \frac{\dot{m}}{\rho_l} \mathbf{n}_{\Gamma} = \mathbf{u}_{E,f} - \mathbf{u}_{pc} \quad (20) $$

where $\mathbf{u}_{E,f}$ is a divergence-free liquid velocity, computed as explained in Section 3.3, while $\mathbf{u}_{pc}$ is the velocity contribution due to the phase change, which shrinks or expands the liquid phase. Introducing this definition of interface velocity allows the phase change source term of Equation 6 to be directly included in the advection step. This operation is not straight-forward, because the interface gradients are computed at cell centers, and so is the vaporization rate $\dot{m}$. Therefore, interpolations are needed to transform the vaporization rate in an interface regression velocity $\mathbf{u}_{pc}$, defined on the cell faces. The proposed interpolation scheme is depicted in Figure 4:

1. If the face connects an interfacial cell with a pure cell, the interface regression velocity on that face is taken from the interfacial cell, weighted by the corresponding interface normal component.
<figure><img src="image_4.png" /><figcaption>Figure 4: Interpolation of the interface regression velocity, from cell-centers to cell-faces (http://basilisk.fr/sandbox/ecipriano/test/interfaceregression.c).</figcaption></figure>

2. If the face connects two interfacial cells, the interface regression velocity is computed from a linear interpolation between the two consecutive vaporization rates, weighted by the interface normal.

Using this approach, the resulting interface advection velocity $\mathbf{u}_\Gamma$ transports the volume fraction in a non-conservative way, but the amount of liquid volume added or removed matches the effect of the phase change phenomena. The advantage of this approach with respect to other literature methods [18, 31] is that advection and phase change are applied in a single step, limiting the possibility of under- and over-shoots that easily arise when the phase change is applied as an explicit source term.

### 3.2. Pressure-Velocity Coupling

The Navier-Stokes equations are discretized using an approximate projection method, where velocity and pressure are collocated at cell-centers. This method and its numerical implementation in a quad/octree grid was extensively described by Popinet in [26] and [27]. During the *Predictor Step*, the solution of the advection-diffusion part of the momentum equation is performed by splitting the advection from the diffusion step. The advection part is estimated using the Bell-Colella-Glaz approach [28], which is a second order unsplit upwind scheme that is stable for CFL condition smaller than one [27]. The diffusion part is instead re-arranged into a Poisson-Helmholtz equation and solved in an implicit manner using a multigrid solver [32]. This method is combined with a Crank-Nicholson discretization which is second
order accurate and stable also for CFL number above one [27]. The time staggered scheme reduces to the following discretization:

$$ \rho^{n+1/2} \left( \frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} + \mathbf{u}^{n+1/2} \cdot \nabla \mathbf{u}^{n+1/2} \right) = \nabla \cdot (2\mu^{n+1/2} \mathbf{D}^*) \quad (21) $$

$$ \mathbf{u}_f^* = \mathbf{u}_{c \to f}^* + \frac{\Delta t}{\rho_f^{n+1/2}} (\sigma \kappa \nabla c)_f^{n+1/2} \quad (22) $$

where the subscript $c \to f$ indicates the linear interpolation of the variable from the cell-centers to the cell-faces and vice-versa ($f \to c$).

Once the temporary velocity $\mathbf{u}^*$ has been computed, the *Projection Step* allows a pressure field to be obtained, which guarantees that the divergence at time $n+1$ respects the continuity equation. In this work, the projection step was modified including the volume expansion term, which makes the divergence of the velocity field non-null (Equation 8). The projection step reduces to a Poisson equation which is solved using a multigrid solver [32]:

$$ \nabla \cdot \left( \frac{1}{\rho^{n+1/2}} \nabla p^{n+1} \right) = \frac{\nabla \cdot \mathbf{u}_f^*}{\Delta t} - \left[ \frac{\dot{m}}{\Delta t} \left( \frac{1}{\rho_g} - \frac{1}{\rho_l} \right) \delta_{\Gamma} \right]^{n-1/2} \quad (23) $$

$$ \mathbf{u}_f^{n+1} = \mathbf{u}_f^* - \frac{\Delta t}{\rho_f^{n+1/2}} \nabla p^{n+1} \quad (24) $$

where $\delta_{\Gamma}$ is approximated as the ratio between the discrete interfacial area (obtained from the PLIC reconstruction) and the volume of the cell. After the projection step, the collocated velocity is reconstructed, adding the acceleration terms and the pressure gradient:

$$ \mathbf{u}^{n+1} = \mathbf{u}^* + \frac{\Delta t}{\rho^{n+1/2}} \left[ (\sigma \kappa \nabla c)_f^{n+1/2} - \nabla p^{n+1} \right]_{f \to c} \quad (25) $$

The surface tension force $\mathbf{f}_{\sigma} = \sigma \kappa \nabla c$ is computed using a combination of the height-functions method for the curvature $\kappa$ and a well-balanced discretization of the surface tension force and the pressure gradient, which allows the equilibrium solution for a static droplet (Laplace equation) to be recovered [27, 33]. In this study we neglected the additional pressure jump caused by the phase change (recoil pressure). According to Mialhe et al. [34], this term enforces momentum conservation by counteracting the momentum loss due to the phase change. While this term is important in a momentum conserving scheme, its contribution does not affect the results discussed in this paper.
### 3.3. Advection Velocity

One of the main problems in modelling the phase change using VOF is obtaining an extended velocity, representative of the velocity of the liquid phase. When the momentum and continuity equations are discretized using a one-field formulation, the interfacial cells are affected by the Stefan flow contribution. This implies that the gas phase velocity interferes with the transport of the liquid volume fraction in the vicinity of the interface. To overcome this problem, different solutions were proposed in the literature [14, 15, 17, 18, 19, 31], although a unique solution valid for every system with phase change, has not been found yet. The most recent approaches are based on the solution of an additional Poisson equation for a velocity potential $\phi$:

$$ \nabla \cdot \left( \frac{1}{\rho} \nabla \phi \right) = \frac{\dot{m}}{\Delta t} \left( \frac{1}{\rho_g} - \frac{1}{\rho_l} \right) \delta_{\Gamma} \quad (26) $$

The potential is used to isolate the Stefan velocity $\mathbf{u}_S$ contribution:

$$ \mathbf{u}_S = -\frac{\Delta t}{\rho} \nabla \phi \quad (27) $$

Thus, the extended velocity field can be obtained by subtracting the Stefan velocity from the field velocity, resulting in an extended velocity field which is divergence-free by construction [18]:

$$ \mathbf{u}_E = \mathbf{u} - \mathbf{u}_S \quad (28) $$

Different variants of this method were adopted. Scapin et al. [17] solves Equation 26 over the whole domain, while Malan et al. [18] solves the same equation in a narrow band of cells close to the interface. Palmore et al. [19] combines PDE based Aslam extrapolation techniques [35], to enforce the continuity of the liquid velocity, with the solution of Equation 26 to correct errors in the velocity divergence after the extrapolation procedure. A different and more recent approach was proposed by Gennari et al. [31], who shift the expansion term in the Projection step toward the pure liquid or pure gas cells close to the interface. This approach does not require the solution of an additional Poisson equation, and it proved good performances for boiling simulations [31].

A common problem of all these methods is that they lead to nonphysical interface deformation for static droplet evaporation problems with strong
density ratio. To overcome this problem, we propose a new method based on a combination between the shifting approach proposed by Gennari et al. [31] and a novel Double Pressure-Velocity Coupling technique. According to this approach, the expansion term in Equation 8 is shifted toward the closest pure cells, such that the resulting velocity **u** contains a Stefan flow, which is continuous across the interfacial cells. After this step, a second set of Navier-Stokes equations is solved, without including the expansion term:

$$ \frac{\partial}{\partial t} \int_V \rho \mathbf{u}_E dV + \oint_S \rho \mathbf{u}_E (\mathbf{u}_E \cdot \mathbf{n}) dS = \oint_S \boldsymbol{\tau}_E \cdot \mathbf{n} dS + \int_V \rho \mathbf{g} dV + \int_V \mathbf{f}_\sigma dV \quad (29) $$

$$ \oint_S \mathbf{u}_E \cdot \mathbf{n} dS = 0 \qquad (30) $$

This set of Navier-Stokes equations is solved to find the extended velocity $\mathbf{u}_E$, which is used for the interface advection process. These equations are solved over the entire domain and using the same boundary conditions used for the velocity and pressure $\mathbf{u}$ and $p$. The same time-staggered discretization approach is adopted:

$$ \rho^{n+1/2} \left( \frac{\mathbf{u}_E^* - \mathbf{u}_E^n}{\Delta t} + \mathbf{u}_E^{n+1/2} \cdot \nabla \mathbf{u}_E^{n+1/2} \right) = \nabla \cdot \left( 2\mu^{n+1/2} \mathbf{D}_E^* \right) \quad (31) $$

$$ \mathbf{u}_{E,f}^* = \mathbf{u}_{E,c \to f}^* + \frac{\Delta t}{\rho_f^{n+1/2}} (\sigma \kappa \nabla c)_f^{n+1/2} \quad (32) $$

$$ \nabla \cdot \left( \frac{1}{\rho^{n+1/2}} \nabla p_E^{n+1} \right) = \frac{\nabla \cdot \mathbf{u}_{E,f}^*}{\Delta t} \quad (33) $$

$$ \mathbf{u}_{E,f}^{n+1} = \mathbf{u}_{E,f}^* - \frac{\Delta t}{\rho_f^{n+1/2}} \nabla p_E^{n+1} \quad (34) $$

$$ \mathbf{u}_E^{n+1} = \mathbf{u}_E^* + \frac{\Delta t}{\rho^{n+1/2}} \left[ (\sigma \kappa \nabla c)_f^{n+1/2} - \nabla p_E^{n+1} \right]_{f \to c} \quad (35) $$

At the beginning of the simulation, $\mathbf{u}$ and $\mathbf{u}_E$ are initialized to the same values. During the time integration loop, the two velocities are resolved independently but using the same discretization methods. The only difference is in the projection step: while the projection for the velocity $\mathbf{u}$ contains the volume expansion term, shifted toward the closest pure cells (Equation 23), the projection step for the velocity $\mathbf{u}_E$ enforces the divergence-free constraint
(Equation 33). The integration order is summarized in Algorithm 1. With this approach, we obtain two different velocity fields, **u**, and **u**<sub>E</sub>, both continuous and divergence-free across the interface, which are used for the interface advection and for the transport of the scalar fields. This method decouples the effect of the Stefan flow from the interface advection, and it proved to work well in situations where the Stefan velocity is much larger than the background velocity (i.e. static droplet cases with strong density ratio).

### 3.4. Transport of Scalar Fields

The transport equations for the scalar fields are solved using an operator splitting approach. First, the solution of the advection process is performed, then the diffusion is solved by taking into account the phase change contribution. According to the two-field approach implemented in this work, it is convenient to define the tracer form (tr) of the scalar fields (s) in the liquid and in the gas phase: tr<sub>l</sub> in s<sub>l</sub> and tr<sub>g</sub> in s<sub>g</sub>, respectively.

$$ 
\begin{aligned}
tr_l &= s_l c \\
tr_g &= s_g (1 - c)
\end{aligned} 
\qquad (36) $$

The values of the tracer fields are exactly the same as the scalar field values except at the interfacial cells. Using two separate tracers, numerical diffusion is limited and it is easy to switch from the scalar to the tracer field form and vice-versa. During this conversion, cells with small volume fraction values (c < ϵ) are considered to have zero volume fraction and the associated tracer fields are set to zero in those cells as well; the same tolerance used by Wenzel and Arienti [36] (ϵ < 1 × 10<sup>-10</sup>), was used here as this value is sufficiently small to avoid mass conservation issues. For post-processing purposes, a one-field formulation can be simply recovered as the sum of the tracer fields in the two phases: tr = tr<sub>l</sub> + tr<sub>g</sub>.

#### 3.4.1. Convective Transport

The convective part of Equation 11 and 14 is resolved by integrating the equation over a single cell of the domain, assuming constant physical properties and using the volume fraction c to correct the integration volume, in order to consider only the volume of the phase k of the scalar field which is being solved:

$$ \frac{\partial}{\partial t} \int_V c_k s_k dV = - \oint_S c_k s_k \mathbf{u}_k \cdot \mathbf{n} dS \qquad (37) $$
This equation is discretized using the same directionally-split procedure used for the volume fraction advection equation (Section 3.1). Following this approach, the value of the scalar field is updated considering the convective transport using the following equations:

$$ \frac{c_k^* - c_k^{n-1/2}}{\Delta t} = -\frac{1}{V} \sum_{f=1}^{NF} F_f(\mathbf{u}_k^n) \quad (38) $$

$$ \frac{(c_k s_k)^* - (c_k s_k)^{n-1/2}}{\Delta t} = -\frac{1}{V} \sum_{f=1}^{NF} s_{k,f}^n F_f(\mathbf{u}_k^n) \quad (39) $$

where $F_f$ is the volume flux of the phase being solved across the generic face $f$ and $s_{k,f}^n$ is the value of the scalar field on the generic face of the cell. Such value is obtained using the second-order upwind Bell-Colella-Glaz scheme [28]. The volume fluxes are computed using the same geometric fluxes calculation method used for the VOF field (Section 3.3). For droplet evaporation problems, the gas phase fluxes are computed using the velocity $\mathbf{u}$, which contains the Stefan flow, while the liquid phase fluxes exploit the extended velocity $\mathbf{u}_E$. After the advection step, the value of the scalar field is reconstructed in each cell as:

$$ s_k^* = \frac{(c_k s_k)^*}{c_k^*} \quad (40) $$

### 3.4.2. Diffusive Transport and Phase Change

The diffusion process is solved including the source terms that account for the phase change. In the context of the finite volume discretization, the diffusion part of the species equation (Equation 11) is re-written in integral formulation, neglecting the advection term and assuming that the material properties are constant:

$$ \frac{\partial}{\partial t} \int_V c_k \omega_{i,k} dV = \oint_S c_k \mathcal{D}_{i,k} \nabla \omega_{i,k} \cdot \mathbf{n} dS - \int_V \frac{\dot{m}_i}{\rho_k} \delta_\Gamma dV \quad (41) $$

This equation is discretized using a Backward Euler scheme for the time derivative, while the face gradients are approximated using the finite difference between two consecutive cells sharing the same face:

$$ \frac{c_k \omega_{i,k}^{n+1/2} - c_k \omega_{i,k}^*}{\Delta t} = \frac{1}{V} \sum_{f=1}^{NF} c_k \mathcal{D}_{i,k} \nabla \omega_{i,k} \Big|_f^{n+1/2} \cdot \mathbf{n}_f S_f - \left( \frac{\dot{m}_i}{\rho_k} \delta_\Gamma \right)^{n-1/2} \quad (42) $$
<figure><img src="image_5.png" /><figcaption>Figure 5: Discretization of the diffusive term with phase change.</figcaption></figure>

This equation is rearranged into a Poisson-Helmholtz equation and is solved in an implicit manner using a multigrid solver [32], as follows:

$$c_k^{n+1/2} \frac{\omega_{i,k}^{n+1/2} - \omega_{i,k}^*}{\Delta t} = \frac{1}{V} \sum_{f=1}^{NF} c_k \mathcal{D}_{i,k} \nabla \omega_{i,k}|_f^{n+1/2} \cdot \mathbf{n}_f S_f + \left(\frac{\dot{m}}{\rho_k \delta_\Gamma}\right)^{n-1/2} \omega_{i,k}^{n+1/2} - \left(\frac{\dot{m}_i}{\rho_k \delta_\Gamma}\right)^{n-1/2} \quad (43)$$

where the phase change contribution is included in the last two terms on the RHS. The source terms describe the amount of species i that evaporates $\dot{m}_i$ (explicit source), and the volume variations due to the phase change (implicit source). A limiting time step may be required for the stability because the solution of the diffusion equation is not fully implicit due to the explicit source term. The accurate solution of Equation 43 requires the surface fraction value $c_k|_f$, which is geometrically calculated from the PLIC interface reconstruction [37]. Using this approach the diffusion process is confined to the specific phase being solved and the diffusive flux across empty faces is simply neglected without introducing artificial diffusion (Figure 5).

The same approach is used for the discretization of the diffusion part of the temperature equation (Equation 14), resulting in the following form:

$$c_k^{n+1/2} \frac{T_k^{n+1/2} - T_k^*}{\Delta t} = \frac{1}{V} \left[ \sum_{f=1}^{NF} c_k \frac{k_f}{\rho_k c_{p,k}} \nabla T_k|_f^{n+1/2} \cdot \mathbf{n}_f S_f + \left( \frac{k_f}{\rho_k c_{p,k}} \frac{\partial T_k}{\partial \mathbf{n}_f} \Big|_k \cdot \mathbf{n}_f S_\Gamma \right)^{n-1/2} \right] \quad (44)$$

which can be solved given the interface temperature and the temperature gradients along the interface, computed as explained in the next section.
### 3.5. Vaporization Rate from the Interface Jump Condition

The interface properties (temperature and mass fractions) and the evaporation rate for every chemical species can be quantified by solving the interface jump condition via a non-linear system of equations (Equation 15) in every interfacial cell. A strategy for the decoupled solution of the non-linear system was developed for this purpose, in order to find reliable first guess values, which are then used by the root finding algorithm to solve the non-linear system. The system can be resolved in a decoupled manner by assuming that the evolution of one of the two phases is typically slower. During the evaporation of a liquid droplet, the evolution of the species mass fractions in the liquid phase is slower than those of the gas phase, due to the lower value of diffusivity. Therefore, the liquid interface mass fractions and temperature are assumed to be equal to the value of the field at the interfacial cells. The procedure for the decoupled solution of the system is reported below. The vaporization rate for each chemical species and the interfacial properties are computed at the beginning of the time step, thus using the volume fraction and scalar fields at time level $n - 1/2$.

1. In each interfacial cell, the interface temperature and liquid mass fractions are computed as:

$$ \hat{\omega}_{i,l} = \frac{\sum_{stencil}^{3x3} c\omega_{i,l}}{\sum_{stencil}^{3x3} c} \qquad (45) $$

$$ \hat{T} = cT_l + (1-c)T_g \qquad (46) $$

where $\hat{\omega}_{i,l}$ is the average, weighted on the volume fraction field, of the liquid mass fraction in the 3×3 stencil around the interfacial cell.

2. The interface mole fractions $\hat{x}_{i,l}$ are computed from the mass fractions, and the thermodynamic equilibrium constant is computed using the interface temperature from the previous point, obtaining the interface mole fractions in the gas phase:

$$ \hat{x}_{i,g} = K_{eq}(\hat{T}) \hat{x}_{i,l} \qquad (47) $$

Finally, the interface mass fractions in gas phase $\hat{\omega}_{i,g}$ are reconstructed.

3. Once the interface mass fractions are known, the total evaporation rate can be computed from mass balances at the gas-liquid interface,
summing over the total number of species in the liquid phase (NLS):

$$ \sum_{i=1}^{NLS} \dot{m} \hat{\omega}_{i,l} - \sum_{i=1}^{NLS} \rho_l \mathcal{D}_{i,l} \left. \frac{\partial \omega_{i,l}}{\partial \mathbf{n}_{\Gamma}} \right|_l = \sum_{i=1}^{NLS} \dot{m} \hat{\omega}_{i,g} - \sum_{i=1}^{NLS} \rho_g \mathcal{D}_{i,g} \left. \frac{\partial \omega_{i,g}}{\partial \mathbf{n}_{\Gamma}} \right|_g \quad (48) $$

where the sum over all the liquid species, of the diffusive fluxes in the liquid phase must be equal to zero by definition [24], while the sum of the mass fractions in the same phase must close to 1. Therefore, an explicit equation for the total evaporation rate is obtained:

$$ \dot{m} = - \frac{\sum_{i=1}^{NLS} \rho_g \mathcal{D}_{i,g} \left. \frac{\partial \omega_{i,g}}{\partial \mathbf{n}_{\Gamma}} \right|_g}{1 - \sum_{i=1}^{NLS} \hat{\omega}_{i,g}} \quad (49) $$

which allows the explicit calculation of the vaporization rate for each chemical species:

$$ \dot{m}_i = \dot{m} \hat{\omega}_{i,g} - \rho_g \mathcal{D}_g \left. \frac{\partial \omega_{i,g}}{\partial \mathbf{n}_{\Gamma}} \right|_g \quad (50) $$

4. Finally, the mass fraction of the inert (gas-only) species is adjusted and the interface temperature value can be refined using a root finding algorithm, from the knowledge of the evaporation rate:

$$ f(\hat{T}) = \sum_{i=1}^{NLS} \dot{m}_i \Delta h_{ev,i} - k_l \left. \frac{\partial T_l}{\partial \mathbf{n}_{\Gamma}} \right|_l - k_g \left. \frac{\partial T_g}{\partial \mathbf{n}_{\Gamma}} \right|_g = 0 \quad (51) $$

Using this procedure, the interface jump condition values are obtained. This first guess solution is then provided to the non-linear system of equations, which refines these values with the coupled solution. This approach allows the transport phenomena at the interface to be decoupled from the rest of the domain. The information obtained from the jump condition is then included in the solution of the transport equations in the form of source terms.

The interface gradient calculation must be carried out accurately for the correct solution of the interface jump condition. In this work, the interface gradients are computed exploiting methods developed for complex solid boundaries in Cartesian grids, like the Embedded Boundary Method by Johansen and Colella [38]. Following this approach, the gas-liquid interface is approximated as a solid boundary, taking advantage of the PLIC interface
<figure><img src="image_6.png" /><figcaption>Figure 6: Embedded Boundary Method for the calculation of the interface gradients. Adapted from [38].</figcaption></figure>

reconstruction, and the interface gradient is computed using a VOF-averaged 2nd order scheme, as proposed by Fleckenstein and Bothe [39]:

$$ \frac{\partial f}{\partial \mathbf{n}_{\Gamma}} = \left( c \frac{\hat{f} - f_1}{d_1} + (1-c) \frac{\hat{f} - f_2}{d_2} \right) \qquad (52) $$

where $\hat{f}$ is the value a generic scalar field at the interface, while $f_1$ and $f_2$ are the values of the same field in two points along the interface normal (Figure 6), interpolated from neighboring cells using bi-quadratic interpolation [31].

### 3.6. Time Stability and Solution Overview

The numerical solution algorithm is summarized in Algorithm 1. First, the vaporization rate is computed solving the interface jump condition at the initial time-step, in order to obtain $\dot{m}_i^{n-1/2}$ and every other interface property. Using that information, and with the velocity field known at time $n$, the VOF interface advection is performed in order to advance the interface to the next time level $c^{n+1/2}$, and to update the material properties by means of an arithmetic average. The species and temperature equations are solved using a two-field approach, and after the solution of the scalar transport equations, the velocity and pressure fields are obtained. The time step for the numerical solution is controlled by the time step required for the stability of the time-explicit discretization of the surface tension force [33]. However,
due to the presence of explicit source terms in the species and temperature equations, it may be necessary to further reduce the time step in case of strong evaporation rates.

<caption>**Algorithm 1: Summary of the solution algorithm**</caption>

Initialization of liquid volume fraction $c$, and velocity $\mathbf{u}$, $\mathbf{u}_E = \mathbf{u}$
Initialization of material properties $\rho$, $\mu$;
Initialization of tracer fields: $c_l = c$, $c_g = 1 - c$, $(c_l\omega_{i,l})$, $(c_g\omega_{i,g})$, $(c_lT_l)$, $(c_gT_g)$;
**while** $t < t_{end}$ **do**

// 1. — Interface Jump Condition
Gather $c^{n-1/2}$, $\omega_{i,l}^{n-1/2}$, $\omega_{i,g}^{n-1/2}$, $T_l^{n-1/2}$, $T_g^{n-1/2}$;
**foreach** *interfacial cell do*
    Solve interface jump condition → $\hat{\omega}_{i,l}^{n-1/2}$, $\hat{\omega}_{i,g}^{n-1/2}$, $\hat{T}_l^{n-1/2}$, $\hat{m}_i^{n-1/2}$ (Eq. 15);
    Update the total vaporization rate → $\dot{m}^{n-1/2}$;
    Update the volume expansion term → $(\dot{m}(1/\rho_l - 1/\rho_g))^{n-1/2}$;
**end**
Shift the expansion term according to [31];

// 2. — Scalar Transport Equations
Set $\Delta t$;
Compute the interface advection velocity → $\mathbf{u}_\Gamma^n$ (Eq. 20);
Solve VOF transport using $\mathbf{u}_\Gamma^n \rightarrow c^{n+1/2}$ (Eq. 18, 19);
Advection step in liquid phase using $\mathbf{u}_E^n \rightarrow c_l^*$, $(c_l\omega_{i,l})^*$, $(c_lT_l)^*$ (Eq. 38, 39);
Advection step in gas phase using $\mathbf{u}^n \rightarrow c_g^*$, $(c_g\omega_{i,g})^*$, $(c_gT_g)^*$ (Eq. 38, 39);
Reset $c_l^{n+1/2} = c^{n+1/2}$ and $c_g^{n+1/2} = 1 - c^{n+1/2}$;
Diffusion step with phase change → $\omega_{i,l}^{n+1/2}$, $\omega_{i,g}^{n+1/2}$, $T_l^{n+1/2}$, $T_g^{n+1/2}$ (Eq. 43, 44);

// 3. — Navier-Stokes Equations
Update one-field properties → $\rho^{n+1/2}$, $\mu^{n+1/2}$ (Eq. 2);
Predictor step for $\mathbf{u}^*$ (Eq. 21) and for $\mathbf{u}_E^*$ (Eq. 31);
Projection step for $p^{n+1/2}$ (Eq. 23) and $p_E^{n+1/2}$ (Eq. 33);
Update velocity $\mathbf{u}^{n+1}$ and $\mathbf{u}_E^{n+1}$ (Eq. 25, 35);
**end**

## 4. Tests and Results

The numerical model described in the previous section was implemented in the Basilisk framework [21]. The code was selected because it implements an efficient adaptive mesh refinement system, and a well-balanced discretization of the surface tension, which is important for small droplets and bubbles. The model developed in this work was verified using benchmark test cases in simple configurations, which allow the numerical results to be compared with the analytical solutions. More complex configurations are compared
with benchmark numerical simulations, discussing the qualitative behavior of the system and the mass conservation. The different test cases are reported in the next sections, while all the codes developed in this work and the simulation setups are freely available in the Basilisk sandbox [40].

### 4.1. Fixed Flux Evaporation

A 2D liquid droplet is placed at the center of a square domain. The droplet evaporates with a fixed vaporization rate, and the simulation proceeds until the complete consumption of the liquid. A simple mass balance on the liquid phase allows the analytical solution of the droplet radius in time to be computed:

$$ \frac{dR}{dt} = -\frac{\dot{m}}{\rho_l} \qquad (53) $$

In this test case $\dot{m}$ is imposed such that the droplet is completely consumed in 4 s. The physical properties adopted for this test case are those reported by Malan et al. [18]: $\rho_l/\rho_g = 2$, $\dot{m}/\rho_l = 0.05$ m s⁻¹, $R_0 = 0.23$ m. Differently from Malan et al.[18], the surface tension force was removed completely in this simulation, to focus on the phase change and to avoid adjustments to the droplet shape. The square domain has unit length, and outflow boundary conditions are imposed, while gravity is neglected. The simulation was performed at five different levels of refinement, from 4 to 8, meaning that, for the case with level 4, the number of cells along each domain direction is equal to $2^4$. Figure 7 (a) shows the trend of the droplet volume in time. The simulation results provide a good approximation of the analytical solution, except for the coarsest level of refinement. A convergent trend is recovered, and the relative error, computed on the droplet volume, shows second order convergence rate. Figure 7 (c) shows the comparison between the PLIC interfaces and the analytical solutions of the droplet shape for three different simulation times (1 s, 2 s, 3 s). Despite the absence of the surface tension force, the liquid droplet is kept spherical throughout the entire lifetime. From this test case, we can conclude that the use of a double pressure-velocity coupling provides a reliable velocity for the VOF transport equation. This velocity does not contain discontinuities and, therefore, it does not lead to the interface smearing. Moreover, the introduction of the interface regression velocity (Equation 20) allows the droplet shrinking process to be simulated with second order accuracy in a single VOF-advection step. Malan et al. [18] do not report the convergence rate for this test case. However, both models work
<figure><img src="image_7.png" /><figcaption>Figure 7: Fixed flux droplet evaporation results. Temporal evolution of the droplet volume: comparison between numerical simulations at different levels of refinement and the analytic solution (a); convergence rate, computed on the relative error of the liquid volume (b); droplet sphericity at 1 s, 2 s, and 3 s of simulation time (c) (http://basilisk.fr/sandbox/ecipriano/run/fixedflux.c).</figcaption></figure>

well with this simulation, obtaining results that approximate the theoretical
solution well.

4.2. Stefan Problem

The Stefan problems are classic benchmark test cases that different au-
thors have used for the validation of their phase change models. The Stefan
problem considered here consists of the evaporation of a liquid plane, pro-
moted by a temperature gradient between the gas phase and the liquid phase
that remains at the saturation temperature. At the beginning of the simu-
lation, a superheated wall heats up the vapor layer, which becomes hotter
than the liquid phase, leading to phase change. The analytical solution for
this system was reported by [18], and it describes the thickness of the vapor
layer $\delta(t)$ in time, as well as the temperature profile. For this simulation, the
same physical properties used by Malan et al. [18] were adopted: $\rho_l = 958$ kg
m$^{-3}$, $\rho_g = 0.6$ kg m$^{-3}$, $\sigma = 0.059$ N m$^{-1}$, $\mu_l = 2.82 \times 10^{-4}$ Pa s, $\mu_g = 1.23 \times 10^{-5}$
Pa s, $k_l = 0.68$ W m$^{-1}$ K$^{-1}$, $k_g = 0.025$ W m$^{-1}$ K$^{-1}$, $cp_l = 4216$ J kg$^{-1}$ K$^{-1}$,
$cp_g = 2080$ J kg$^{-1}$ K$^{-1}$, $\Delta h_{ev} = 2.256 \times 10^6$ J kg$^{-1}$. The simulation is charac-
terized by a Jacob number equal to 29.84 and it is performed using a square
domain, with length 0.01 m (Figure 8). The initial value of the gas layer
thickness is 322.5 $\mu$m, which is selected such that at least one layer of vapor
cells is present at the beginning of the simulation, even for the coarsest level
of refinement. Symmetry boundary conditions are imposed to the top and
<figure><img src="image_8.png" /><figcaption>Figure 8: Temperature field map and detail of the velocity jump across the interface at 10s. The white line is the 0.5 contour of the volume fraction field.</figcaption></figure>

bottom sides of the domain; an inlet zero velocity condition is imposed on the right side of the domain, where the temperature is fixed to 383 K; while on the left wall outflow boundary conditions are imposed. The simulation is performed at three different levels of refinement: 5, 6, and 7.

Figure 8 shows the temperature field in the domain after a total simulation time of 10 s. The temperature profile varies in the gaseous layer from the boundary condition imposed on the superheated wall to the saturation temperature at the gas-liquid interface. The zoom on a domain region close to the interface displays the strong velocity jump caused by the local volume expansion. This velocity jump is correctly managed by the double pressure-velocity coupling technique, which allows the scalar fields and the volume fraction to be transported by velocities which are divergence-free at the interface. The profiles reported in Figure 9 show the comparison between the simulation results and the analytical solution at three different levels of refinement. The thickness of the vapor layer approximates the analytical solution well even for the coarsest grid. The analysis of the relative error on the vapor layer thickness for the different levels of refinement shows a relative error around 0.05% already at the coarsest level of refinement, and it tends to converge approximately with second order accuracy with increasing mesh resolution. The comparison with Malan et al. [18] shows that the relative errors obtained using the model described in this work are smaller and converge faster to the analytic solution, compared to the approximately first order convergence rate reported in [18], whose relative errors are: 0.64%, 0.39%, 0.23% for the levels of refinement 6, 7, and 8, respectively. The faster convergence rate can be
<figure><img src="image_9.png" /><figcaption>Figure 9: Thickness of the vapor interface (a); relative errors on the interface thickness (b); temperature profile in the gas layer (c) at different levels of refinement (http://basilisk.fr/sandbox/ecipriano/run/stefanproblem.c).</figcaption></figure>

due to the combination between the calculation of the advection velocity and a more accurate evaluation of the interface gradients.

### 4.3. Epstein-Plesset Problem

This test case describes the dissolution of a spherical bubble in an undersaturated liquid environment. The phase change process is driven by a gradient of chemical species concentration between the bubble interface and the liquid bulk. Pure diffusive conditions are considered, meaning that no convective transport in the gas phase due to the Stefan flow is taken into account. Epstein and Plesset [41] obtained the analytic solution for the bubble radius in time, under quasi-static conditions, meaning that the vaporization rate is evaluated neglecting the effect of the interface motion. Under this assumption, the concentration field satisfies the steady-state diffusion equation at any time instant. This approximation is justified if the characteristic time of the diffusion process is much shorter than that of the interface regression [41]. From a mass balance, the radius of the bubble can be obtained as:

$$ \frac{dR}{dt} = -MW \frac{\mathcal{D}(\hat{C} - C_{\text{bulk}})}{\rho_g} \left( \frac{1}{R} + \frac{1}{\sqrt{\pi D t}} \right) \quad (54) $$

where $\hat{C}$ and $C_{\text{bulk}}$ are the interface concentrations of the chemical species and the liquid bulk phase concentration, respectively. The concentration profile can be obtained from the analytical solution of a species diffusion equation in radial coordinates:

$$ C(r,t) = C_{\text{bulk}} + (\hat{C} - C_{\text{bulk}}) \frac{R(t)}{r} \text{erfc} \left( \frac{r - R(t)}{\sqrt{4Dt}} \right) \quad (55) $$
The simulation setup was borrowed from Farsoiya et al. [42], the density of the two phases is set to the same value equal to 1 kg m⁻³ in order to remove the volume expansion term from the continuity equation, and, therefore, the Stefan convection. The chemical species diffusivity is set to 1 m²s⁻¹, while the surface tension σ is equal to 0.1 N m⁻¹. At the gas phase side of the interface, the concentration is set to a fixed value equal to 0.8 mol m⁻³, and the Stanton number is equal to 8 × 10⁻⁴ [42]. The simulation was performed using a 2D axial-symmetric domain, initializing the bubble at the lower-left corner of the domain, using an adaptive grid at maximum level of refinement 10. The concentration is set to the bulk value (C<sub>bulk</sub> = 0) at the boundaries of the liquid phase. The results in Figure 10 (a) show the consumption dynamics of the vapor bubble, which is expressed as the normalized radius of the bubble in time. This trend approximates well the Epstein-Plesset solution during the entire bubble lifetime. The concentration profile in time, Figure 10 (b), is plotted by sampling the value of chemical species concentration in a single point of the domain, corresponding to R₀ + 0.2 m. From this profile, the validity of the quasi-static approximation can be observed, since the concentration in the liquid phase goes from the bulk value to the concentration value established by the phase change phenomena, in a very small transient. The comparison between this work and the simulation reported in [42] shows that both models are able to accurately predict the theoretical profiles for the bubble radius and the chemical species concentration, except for small deviations due to the square domain and the quasi-static assumption. The model proposed by Farsoiya et al. [42] proved to have better performance in terms of computational time, while the model described in this work required a smaller time step for the stability of the diffusion step (Section 3.4). Despite this drawback, the model developed here is a more general formulation: it can be used both for evaporation problems in which the Stefan convection plays a major role, and also in pure diffusive conditions, like for example this test case.

### 4.4. *Scriven Problem*

The Scriven test case analyzes a bubble expanding in a superheated environment. The driving force for the phase change process is a temperature gradient between the interface of the bubble and the bulk phase. The primary difference with respect to the Epstein-Plesset case, is that the solution proposed by Scriven [43] includes the Stefan convection. The analytical solution is reported in [43, 44] and it describes the evolution of the bubble
<figure><img src="image_10.png" /><figcaption>Figure 10: Normalized diameter profile (a), and chemical species concentration at the coordinate $R_0+0.2m$ (b) for the Epstein Plesset problem (http://basilisk.fr/sandbox/ecipriano/run/epsteinplesset.c).</figcaption></figure>

radius in time, as well as the temperature field.

This test case was simulated on a 2D axial-symmetric configuration, where the bubble with initial radius 1 mm, is placed on a square domain with length 12 mm (Figure 11). The simulation was run at four different levels of refinement, from 6 to 9. The adaptive mesh refinement is used to minimize the computational cost of the simulation, which was run for 0.5 s, corresponding to a final bubble radius that is equal to twice the initial radius. The initial radius cannot be imposed to zero, as in the analytic solution, because it would require the ability to simulate the bubble nucleation process, which is beyond the aim of this work. The simulation setup was borrowed from [44], and it is characterized by a large density ratio, $ρ_l/ρ_g = 1623$, and a Jacob number equal to 3. The following physical properties were used: $ρ_l = 958$ kg m$^{-3}$, $ρ_g = 0.59$ kg m$^{-3}$, $σ = 0.001$ N m$^{-1}$, $μ_l = 2.81 × 10^{-4}$ Pa s, $μ_g = 1.26 × 10^{-6}$ Pa s, $k_l = 0.6$ W m$^{-1}$ K$^{-1}$, $k_g = 0.026$ W m$^{-1}$ K$^{-1}$, $cp_l = 4216$ J kg$^{-1}$ K$^{-1}$, $cp_g = 2034$ J kg$^{-1}$ K$^{-1}$, $Δh_{ev} = 2.256 × 10^6$ J kg$^{-1}$.

Figure 11 (a) shows that, despite the displacement between the reference solution and the numerical results for the coarsest level of refinement, the bubble radius dynamics tend to converge to the analytic solution when increasing the mesh resolution. The relative error between the numerical results and the Scriven solution for the bubble radius converges with approximately second order accuracy. The temperature profile along the radial coordinate is plotted in Figure 12 (c), and it shows that the profile obtained from the numerical solution converges to the analytical profile as well. As reported by
<figure><img src="image_11.png" /><figcaption>Figure 11: Map of the temperature field and the velocity field for the Scriven test case. Results at the maximum level of refinement and at the last time shot of the simulation (http://basilisk.fr/sandbox/ecipriano/run/scrivenproblem.c).</figcaption></figure>

Tanguy et al. [44], obtaining a monotonically convergent solution for this test case depends on the method used to obtain the extended velocity for transporting the interface and on the correct solution of the temperature field. The errors obtained by Tanguy et al. [44] using the Level Set approach for their most effective method are equal to 9.5%, 1.4%, and 1.0%, for levels 7, 8, and 9, respectively. The comparison between those results and the results obtained using the model described in this work show that the relative errors on the bubble radius are comparable. Obtaining convergence on this test case without using computationally expensive extrapolation procedures is a notable achievement: the bubble is always maintained spherical, despite the low surface tension value, and a spherically-symmetric temperature profile as well as a smooth velocity field are obtained.

### 4.5. Isothermal Droplet Evaporation

The evaporation of a multicomponent liquid droplet is tested against the numerical results obtained by Pathak et al. [45]. These results consider a pure liquid droplet, which evaporates due to a gradient of chemical species concentration between the interface and the gas phase. The interface temperature is taken to be constant, and therefore the interface mass fraction and the vapor pressure at the gas-liquid interface are constant as well. The benchmark model proposed in [45] utilizes spherical symmetry and unsteady conditions. The gas and the liquid phase are described using two body-fitted
<figure><img src="image_12.png" /><figcaption>Figure 12: Plot of the bubble radius in time (a); relative error on the bubble radius (b); temperature profile along the radius of the domain (c) at four different levels of refinement (http://basilisk.fr/sandbox/ecipriano/run/scrivenproblem.c).</figcaption></figure>

grids, which are connected by the interface point. The additional transport due to the Stefan flow is included, and no *quasi-static* approximations are used.

A droplet with starting diameter of 0.4 mm is initialized at the lower-left corner of a square domain, with dimensions equal to twice the droplet diameter (Figure 13). This simulation was performed using the following physical properties: $ρ_l = 10 \text{ kg m}^{-3}$, $ρ_g = 1 \text{ kg m}^{-3}$, $σ = 0.01 \text{ N m}^{-1}$, $μ_l = 1 \times 10^{-4} \text{ Pa s}$, $μ_g = 1 \times 10^{-5} \text{ Pa s}$, $D_l = 2 \times 10^{-3} \text{ m}^2 \text{ s}^{-1}$. Although this test case was conceived for pure liquid droplets, in this work it was run using the multicomponent model with four different pseudo-species. Following this approach, each chemical species is initialized with mass fraction 0.25 in the liquid phase, and 0 in the gas phase. Assuming that all the species have the same physical properties, the behavior of the liquid droplet should be equivalent to that of a pure droplet.

The simulation setup included symmetry boundary conditions at the boundaries in contact with the liquid droplet. On the top and right boundaries, zero-pressure outflow conditions are used, while the chemical species' mass fractions are solved with homogeneous Dirichlet boundary conditions. The mass fraction at the gas-phase side of the interface is constant and equal to 0.667, and the simulation was run at four different levels of refinement, from 5 to 8. Figure 14 shows the comparison between the benchmark results and the numerical model described in this work. The results show that the dynamics of the square diameter decay differ substantially from the steady behavior described by the $d^2$ law, which indicates that the square diameters
<figure><img src="image_13.png" /><figcaption>Figure 13: Map of the mass fraction profiles at three different simulation time instants. The white line, representing the interface, is the 0.5 iso-contour of the volume fraction field (http://basilisk.fr/sandbox/ecipriano/run/pureisothermal.c).</figcaption></figure>

ter decay variation is constant. Three different regions can be identified: during the first transient ($t < 0.15 \times 10^{-4}$ s), the droplet experiences fast consumption, due to the strong gradient between the interface and the gas phase at the beginning of the simulation. Between $t > 0.15 \times 10^{-4}$ s and $t < 0.7 \times 10^{-4}$ s the concentration of mass fraction in gas phase increases decreasing the consumption rate, until a steady-state vaporization rate is approached at $t > 0.7 \times 10^{-4}$ s. The displacement between the results obtained with this model and the benchmark solution increases as the liquid droplet shrinks. This is not surprising considering that the number of cells per diameter of the droplet decreases as the droplet is consumed. In fact, the higher the level of refinement, the more the square diameter decay converges to the benchmark solution. The agreement between the results of this work and the profile obtained by Pathak et al. [45] is very good, and the transient of the droplet consumption is well-captured. Figure 14 (b) reports the species mass fraction profiles for the highest level of refinement at three different simulation time instants. The trend obtained from the numerical simulations is close to the profiles obtained by Pathak et al. [45]. Since the multicomponent model with pseudo-species is used, the mass fractions profiles plotted in Figure 14 (b) are the sum of the gas phase mass fractions of each evaporating species. Using this approach, the total mass fractions are consistent with those of the pure species.

### 4.6. Non-Isothermal Droplet Evaporation

A liquid droplet is simulated in a non-isothermal environment, meaning that the interface mass fraction is not constant as in Section 4.5, but it is
<figure><img src="image_14.png" /><figcaption>Figure 14: Temporal evolution of normalized droplet squared diameter (a); radius profile of the gas phase mass fraction field at three different simulation time instants (b). The lines are the numerical results from the model described in this work, while the markers are the benchmark results from Pathak et al. [45] (http://basilisk.fr/sandbox/ecipriano/run/pureisothermal.c).</figcaption></figure>

a function of the interface temperature. The droplet is pure, the system pressure is 28.6 bar, and gravity is neglected. The initial diameter of the droplet is 5 µm, while the domain length (2D axial-symmetric) is twice the droplet diameter, following the same setup reported in [45]. The complexity of the problem and the multitude of physics phenomena involved do not allow an analytic solution. Therefore, the simulation results using the model developed in this work were compared with the numerical benchmark results provided by Pathak et al. [45]. The simulation was run at three different levels of refinement, from 6 to 8. The physical properties selected for this simulation are chosen in order to mimic the properties of n-heptane in nitrogen at high pressure: $ρ_l = 626.7$ kg m⁻³, $ρ_g = 17.51$ kg m⁻³, $σ = 0.01$ N m⁻¹, $μ_l = 1 × 10⁻⁴$ Pa s, $μ_g = 1 × 10⁻⁵$ Pa s, $k_l = 0.1121$ W m⁻¹ K⁻¹, $D_g = 6.77 × 10⁻⁷$ m² s⁻¹, $k_g = 04428$ W m⁻¹ K⁻¹, $cp_l = 2505$ J kg⁻¹ K⁻¹, $cp_g = 1053$ J kg⁻¹ K⁻¹, $Δh_{ev} = 3.23 × 10⁵$ J kg⁻¹. The molecular weight of n-heptane and nitrogen are 100 and 29 kg kmol⁻¹, respectively. The initial temperature of the liquid droplet is 363 K, while the environment temperature is 563 K. N-heptane mass fraction is equal to 1 inside the droplet and zero in the gas-phase, where pure nitrogen is present at the beginning of the simulation. Symmetry boundary conditions are imposed at the left and bottom sides of the domain, while outflow boundary conditions for velocity and pressure are imposed on the top and right sides. The conditions for
<figure><img src="image_15.png" /><figcaption>Figure 15: Maps of the n-heptane mass fraction and temperature field (http://basilisk.fr/sandbox/ecipriano/run/c7pathak.c).</figcaption></figure>

n-heptane mass fraction on the top and right sides are Dirichlet boundary conditions, with mass fraction set to 0, while the gas phase temperature is set to a value which is constant and equal to the initial gas temperature. The thermodynamic equilibrium at the gas-liquid interface is computed using An-toine equation, with coefficients from the NIST database [46, 47], in order to be consistent with the thermodynamics used in [45] to evaluate vapor-liquid equilibrium. Figure 15 shows the evolution of the n-heptane mass fractions and the temperature field. At the beginning of the simulation, the vaporization rate is small because the liquid droplet is cold. After a short time, the heat exchanged between the environment and the liquid droplet heats up the liquid phase increasing the temperature and therefore the vaporization rate of the droplet, which starts to be consumed. Figure 16 shows the square diameter decay and the interface temperature profile throughout the simulation. The square diameter profile does not suddenly decay at the beginning of the simulation, as in the isothermal evaporation case (Figure 14), because in this case the interface mass fraction is not fixed and it increases with the droplet heating from the hot environment. As the simulation proceeds, the
<figure><img src="image_16.png" /><figcaption>Figure 16: Plots of the square diameter decay in time (a); the interface temperature in time (b); temperature profiles along the domain radius (c); mass fraction profiles along the domain radius (d) (http://basilisk.fr/sandbox/ecipriano/run/c7pathak.c).</figcaption></figure>

interface temperature increases until approaching the wet bulb temperature of the liquid drop, which results from the interplay between the heating from the environment and the evaporation process that cools down the interface. The simulation results obtained with the model developed in this work converge to the numerical results provided by Pathak et al. [45], for both the square diameter decay and the interface temperature. Figure 16 also reports the radial profiles of temperature and n-heptane mass fractions at three different simulation time-steps. Apart from the good agreement between the simulation and the benchmark results, it is observed that the temperature of the liquid droplet increases and it becomes constant throughout the radius of the droplet. At the same time, n-heptane interface mass fraction increases in time and it matches the benchmark numerical values, starting from the mass fraction inside the liquid droplet, which is equal to 1, the mass fraction drops to the interface mass fraction and it follows the benchmark profile in the gas phase.
### 4.7. *Static Multicomponent Droplet Evaporation*

The isothermal evaporation of a binary droplet is studied, assuming the droplet is made of two chemical species with the same physical properties but with different volatility, in order to focus on the different vaporization rates of the two species. Initially, the mass fractions of the two chemical species are equal to 0.5, while an inert compound is present in the gas phase. The initial droplet diameter is 400 µm, the domain length is 4 times the droplet diameter, gravity is neglected and the following physical properties were used: $ρ_g = 1$ kg m⁻³, $σ = 0.03$ N m⁻¹, $μ_l = 1 × 10^{-4}$ Pa s, $μ_g = 1 × 10^{-5}$ Pa s, $D_l = 4 × 10^{-6}$ m² s⁻¹, $D_g = 8 × 10^{-5}$ m² s⁻¹. The relative volatility of the two chemical species in the liquid phase is equal to 2. In particular, the ratio between the vapor pressure and the thermodynamic pressure is 0.8 and 0.4 for the light and for the heavy components, respectively. The density ratio is equal to 10, to minimize the computational time required for the simulation, allowing the study of the droplet consumption trend and the mass conservation at three different levels of refinement: 7, 8, and 9. The diffusivity in liquid phase is adjusted such that the product between $ρ_l$ and $D_l$ is representative of the value of a hydrocarbon like n-heptane. Symmetry boundary conditions are imposed on the left and bottom edges of the domain, while outflow boundary conditions for velocity, pressure, and mass fractions are imposed on the top and right sides.

Figure 17 shows the evolution of the light, heavy, and inert chemical species mass fractions in time. The light component starts to evaporate first, and its mass fraction inside the liquid droplet decreases preferentially close to the gas-liquid interface. The heavy species accumulates, increasing its mass fraction inside the liquid phase and, therefore, its vaporization rate increases as well in time. Unlike the pure isothermal droplet (Section 4.5), the interface mass fractions for the binary droplet are not constant, due to the variable composition of the liquid phase.

Figure 18 reports the droplet consumption dynamics. The square diameter decay evolves similarly to the isothermal pure droplet evaporation test case. The numerical solution, at three different levels of refinement, converges to the same solution. The markers in Figure 18 (a) do not represent the analytical solution, which cannot be obtained for this test case; instead, they plot the results obtained from the model developed in this work, but without the solution of the non-linear system of equations. In these conditions, it can be observed that the procedure to obtain first-guess values for the non-linear system gives almost the same solution obtained solving the non-linear system.
<figure><img src="image_17.png" /><figcaption>Figure 17: Evolution of the mass fractions and the grid refinement for the static binary droplet evaporation: light component (top-left), heavy component (top-right), inert species (bottom-right) (http://basilisk.fr/sandbox/ecipriano/run/staticbi.c).</figcaption></figure>

system at the interface. This demonstrates the ability of the algorithm to find the first-guess interface values to provide reliable results leading to the same consumption dynamics of the liquid droplet. Figure 18 also reports the mass conservation test. The evaporation model is said to be conservative if the amount of each chemical species that evaporates is recovered in the gas phase and removed from the liquid phase. In practice, the total amount of evaporated species:

$$m_{Evap,i} = \int_{t_0}^{t_F} dt \int_V \dot{m}_i \delta_{\Gamma} dV \qquad (56)$$

is compared with the mass of the same species in liquid phase:

$$m_{Liq,i} = \int_V \rho_l c \omega_{l,i} dV \qquad (57)$$

and with the mass of the species in gas phase:

$$m_{gas,i} = \int_V \rho_g (1-c) \omega_{g,i} dV + \int_{t_0}^{t_F} dt \oint_S (\rho_g \omega_{g,i} \mathbf{u}_f - \rho_g D_i \nabla \omega_{g,i}) \cdot \mathbf{n} dS \quad (58)$$

which comprises the integration of the convective and diffusive mass fluxes across the boundaries. Figure 18 (b) and (c) reports the ratio between the mass of the evaporated component i and the initial mass of the same component in the gas and in the liquid phase respectively. The plots show that the light component increases its mass in the gas phase quickly with respect to the heavy component, and the consumption from the liquid phase follows
<figure><img src="image_18.png" /><figcaption>Figure 18: Square diameter decay of the static binary droplet (a); ratio between the mass of each component in gas phase and the initial mass of the same species in liquid phase (b); ratio between the mass of each component in liquid phase and the initial mass of the same species in liquid phase (c). Where not specified, the numerical results are obtained at the highest level of refinement (http://basilisk.fr/sandbox/ecipriano/run/staticbi.c).</figcaption></figure>

the same dynamics. The difference between the evaporated mass $m_{Evap}$ and $m_{Liq}$ for the light component is estimated using a relative error, defined as follows.

$$ \epsilon_{rel} = \frac{|m_{Evap,light} - m_{Liq,light}|}{m_{Evap,light}} \qquad (59) $$

The mass conservation errors for the three levels of refinement are 0.186, 0.047, and 0.017. They improve with increasing mesh resolution.

## 4.8. Convective Multicomponent Droplet Evaporation

A binary droplet in forced convective conditions is studied. The droplet has initial diameter 400µm and it is placed on a square domain (2D) with length 36 times the droplet diameter. Initially, the droplet is placed at a distance of 6 times the initial diameter from the left wall, and 18 times the initial diameter from the bottom. The droplet is made of two chemical species with the same physical properties but with different volatility, and the same relative volatility used in Section 4.7 is adopted. The simulation was performed at four different levels of refinement from 9 to 12, using the quadtree discretization, which is crucial for the optimization of the computational time of this test case. Gravity is neglected, but a gas flowrate is injected from the left side of the domain, in order to have $Re = 160$ with the initial properties of the simulation. Therefore, inlet boundary conditions for velocity and pressure are imposed on the left side, while outflow boundary
<figure><img src="image_19.png" /><figcaption>Figure 19: Evolution of the light component mass fraction, and detail of the light and heavy component mass fraction inside the liquid droplet at three different time snapshots (http://basilisk.fr/sandbox/ecipriano/run/forcedbi.c).</figcaption></figure>

conditions are used for the right wall. This simulation was performed using the following physical properties: $\rho_l = 800 \text{ kg m}^{-3}$, $\rho_g = 5 \text{ kg m}^{-3}$, $\sigma = 0.073 \text{ N m}^{-1}$, $\mu_l = 1.138 \times 10^{-3} \text{ Pa s}$, $\mu_g = 1.78 \times 10^{-5} \text{ Pa s}$, $D_l = 1.4 \times 10^{-7} \text{ m}^2 \text{ s}^{-1}$, $D_g = 1.25 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$. It was run in a static reference frame until the liquid droplet reaches the right side of the domain (corresponding to almost 20% consumption of the initial droplet volume).

Figure 19 (a) shows the mass fraction of the light components at three different time instants. When the simulation starts, the light component evaporates faster and the evaporated species is transported by the forced convective flux toward the right side of the domain. The selected Reynolds number regime eventually leads to the formation of Von Kármán-like streets which can be observed from the transported mass fraction field. Figure 19 shows the evolution of the mass fraction fields for both species inside the liquid droplet. The non-spherically symmetric distribution of the chemical species can be predicted by the numerical model, observing vortical structures that induce recirculation inside the liquid phase. This phenomenon has a direct influence on the heat and mass transfer and on the consumption time of the droplet. This configuration does not have analytical or numerical benchmark solutions, but the trend of the square diameter decay in time (Figure 20 (a)) shows that the consumption profile tends to the same result with increasing mesh resolution. The mass conservation test is reported in Figure 20 (b), where the variation of the liquid mass in the domain, computed
<figure><img src="image_20.png" /><figcaption>Figure 20: Evolution of the square diameter decay at four different levels of refinement (a); comparison between the variation of liquid mass in the domain for each species (solid line) and the total amount of vaporizing species (markers) normalized by the initial liquid mass (b) at the maximum level of refinement (http://basilisk.fr/sandbox/ecipriano/run/forcedbi.c).</figcaption></figure>

as explained in the previous section, is compared with the total amount of vaporizing mass. This analysis is performed for the two chemical species in the liquid phase and for the total liquid mass, obtaining good mass conservation for all three cases.

## 5. Conclusions

This paper presents a novel numerical model for multicomponent phase change in VOF simulations. The model describes systems with multiple chemical species in non-isothermal environments, considering the Stefan flow. We introduced a novel approach for the solution of the interface jump condition, and geometric discretization of the scalar transport equations. The volume expansion due to the phase change is tackled by proposing a combination between the expansion term shifting and a Double Pressure Velocity Coupling. The implemented numerical methods were validated using classic phase change problems: fixed flux evaporation, Stefan problem, and Scriven problem, obtaining second order convergence rates. The model is general, and it allows the simulation of pure diffusive systems without the Stefan flow, like the Epstein-Plesset test case. More complex configurations, like the evaporation of multicomponent droplets, reproduce the benchmark re-
sults from 1D numerical models well. The mass conservation error is small,
and it convergences with grid refinement.

The advantage of this model is the ability to simulate evaporation of
mixtures including thermal effects, and to correctly manage the interface
velocity jump also for static droplets with strong density ratio. The code
developed in this work and the simulation setups are documented and re-
leased on the Basilisk sandbox [40], allowing other researchers interested in
phase change simulations to develop and test new methods starting from an
existing framework.

References

[1] B. Abramzon, W. A. Sirignano, Droplet vaporization model for spray combustion calculations, International journal of heat and mass transfer 32 (9) (1989) 1605–1618.

[2] C. Law, Recent advances in droplet vaporization and combustion, Progress in energy and combustion science 8 (3) (1982) 171–201.

[3] W. A. Sirignano, Fluid dynamics and transport of droplets and sprays, Vol. 2, Cambridge university press Cambridge, 2000.

[4] A. Cuoci, M. Mehl, G. Buzzi-Ferraris, T. Faravelli, D. Manca, E. Ranzi, Autoignition and burning rates of fuel droplets under microgravity, Combustion and Flame 143 (3) (2005) 211–226.

[5] C. W. Hirt, B. D. Nichols, Volume of fluid (vof) method for the dynamics of free boundaries, Journal of computational physics 39 (1) (1981) 201–225.

[6] L. Deike, W. K. Melville, S. Popinet, Air entrainment and bubble statistics in breaking waves, Journal of Fluid Mechanics 801 (2016) 91–129.

[7] D. Fuster, A. Bagué, T. Boeck, L. Le Moyne, A. Leboissetier, S. Popinet, P. Ray, R. Scardovelli, S. Zaleski, Simulation of primary atomization with an octree adaptive mesh refinement and vof method, International Journal of Multiphase Flow 35 (6) (2009) 550–565.

[8] D. Fuster, G. Agbaglah, C. Josserand, S. Popinet, S. Zaleski, Numerical simulation of droplets, bubbles and waves: state of the art, Fluid dynamics research 41 (6) (2009) 065001.
[9] S. W. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, Journal of computational physics 160 (2) (2000) 662–682.

[10] D. Juric, G. Tryggvason, Computations of boiling flows, International journal of multiphase flow 24 (3) (1998) 387–410.

[11] G. Son, V. Dhir, Numerical simulation of saturated film boiling on a horizontal surface, Journal of Heat Transfer 119 (3) (1997) 525–533.

[12] A. Esmaeeli, G. Tryggvason, Computations of film boiling. part i: numerical method, International journal of heat and mass transfer 47 (25) (2004) 5451–5461.

[13] A. Esmaeeli, G. Tryggvason, Computations of film boiling. part ii: multi-mode film boiling, International journal of heat and mass transfer 47 (25) (2004) 5463–5476.

[14] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source terms, Journal of Computational Physics 227 (11) (2008) 5871–5895.

[15] C. Kunkelmann, P. Stephan, Cfd simulation of boiling flows using the volume-of-fluid method within openfoam, Numerical Heat Transfer, Part A: Applications 56 (8) (2009) 631–646.

[16] J. Schlottke, B. Weigand, Direct numerical simulation of evaporating droplets, Journal of Computational Physics 227 (10) (2008) 5215–5237.

[17] N. Scapin, P. Costa, L. Brandt, A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows, Journal of Computational Physics 407 (2020) 109251.

[18] L. Malan, A. G. Malan, S. Zaleski, P. Rousseau, A geometric vof method for interface resolved phase change and conservative thermal energy advection, Journal of Computational Physics 426 (2021) 109920.

[19] J. Palmore Jr, O. Desjardins, A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows, Journal of Computational Physics 399 (2019) 108954.
[20] S. Zhao, J. Zhang, M.-J. Ni, Boiling and evaporation model for liquid-gas flows: A sharp and conservative method based on the geometrical vof approach, Journal of Computational Physics 452 (2022) 110908.

[21] S. Popinet, Basilisk c reference manual., http://basilisk.fr/Basilisk%20C (2013).

[22] J. Lopez-Herrera, A. Ganan-Calvo, S. Popinet, M. Herrada, A vof numerical study on the electrokinetic effects in the breakup of electrified jets, Int. J. Multiph. Flow 71 (2015) 14-22.

[23] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of gas-liquid multiphase flows, Cambridge university press, 2011.

[24] R. B. Bird, Transport phenomena, Appl. Mech. Rev. 55 (1) (2002) R1-R4.

[25] R. Rota, Fondamenti di termodinamica dell'ingegneria chimica, Pitagora, 2015.

[26] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries, Journal of computational physics 190 (2) (2003) 572-600.

[27] S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, Journal of Computational Physics 228 (16) (2009) 5838-5866.

[28] J. B. Bell, P. Colella, H. M. Glaz, A second-order projection method for the incompressible navier-stokes equations, Journal of computational physics 85 (2) (1989) 257-283.

[29] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual review of fluid mechanics 31 (1) (1999) 567-603.

[30] G. D. Weymouth, D. K.-P. Yue, Conservative volume-of-fluid method for free-surface simulations on cartesian-grids, Journal of Computational Physics 229 (8) (2010) 2853-2865.

[31] G. Gennari, R. Jefferson-Loveday, S. J. Pickering, A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows, Chemical Engineering Science 259 (2022) 117791.
[32] S. Popinet, A quadtree-adaptive multigrid solver for the serre-green-naghdi equations, Journal of Computational Physics 302 (2015) 336-358.

[33] S. Popinet, Numerical models of surface tension, Annual Review of Fluid Mechanics 50 (2018) 49-75.

[34] G. Mialhe, S. Tanguy, L. Tranier, E.-R. Popescu, D. Legendre, An extended model for the direct numerical simulation of droplet evaporation. influence of the marangoni convection on leidenfrost droplet, Journal of Computational Physics 491 (2023) 112366.

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

[36] E. A. Wenzel, M. Arienti, A conservative, interface-resolved, compressible framework for the modeling and simulation of liquid/gas phase change, Journal of Computational Physics 477 (2023) 111957.

[37] Q. Magdelaine-Guillot de Suduiraut, Hydrodynamique des films liquides hétérogènes, Ph.D. thesis, Sorbonne université (2019).

[38] H. Johansen, P. Colella, A cartesian grid embedded boundary method for poisson's equation on irregular domains, Journal of Computational Physics 147 (1) (1998) 60-85.

[39] S. Fleckenstein, D. Bothe, A volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes, Journal of Computational Physics 301 (2015) 35-58.

[40] E. Cipriano, Code repository., http://basilisk.fr/sandbox/ecipriano/ (2023).

[41] P. S. Epstein, M. S. Plesset, On the stability of gas bubbles in liquid-gas solutions, The Journal of Chemical Physics 18 (11) (1950) 1505-1509.

[42] P. K. Farsoiya, Q. Magdelaine, A. Antkowiak, S. Popinet, L. Deike, Direct numerical simulations of bubble-mediated gas transfer and dissolution in quiescent and turbulent flows, Journal of Fluid Mechanics 954 (2023) A29.
[43] L. Scriven, On the dynamics of phase growth, Chemical engineering science 10 (1-2) (1959) 1-13.

[44] S. Tanguy, M. Sagan, B. Lalanne, F. Couderc, C. Colin, Benchmarks and numerical methods for the simulation of boiling flows, Journal of Computational Physics 264 (2014) 1-22.

[45] A. Pathak, M. Raessi, Steady-state and transient solutions to drop evaporation in a finite domain: Alternative benchmarks to the d2 law, International Journal of Heat and Mass Transfer 127 (2018) 1147-1158.

[46] G. F. Carruth, R. Kobayashi, Vapor pressure of normal paraffins ethane through n-decane from their triple points to about 10 mm mercury, Journal of Chemical and Engineering Data 18 (2) (1973) 115-126.

[47] D. Rossini, Vapor pressures and boiling points of some paraffin, alkylcyclopentane, alkylcyclohexane, and alkylbenzene hydrocarbons, Journal of Research of the National Bureau of Standards 35 (1945).