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

HAL
open science

# A Vaporization Model for Continuous Surface Force Approaches and Subcooled Configurations

Charles Brissot, Léa Cailly-Brandstätter, Elie Hachem, Rudy Valette

► To cite this version:

Charles Brissot, Léa Cailly-Brandstätter, Elie Hachem, Rudy Valette. A Vaporization Model for Continuous Surface Force Approaches and Subcooled Configurations. Fluids, 2023, 8 (8), pp.233. 10.3390/fluids8080233. [hal-04384092](https://hal-04384092)

HAL Id: hal-04384092

https://minesparis-psl.hal.science/hal-04384092v1

Submitted on 8 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Vaporization Model for Continuous Surface Force Approaches and Subcooled Configurations

Charles Brissot, Léa Cailly-Brandstätter, Elie Hachem and Rudy Valette *

MINES ParisTech, PSL Research University, CEMEF—Centre for Material Forming, CNRS UMR 7635, CS 10207
Rue Claude Daunesse, 06904 Sophia-Antipolis Cedex, France

* Correspondence: rudy.valette@minesparis.psl.eu

**Abstract:** The integration of phase change phenomena through an interface is a numerical challenge that requires proper attention. Solutions to properly ensure mass and energy conservation were developed for finite difference and finite volume methods, but not for Finite Element methods. We propose a Finite Element phase change model based on an Eulerian framework with a Continuous Surface Force (CSF) approach. It handles both momentum and energy conservation at the interface for anisotropic meshes in a light an efficient way. To do so, a model based on the Level Set method is developed. A thick interface is considered to fit with the CSF approach. To properly compute the energy conservation, heat fluxes are extended through this interface thanks to the resolution of a transport equation. A dedicated pseudo compressible Navier-Stokes solver is added to compute velocity jumps with a source term at the interface in the velocity divergence equation. Several 1D and 2D benchmarks are considered with increasing complexity to highlight the performances of each feature of the framework. This stresses the capacity of the model to properly tackle phase change problems.

**Keywords:** phase change; continuous surface force; level set; heat flux jump computation; velocity jump; subcooling; Stefan problem

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

check for
updates

**Citation:** Brissot, C.;
Cailly-Brandstätter, L.; Hachem, E.;
Valette, R. A Vaporization Model for
Continuous Surface Force Approaches
and Subcooled Configurations. *Fluids
2023*, *8*, 233. https://doi.org/
10.3390/fluids8080233

Academic Editors: Artur V. Dmitrenko
and Nilanjan Chakraborty

Received: 6 July 2023
Revised: 29 July 2023
Accepted: 8 August 2023
Published: 19 August 2023

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

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

## 1. Introduction

Boiling phenomena are a key subject for many industrial domains in need of efficient
heat exchangers. Most of the time, boiling is used as a great energy sink. The characteriza-
tion of heat transfer between solids and fluids is of high interest for all applications.

Theoretical models and empirical laws have been developed to assess this quantity. In
1934, Nukiyama [1] was the first to describe, as a whole, all the different modes of boiling
by plotting heat fluxes versus temperature, thanks to a heat wire that was plunged into
a bath full of water. Nucleation, transition boiling and film boiling were observed. Then,
many experiments using various fluids, conditions and parameters were made to better
understand each mode. Correlations based on physical arguments were developed to
estimate heat flux laws by pioneer works [2–7] as well as more recent investigations [8,9].
Some simple analytical models have also been developed [10]. All these approaches are very
interesting and provide analytical laws but are often case-dependent, and their precision
remains limited when confronted with industrial configurations.

Computational fluid dynamics (CFD) models have been proposed to better under-
stand the physics of phase change. CFD allows the reproduction and evaluation of local
phenomena that are difficult to estimate during tests. Once numerical models are validated,
they can even be used as a powerful experimental laboratory. However, the development of
a numerical model that could predicatively simulate all boiling modes remains a challenge.

Pioneer works have been conducted to simulate two-phase flows with phase change
in 2D. Lee [11] developed a semi-implicit two-fluid phase change model with a sharp
interface. A projection technique with an iterative scheme was used to compute the velocity
jump, and the mass transfer was computed thanks to a proportional law. This enabled the simulation of a steam-water jet impingement. Using the volume of fluid (VOF) method, this work was extended by Rattner et al. [12] for film condensation. The mass transfer rate previously chosen was arbitrarily assessed through enthalpy considerations. This model was embedded for a finite volume approach within the Open-FOAM library. More recently, Giustini et al. [13] proposed a similar approach with the added value that their framework could fit arbitrary meshes without the need to introduce smearing of the mass transfer term. Son et al. [14] combined this second-order projection method with a Level Set method to solve saturated 2D horizontal film boiling. Juric et al. [15] also simulated this problem with a single-field formulation and a front tracking method implemented with finite differences. Welch et al. [16] performed the same but with a volume-of-fluid-based interface tracking method. Once again, the projection method was used, and the reconstructed interface led to the computation of temperature gradients from both sides. Kang et al., Gibou et al. and Tanguy et al. [17–19] implemented a Ghost Fluid Method to ease the introduction of interfacial terms within an Eulerian finite difference method. A similar approach was conducted by Esmaeeli et al. [20], who eliminated the iterative procedure with a predictor-corrector algorithm. Using the interface tracking and projection methods with a sharp interface approach, Sato et al. [21] implemented a staggered finite volume vaporization solver and computed proper bubble nucleation and growth. Khalloufi et al. [22] worked on the complete simulation of quenching processes in 3D for simple geometries thanks to an Eulerian framework with Level Set methods and Finite Elements. Other innovative methods have been developed, like the coupled Level Set and volume of fluid method [23], aiming to improve the interface definition and mass change rate computation.

These recent developments, however, require a large number of numerical computations over the interface. These could be troublesome when dealing with industrial applications. This is the reason why the choice of the Continuous Surface Force (CSF) approach is often preferred, as it guarantees an integrated implementation that is very light. This feature is especially appreciated when dealing with industrial applications. In the present work, we propose a phase change model within an Eulerian framework based on a stabilized Finite Element method. These methods have been proven to support the Navier-Stokes equations as well as convection diffusion reaction numerical schemes, allowing them to overcome challenging test cases [24]. The Level Set method is used to track the interface. A CSF approach is applied to deal with stiff differences in the physical properties, with a dedicated method to properly handle thermal variations. This consideration has been shown to bring stability to multiphase flow solvers with a simple implementation [25–27].

Most of the time, when boiling takes place in industrial processes, the liquid located far from the vaporizing front is below its saturation temperature. The difference in temperature (usually called subcooling) has been shown to significantly improve heat transfers [9,28–33]. The entire boiling curve is shifted at a higher Leindenfrost temperature the minimum temperature for the stabilization of a vapor film. More generally, the physics of boiling drastically changes from the saturated case. This is the reason why proper consideration of heat flux jumps at the vapor/liquid interface is crucial.

Moreover, as the vapor phase is usually much lighter than the liquid phase, vaporization entails an important dilatation at the interface. This is numerically challenging and requires specific attention adapted to the considered framework. To do so, a dedicated pseudo-compressible Navier-Stokes solver is added, based on the work of Khalloufi et al. [22]. It computes velocity jumps with a volume source term at the interface in the velocity divergence equation. This ensures mass and momentum conservation at the interface in a light and efficient way, with consideration of volume mass transfer distribution. To reduce the computational time in prediction of future cases, it is coupled with a remeshing algorithm (see [34] for further details on the implementation of this algorithm).

To the authors' knowledge, there is not any existing method that accurately combines the CSF approach with the heat flux jump computation and a volume mass transfer distribution. In the present work, this is conducted with the consideration of a complete fictitious
thick interface for the mixing laws described in Section 2. The phase change distribution is also smoothed and the temperature profile is set to the phase change temperature at the interface. Thus, the jump is not computed on one straight surface that usually features the interface, but on the border of this thick interface. Such consideration is new and requires the design of a dedicated method that will be described in details in Sections 3 and 4. To validate the presented method for diffusion-only cases, simple 2D test cases are studied in Section 5. Then test case validating the pseudo-compressible approach are presented in Section 6. The complete vaporization framework is finally tested and validated in Section 7 on the compressible subcooled Stefan problem and on an adiabatic problem.

## 2. Phase Change Modeling

### 2.1. Mathematical Framework — Level Set Method

Two-phase flows with a mixture of liquid and vapor are considered. In order to separate the two different phases, the Level Set formalism is used. For $\vec{x}$ the position vector, a signed distance function $(\vec{x}, t) \mapsto \alpha$ describes the interface $I$ between the two phases:

$$ \alpha(\vec{x}, t) = \begin{cases} d(\vec{x}, I) & \text{in the vapor phase} \\ -d(\vec{x}, I) & \text{in the liquid phase} \\ 0 & \text{on the interface} \end{cases} \qquad (1) $$

This way, the gradient of the distance function $\vec{\nabla}\alpha$ defines the normal of the interface. The normal of the interface $\vec{n}$ is defined as $\frac{\vec{\nabla}\alpha}{|\vec{\nabla}\alpha|}$.

The associated characteristic function $\alpha \mapsto H_{\alpha}$ allows to determine the fraction of vapor and liquid at any point. To avoid numerical instabilities, an interface thickness $2\epsilon$ is considered. Different functions can be used to smooth this Heavyside function, such as a tangent or a sine function. This latter has been chosen in this work:

$$ H_{\alpha} = \begin{cases} 1 & \text{if } \alpha > \epsilon \\ \frac{1}{2}\left(1 + \frac{\alpha}{\epsilon} + \frac{1}{\pi}\sin\left(\frac{\pi\alpha}{\epsilon}\right)\right) & \text{if } |\alpha| \le \epsilon \\ 0 & \text{if } \alpha < -\epsilon \end{cases} \qquad (2) $$

The reasons behind this need of a smoothing function and a complete justification of the consistency of this solution can be found in [34]. This allows us to determine the characteristics of a fluid particle regarding the characteristics of the two phases.

Finally, to determine properties at the interface, the smoothed Dirac function $\alpha \mapsto \delta_{\alpha}$ is defined as the derivative of $H_{\alpha}$ with respect to $\alpha$. Thus:

$$ \delta_{\alpha} = \begin{cases} 0 & \text{if } |\alpha| > \epsilon \\ \frac{1}{2\epsilon}\left(1 + \cos\left(\frac{\pi\alpha}{\epsilon}\right)\right) & \text{if } |\alpha| \le \epsilon \end{cases} \qquad (3) $$

Two relations can be noticed:

$$ \vec{\nabla} H_{\alpha} = (\vec{\nabla}\alpha)\delta_{\alpha} \qquad (4) $$

$$ \frac{\partial H_{\alpha}}{\partial t} = \frac{\partial \alpha}{\partial t} \delta_{\alpha} \qquad (5) $$

### 2.2. Continuous Description of Phase Change

The mass change is represented by the vector $\vec{m}: |\vec{m}|$ is the absolute surface mass transfer rate from one phase to another. The local direction of $\vec{m}$ indicates if it stands for
vaporization or condensation. $\vec{\dot{m}}$ is assumed to be orthogonal to the surface, thus aligned with the normal vector $\vec{n}$. The total value of mass change over a surface S of the interface is $\int_S \vec{\dot{m}} \cdot \vec{n} \, dS$ (see Figure 1).

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic description of the mass transfer rate vector $\vec{\dot{m}}$. The red dashed line corresponds to the interface position after a time dt. On the **upper part**, the liquid vaporizes and $\vec{\dot{m}}$ points towards the vapor phase. On the **lower part**, the vapor condensates and $\vec{\dot{m}}$ points towards the liquid phase.</figcaption></figure>

A continuous mass change profile [22] is chosen. The mass change is spread on the volume of an interface of thickness $2\epsilon$ thanks to a Dirac distribution $\delta_\alpha$ (see Figure 2). The local infinitesimal mass transfer rate $d\dot{M}$ over an infinitesimal volume $dV$ then reads:

$$d\dot{M} = \vec{\dot{m}} \cdot (\vec{\nabla}\alpha) \delta_{\alpha} dV. \qquad (6)$$

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Modeling of the interface temperature profile adapted to a Continuous Surface Force approach. The temperature is fixed at $T_{sat}$ on the entire layer of thickness $2\epsilon$ to create a well defined heat flux jump. The mass transfer is spread out on the layer with a $\delta_\alpha$ profile. (a) Sharp interface configuration; (b) Smooth interface configuration.</figcaption></figure>

With this configuration, $\vec{\dot{m}}$ is supposed to be locally constant and orthogonal to the interface, and the total value of mass change over a surface S of the interface is conserved as:

$$\int_S \int_{[-\epsilon; \epsilon]} d\dot{M} = \int_S \int_{[-\epsilon; \epsilon]} \vec{\dot{m}} \cdot (\vec{\nabla}\alpha) \delta_{\alpha} d\alpha \, dS = \int_S \vec{\dot{m}} \cdot \vec{n} \, dS \qquad (7)$$

by definition of $\alpha$ and $\delta_\alpha$.
Working with a Continuous Surface Force approach, the notion of heat flux jump is ambiguous, as all quantities are continuous. Moreover, its value must be accessible for each node inside the interface to determine the local constant value of $\vec{m}$, and to be able to compute the velocity jump and the energy jump terms.

To address these issues, the present work proposes to set the entire interface temperature at $T_{sat}$, as shown in Figure 2. The heat flux jump is then considered from the two extremities of the interface. As only conductive heat fluxes are considered in the present study, the condition on $\vec{m}$ reads:

$$|\vec{m}|\mathcal{L} = \left[ k \nabla T \cdot \vec{n} \right]_{-\epsilon}^{\epsilon} \qquad (8)$$

where $k$ is the conductivity.

The extension method and the algorithm relying on it to compute the balance of flux on each points of the interface are presented hereafter.

### 3. Heat Flux Jump Computation

#### 3.1. Description of the Extension Method

This method relies on a common technique explained by Aslam et al. [35]. The basis of this technique is to extend on $\Omega$ a scalar quantity $q_0$ that lives on $\Omega_0$ (generally $\Omega$ and $\Omega_0$ form a partition of the computation space). The extension is conducted by following streamlines of a normalized vector field $\vec{n}$ through an interface $\Gamma$ that separates the two domains $\Omega$ and $\Omega_0$. In other words, the solution $q$ equals $q_0$ on $\Omega_0$, and on $\Omega$, $q$ is set constant on the streamlines of $\vec{n}$, its value being the value of $q_0$ on the intersection of $\Gamma$ with the streamline. This corresponds to the resolution of the following system:

$$\left\{ \begin{array}{ll} q = q_0 & \text{on } \Omega_0 \\ \vec{n} \cdot \nabla q = 0 & \text{on } \Omega \end{array} \right. \qquad (9)$$

$$\left\{ \begin{array}{ll} q = q_0 & \text{on } \Omega_0 \\ \vec{n} \cdot \nabla q = 0 & \text{on } \Omega \end{array} \right. \qquad (10)$$

Equation (9) is solved immediately, and is used as a Dirichlet condition on $\Gamma$ for the resolution of (10). For frontiers of $\Omega$ that are at the vicinity of the calculation domain, a Neumann condition of null flux is used.

The second equation by itself is not numerically unconditionally stable, so a diffusion term is added:

$$\vec{n} \cdot \nabla q - \nabla \cdot (\lambda \nabla q) = 0 \qquad (11)$$

The value of $\lambda$ should be as small as possible in $\Omega$ to limit its impact on the solution, but large enough to stabilize the solver. If $h$ is the mesh size, a dimensional analysis leads to:

$$\frac{|\vec{n}|q}{h} \approx \frac{\lambda q}{h^2} \qquad (12)$$

As $\vec{n}$ is normalized, the condition reads $\lambda \approx h$. $\lambda$ will be chosen at most of the same order of magnitude as $h$.

Before its implementation on heat flux jump computation, this method is tested on a theoretical test case: the extension of a scalar field that has a cylindrical symmetry on a disc. The field $\vec{n}$ is based on the gradient of the Level Set $\alpha$ (the signed distance function in respect to the interface that is positive inside $\Omega$). As a result, $\vec{n}$ is normal to the interface.

#### 3.2. The Disc Case

$\Omega$ is a disc of radius $R = 0.2$ centered on the origin of a $1 \times 1$ square domain (see Figure 3a). The source field $q_0$ is defined in the cylindrical coordinates system $(r, \theta)$ such as:

$$q_0(r, \theta) = -2r \qquad (13)$$

$\vec{n}$ is defined as:

$$\vec{n}(r, \theta) = -e_r \qquad (14)$$

Thus $\vec{n}$ has a singularity at the origin.
The case is implemented on an unstructured mesh of triangles of characteristic size $h = 10^{-2}$ with $\lambda = 10^{-3}$. Results projected on the horizontal axis of the disc are plotted on Figure 3b.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** The disc test case: an initial radial scalar field is extended on a disc (the red zone) following the radial vector field (the black arrows). (a) Schematic description of the disc test case. (b) Source and extended field values on the radial axis $y = 0$.</figcaption></figure>

We can see that the extension of $q_0$ over $\Omega$ is well computed, and the singularity of the vector field $\vec{n}$ does not cause issues to the solver. Errors can be controlled by refining the mesh.

### 3.3. Determination of Heat Flux Jumps over a Smooth Interface

Back to the phase change problem, the idea is to apply this tool on each side of the interface to access the values of the balance of heat fluxes on every point of the interface.

The method is applied with the scalar field $\frac{\vec{\nabla}\alpha}{|\vec{\nabla}\alpha|} \cdot (k\vec{\nabla}T)$ to determine the heat fluxes on the vapor and on the liquid sides.

A first simple approach would be to choose $\Omega$ as the domain of the interface, and $\Omega_0$ the union of the vapor and the liquid domains for both cases (it is useless to further project the heat flux from one phase on the whole domain of the other phase).

However, due to diffusion effects, the transition from $T_{sat}$ on the interface to the temperature profile diffuses on a few elements. A "safety margin" is taken to make sure that the temperature profile is well established. This is the reason why the heat flux value to project is drawn one mesh size away from the isovalue $\alpha = \pm\epsilon$. Moreover, to further reduce the influence of the diffusion term $\lambda$, the end of $\Omega$ is considered three mesh sizes away from the isovalue $\alpha = \pm\epsilon$. Figure 4 summarizes these considerations. Doing so, the heat flux jump is well defined and constant on the domain $\{\alpha \in [\epsilon; \epsilon]\}$.

As an illustration and test purpose, this method is implemented on a simple 2D test case. A 1 × 1 square with a vertical interface $x = 0$ at its center (see Figure 5a).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Modeling of the interface temperature and heat flux profiles with a Continuous Surface Force approach. In order to well capture the heat flux jump, heat flux values are extended from a distance $\varepsilon + h$ of the interface center. The end of the extension domain is located at a distance $\varepsilon + 3h$ of the interface center to reduce diffusion impacts on the extended value. (a) Sharp interface configuration. (b) Smooth interface configuration.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** *Cont.*</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 5. The interface test case—an initial radial scalar field is extended from both sides of a thick interface (the red zone) following the normal vector field (the black arrows). Results projected on the horizontal and vertical green dotted line are shown in (b,c). (a) Schematic description of the Interface Case. (b) Source and extended field values on the horizontal axis y = 0.25. (c) Theoretical and extended field values on the centered vertical axis x = 0.</figcaption></figure>

The field to project from the left side of the interface is $q_{0L} = 5000y(1+x)$, and the associated domain $\Omega_L$ is located between the axes $x = -0.07$ and $x = 0.05$. The field to project from the right side of the interface is $q_{0R} = 10000y(1+x)$, and the associated domain $\Omega_R$ is located between the axes $x = -0.05$ and $x = 0.07$.

$\vec{n}_L$ and $\vec{n}_R$ are taken normal to the interface, but of opposite directions.

The case is implemented on an unstructured mesh of triangles of characteristic size $h = 10^{-2}$ with $\lambda = 10^{-3}$. Results projected on the central vertical axis and on the horizontal axis $y - 0.25$ are plotted on Figure 5b,c.

The extension is well computed, and errors are in the order of magnitude of the mesh precision.

### 3.4. Precision of the Method

As the extended heat flux is located one mesh size away $h$ from the boundary of the interface, an error is made, which is proportional to $h$ (see Figure 6a). Noting $s$ the position of the boundary interface, and taking the distance function $\alpha$, a linear approximation on $q$ around $s$ reads:

$$q(s+h) = q(s) + \frac{\partial q}{\partial \alpha}(s)h + o(h) \quad (15)$$
Extending $q$ from the position $s + h$ instead of $s$ thus leads to an error $\frac{\partial q}{\partial \alpha}(s)$.

Furthermore, unstructured mesh have been used in this study, and the evaluation of points $\varepsilon + h$ away from the interface used as reference points for the extension method leads at most to a doubled “safety margin”, as shown on Figure 6b.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Description of the “safety margin” taken for the extension of the heat flux that avoids errors due to the continuity of the temperature gradient. Heat flux values are taken from points away of at least $\varepsilon + h$ from the interface. As the mesh is not structured, the effective distance varies from one point to another up to $\varepsilon + 2h$. (a) Graphical illustration of the safety margin. (b) Schematic description of the selection of reference points of the extension method.</figcaption></figure>

Considering only conductive heat fluxes $q = \vec{q} \cdot \frac{\vec{\nabla}\alpha}{|\vec{\nabla}\alpha|} = -k \frac{\partial T}{\partial \alpha}$, this leads to a maximum error on $|\vec{m}|$:

$$ \Delta|\vec{m}| \propto \frac{1}{\mathcal{L}} \sum_{L,V} k \left| \frac{\partial^2 T}{\partial \alpha^2} \right| 2h \quad (16) $$

## 4. Two-Phase Thermal Solver

The presented mass transfer modeling is associated with a thermal equation and a Level Set convection equation. This solver is made to solve diphasic and “static” (no fluid velocity) Thermal systems with phase change. Thus, only divergence free cases are tackled for the moment, meaning equal density between phases.

### 4.1. Thermal Equations

As no movements are considered, a purely diffusive system with phase change is solved. According to the Continuous Surface Force approach, the integration in the energy conservation equation of the source term due to phase change from Equation (6) reads:

$$ \rho c_p \frac{\partial T}{\partial t} = \vec{\nabla} \cdot (k \vec{\nabla} T) - \mathcal{L} (\vec{m} \cdot \vec{\nabla} \alpha) \delta_{\alpha} \quad (17) $$

where $\rho$ is the density and $c_p$ the specific heat capacity.
There is an equivalence between the condition $T = T_{sat}$ at the interface and the formulation of the source term $-\mathcal{L}(\vec{m} \cdot \vec{\nabla}\alpha)\delta_{\alpha}$ in addition to the computation of the mass transfer thanks to the Heat Flux Jump Computation $\vec{m} = \left[ \frac{k\vec{\nabla}T \cdot \frac{\vec{\nabla}\alpha}{|\vec{\nabla}\alpha|}}{|\vec{\nabla}\alpha|} \right]$; formulated this way, the system theoretically maintains the interface temperature at $T_{sat}$.

To be consistent regarding the mixing laws of every parameter, a choice is made to consider that $H_{\alpha}$ is a volume ratio. This entails the mixing law of $\rho c_p$. For the mixing laws of $k$, a geometric law is used following the recommendations of former studies [36,37] that shew the importance of the mixing law to properly tackle the flux continuity:

$$ \rho c_p = H_{\alpha}(\rho c_p)_V + (1 - H_{\alpha})(\rho c_p)_L \quad (18) $$

$$ \frac{1}{k} = \frac{H_{\alpha}}{k_V} + \frac{1 - H_{\alpha}}{k_L} \quad (19) $$

As the phases are considered to have the same density $\rho$, the velocity field is zero. The mass conservation then entails the following Level Set convection Equation ($\alpha$ being positive inside the vapor phase):

$$ \frac{\partial \alpha}{\partial t} - \frac{\vec{m}}{\rho} \cdot \vec{\nabla} \alpha = 0 \quad (20) $$

## 4.2. Stabilization

The proposed system is theoretically consistent, as the formulation tends toward the sharp interface formulation for small $\epsilon$. However, the conservation of $T_{sat}$ is unfortunately not ensured as such. Among other reasons, overshooting effects are observed. Moreover the conductivity jump leads to stronger diffusion effects on the liquid side that tends to cool down the interface.

To guaranty the efficiency of the method, a penalty term of the form $A\delta_{\alpha}(T - T_{sat})$ is added, $A$ being a sufficiently large parameter (taken as $1 \times 10^6 \text{ W K}^{-1} \text{ m}^{-2}$, the order of magnitude of $\mathcal{L} \times 1 \text{ K}$). This formulation guarantees the conservation of $T_{sat}$ on the interface.

$\vec{m}$ being computed thanks to the heat flux jump, the final system reads:

$$ \rho c_p \frac{\partial T}{\partial t} + A(T - T_{sat})\delta_{\alpha} = \vec{\nabla} \cdot (k \vec{\nabla} T) - \mathcal{L}(\vec{m} \cdot \vec{\nabla} \alpha)\delta_{\alpha} \quad (21) $$

$$ \frac{\partial \alpha}{\partial t} - \frac{\vec{m}}{\rho} \cdot \vec{\nabla} \alpha = 0 \quad (22) $$

$$ \text{with: } \mathcal{L}\vec{m} = \left[ \frac{k\vec{\nabla}T \cdot \frac{\vec{\nabla}\alpha}{|\vec{\nabla}\alpha|}}{|\vec{\nabla}\alpha|} \right] \quad (23) $$

As a reminder, the advection velocity of Equation (22) represents the displacement of the interface due to the phase change. Furthermore, Equation (23) represents the energy conservation at the interface.

## 4.3. Numerical Method and Mesh

This formulation is composed of convection diffusion equations. These equations are solved using a Finite Element method, and are stabilized thanks to a SUPG-SCPG scheme presented in [38] that presents the solver used in this work. Basically, the Finite Element method consists of discretizing space into small, simple elements and choosing an approximation for the solution within each element using, for example, linear functions. Stabilization schemes such as SUPG-SCPG are necessary to ensure stability (for example in advection-dominated problems) and compatibility of the chosen approximations for each field (velocity, pressure, temperature ...). The anisotropic mesh at the interface is refined
(to improve accuracy of the solution, see [22]) using the gradient of the hyperbolic tangent of the Level Set. The physical properties of each domain are set with appropriate mixing laws. Each mesh is adapted under the constraint of a fixed, case-dependent number of edges. Representative interface regions at the end of the anisotropic adaptation process are shown in Figure 7. In practice, the adaptation process consists of adding nodes locally in the vicinity of the interface. The rest of the elements keep the same background size, that increases with the distance via three successive refinement steps to accurately capture the near-wake region.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Temperature field and mesh of the Case 3 simulation for $h = 2 \times 10^{-6}$ m at time 0.1 s. The mesh is refined around the interface where temperature gradients and physical properties variations are important. Only the left side of the computational domain is plotted.</figcaption></figure>

# 5. The Isochoric Stefan Problem

As a first simple benchmark of this model, the isochoric 1D Stefan problem is solved in a 2D domain and compared with the analytical solution. This benchmark has been first used while working with solid/fluid phase changes (see for exemple [39]). It has been recently used as a basis by Chaurasiya et al. to study such systems with time dependent conductivity and convection [40–43]. It has then been extended to configurations of liquid/gas interactions (see for exemple [19,44]).

Both phases are considered to have the same density $\rho$, though different $\rho c_p$ of a liquid/vapor mixture. Two configurations are studied: the saturated Stefan problem and the subcooled Stefan problem.

The saturated problem has been tackled many times before (see [14,16,21,22,44–46]). However the subcooled Stefan problem is rarely studied to the authors' knowledge.

## 5.1. Reminders of the Problem

The considered Stefan problem is a semi infinite domain (indexed by the $x \in [0, +\infty[$ coordinate) filled with liquid at initial temperature $T_\infty$, and touching a wall ($x = 0$) at temperature $T_w > T_{sat}$ (see Figure 8). At $t > 0$, the wall warms up the liquid and vaporization occurs, creating a moving interface positioned by $s(t)$ that goes away from the wall. The boundary is set to $T_{sat}$.
<figure><img src="image_13.png" /><figcaption>Figure 8. Schematic description of the Stefan problem test case.</figcaption></figure>

The governing equations of this problem read:

$$ \left\{ \begin{array}{ll} (\rho c_p)_V \frac{\partial T}{\partial t} = k_V \frac{\partial^2 T}{\partial x^2} & \text{for } x \in [0, s(t)] \\ (\rho c_p)_L \frac{\partial T}{\partial t} = k_L \frac{\partial^2 T}{\partial x^2} & \text{for } x \in [s(t), +\infty[ \end{array} \right. \quad (24) $$

$$ \left\{ \begin{array}{ll} (\rho c_p)_L \frac{\partial T}{\partial t} = k_L \frac{\partial^2 T}{\partial x^2} & \text{for } x \in [s(t), +\infty[ \\ (\rho c_p)_V \frac{\partial T}{\partial t} = k_V \frac{\partial^2 T}{\partial x^2} & \text{for } x \in [0, s(t)] \end{array} \right. \quad (25) $$

With the boundary conditions:

$$ T(x, t = 0) = T_{\infty} \quad \text{for } x > 0 \quad (26) $$

$$ T(x = 0, t) = T_w \quad \text{for } t > 0 \quad (27) $$

$$ T(x = +\infty, t) = T_{\infty} \quad \text{for } t > 0 \quad (28) $$

$$ T(x = s(t), t) = T_{\text{sat}} \quad \text{for } t > 0 \quad (29) $$

Moreover, the energy jump condition at the interface remaining at $T_{\text{sat}}$ and that governs the behavior of s reads:

$$ \mathcal{L}|\vec{m}| = \rho_V \mathcal{L} \frac{ds}{dt} = -k_V \frac{\partial T}{\partial x} \Big|_{x=s^{-}} + k_L \frac{\partial T}{\partial x} \Big|_{x=s^{+}} \quad (30) $$

The analytical solution classically reads:

$$ T(x, t) = \begin{cases} T_w + \frac{T_{\text{sat}} - T_w}{\text{erf}(\chi)} \text{erf}\left(\chi \frac{x}{s(t)}\right) & \text{for } x \in [0, s(t)] \\ T_{\infty} + \frac{T_{\text{sat}} - T_{\infty}}{\text{erfc}\left(\chi \sqrt{\frac{D_V}{D_L}}\right)} \text{erfc}\left(\chi \frac{x}{s(t)} \sqrt{\frac{D_V}{D_L}}\right) & \text{for } x \in [s(t), +\infty[ \end{cases} \quad (31) $$

where $D_L = \frac{k_L}{(\rho c_p)_L}$ is the liquid diffusivity and $D_V = \frac{k_V}{(\rho c_p)_V}$ is the vapour diffusivity. The evolution of s reads:

$$ s(t) = 2\chi\sqrt{D_V t} \quad (32) $$

$\chi$ is evaluated thanks to the resolution of the following equation:

$$ \rho_V \mathcal{L} \chi \sqrt{D_V} + \frac{k_V (T_{\text{sat}} - T_w) e^{-\chi^2}}{\sqrt{\pi D_V} \text{erf}(\chi)} + \frac{k_L (T_{\text{sat}} - T_{\infty}) e^{-\chi^2 \frac{D_V}{D_L}}}{\sqrt{\pi D_L} \text{erfc}\left(\chi \sqrt{\frac{D_V}{D_L}}\right)} = 0 \quad (33) $$
The temperature gradient reads:

$$ \frac{\partial T}{\partial x}(x,t) = \begin{cases} \frac{1}{\sqrt{\pi D_V t}} \frac{T_{sat} - T_w}{\text{erf}(\chi)} e^{-\left(\chi \frac{x}{s(t)}\right)^2} & \text{for } x \in [0, s(t)] \\ -\frac{1}{\sqrt{\pi D_L t}} \frac{T_{sat} - T_\infty}{\text{erfc}\left(\chi \sqrt{\frac{D_V}{D_L}}\right)} e^{-\left(\chi \frac{x}{s(t)} \sqrt{\frac{D_V}{D_L}}\right)^2} & \text{for } x \in [s(t), +\infty[ \end{cases} \quad (34) $$

One can notice that the temperature gradient (and its derivative) at the interface tends toward infinity as time gets close to zero. This entails an extra complexity for the computation of the Heat Flux Jump at the beginning of the simulation: it is interesting as it allows to test the performances of the Heat Flux Jump Computation feature and the error estimation (16).

All physical properties are those of a water and vapor mixture and summarized in Table 1. Yet a divergence free test case is considered and the density of water is taken as the density of vapor (though $\rho c_p$ values are not changing).

<table><caption>**Table 1.** Physical properties of the considered fluids for the isochoric Stefan problem.</caption><thead><tr><th></th><th>Density ρ<br>kg m<sup>-3</sup></th><th>Volume Heat Capacity ρc<sub>p</sub><br>J m<sup>-3</sup> K<sup>-1</sup></th><th>Thermal Conductivity k<br>W m<sup>-1</sup> K<sup>-1</sup></th><th>Latent Heat L<br>J kg<sup>-1</sup></th></tr></thead><tbody><tr><td>Vapor</td><td rowspan="2">5.97 × 10<sup>-1</sup></td><td>1.12 × 10<sup>3</sup></td><td>2.48 × 10<sup>-2</sup></td><td rowspan="2">2.26 × 10<sup>6</sup></td></tr><tr><td>Liquid</td><td>4.40 × 10<sup>6</sup></td><td>6.79 × 10<sup>-1</sup></td></tr></tbody></table>

## 5.2. Studied Cases

### 5.2.1. Case 1: Without Subcooling

The first case solved is the classical Stefan problem with the liquid being at saturation temperature. The wall is 10 K above the saturation temperature, and $\chi = 0.067$.

The test case is computed for an unstructured 2D mesh of dimension $1 \times 10^{-4}$ by $5 \times 10^{-4}$ m$^2$ with 5 mesh sizes $h$ (in m): $1 \times 10^{-6}$, $2 \times 10^{-6}$, $4 \times 10^{-6}$, $7 \times 10^{-6}$ and $1 \times 10^{-5}$. The time step is set to $1 \times 10^{-5}$ s. The interface thickness is set to $6h$ with an extension zone of $8h$.

The simulation starts with an initial interface position of $2 \times 10^{-5}$ m (identified as the border of the thick interface on the vapor side). The initial temperature profile is computed according to its analytical solution.

### 5.2.2. Case 2: With a Small Subcooling

The second case is similar to the first one, except that the water phase is not at saturation temperature anymore. An initial temperature for the water of $T_{sat} - 1$ K is taken, and $\chi = 0.021$. The domain, mesh sizes, time step and other parameters are the same. The initial position of the interface is also $2 \times 10^{-5}$ m.

### 5.2.3. Case 3: With a High Subcooling

The robustness of the method is assessed in computing a more extreme Stefan Case for higher temperature: the wall is 900 K above the saturation temperature, the liquid is 10 K below, and $\chi = 0.101$. Furthermore, the remeshing methods is applied to test the good behavior of the combined methods.

The test case is computed for an unstructured non uniform 2D mesh of dimension $2 \times 10^{-4}$ by $2 \times 10^{-3}$ m$^2$ with 30,000 elements whose mesh size $h$ at the interface varies among 5 values (in m): $5 \times 10^{-7}$, $1 \times 10^{-6}$, $2 \times 10^{-6}$, $3 \times 10^{-6}$ and $5 \times 10^{-6}$. The time step is still $1 \times 10^{-5}$ s. The initial position of the interface is this time $1 \times 10^{-4}$ m. The interface thickness is set to $6h$ with an extension zone of $8h$.
5.3. Results

An example of temperature field values and mesh for the Case 3 is shown in Figure 7.

Results of interface position and mass transfer rate values versus time are plotted in Figures 9a and 10a for the Case 1, in Figures 9b and 10b for the Case 2, and in Figures 9c and 10c for the Case 3.

Results of Case 1 are satisfactory, as errors on the position and mass transfer rate are small (less than 1%) whatever the mesh size.

The influence of the mesh size appears for the second case, mostly at the first time steps of the simulation where the heat flux jump computation is sensitive. This is explained by a higher heat flux gradient on the liquid size that raises the error on the mass transfer rate described by (16). Otherwise the case remains correctly described by the simulations.

<figure><img src="image_14.png" /><figcaption>Figure 9. Analytic and simulated positions of the interface for different mesh sizes *h*. (a) Case 1. (b) Case 2. (c) Case 3.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 10. *Cont.*</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 10. Analytic and simulated mass transfer rates for different mesh sizes h. (a) Case 1. (b) Case 2. (c) Case 3.</figcaption></figure>

Results for the Case 3 stress even more this observation, as even smaller mesh sizes reveal non negligible errors due to the important thermal constraint at play at the first time steps of the simulation. As soon as the temperature gradient is lower, errors on the mass transfer rate quickly decrease.

5.4. Precision and Accuracy Order

A first estimation of the HFJC method accuracy given by Equation (16) can be confirmed by plotting the errors of the three simulations on a $(h, \Delta|\vec{m}|)$ graph. This is shown on Figure 11.

The lines represent the error estimation, thanks to the computation of $\frac{\partial^2 T}{\partial x^2}(s, t = t_0)$ given by the analytical solution. The error estimation correctly fits with the simulations results. This means that these errors are due to a too coarse grid, and can be controlled and reduced with the mesh size. An optimized approach to properly solve sub-cooled Stefan problems would be to control the mesh refinement according with the time.

To further study the accuracy of the method, a convergence analysis with the mesh size h has been conducted. The error is defined as:

$$err = \frac{1}{N} \sum_{i \in \{1, N\}} |s_{simu}^i - s_{an}^i| \qquad (35)$$

where N is the number of sample steps among each iteration step, taken here as 20, and $s_{simu}^i$ and $s_{an}^i$ are the position of the simulated and analytic interface at sample step i.
Results for the three cases are plotted on Figure 12a-c.

<figure><img src="image_17.png" /><figcaption>**Figure 11.** Relative error of the initial mass transfer rate value versus interface mesh size. Estimations of this error given by Equation (16) are plotted in dashed lines</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 12.** Space-time convergence analysis for the isochoric Stefan Benchmark. (a) Case 1. (b) Case 2. (c) Case 3.</figcaption></figure>

For all cases, the convergence order is around 1 as predicted by (16). The accuracy of the HFJC method is controlled by the value of $\varepsilon$ that should be small enough to properly compute heat fluxes. If this condition is validated, and providing that the mesh is properly refined and validates the Fourier condition far from the interface, the presented method can be applied to more complex cases.
# 6. Multiphase Pseudo-Compressible Framework

To tackle the difference of density between the two phases, the consideration of compression or dilatation shall be implemented. In this framework, thermics drives the phase change, meaning a strong constrain on the mechanics. Continuing with a CSF approach, volume changes are implemented as a source term at the interface. This means the need to implement a Navier-Stokes solver that takes into account a source term on the mass conservation equation.

## 6.1. Pseudo-Compressible Navier-Stokes Solver

The one phase pseudo-compressible Navier-Stokes solver resolves a modified version of the incompressible Navier-Stokes equations:

$$ \left\{ \begin{array}{l} \vec{\nabla} \cdot \vec{u} = \varphi \\ \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla}) \vec{u} \right) = -\vec{\nabla} P + \vec{\nabla} \cdot 2\eta \dot{\epsilon}(\vec{u}) + \rho \vec{g} \end{array} \right. \qquad (36) $$

where $\varphi$ is an imposed volume term, $P$ is the pressure, $\dot{\epsilon}$ is the strain rate tensor and $\vec{g}$ the gravity.

## 6.2. Numerical Method

This Navier-Stokes system is solved with a Finite Element method and stabilized with the Variational MultiScale (VMS) method presented in previous communication of the laboratory who provided this work [34,38,47]. However, the stabilisation scheme has to be modified due to the addition of the source term in the velocity divergence equation. The new pressure residual (formerly $R_p = -\vec{\nabla} \cdot \vec{u}$) now reads:

$$ R_p = \varphi - \vec{\nabla} \cdot \vec{u} \qquad (38) $$

This slightly changes the formulation of the resolution matrix. However, the source term does not change significantly the solver accuracy and stability. The main aspects of the method, as well as the value of the stabilization parameters presented by Hachem et al. [38] are not changed.

## 6.3. A 2D Benchmark—The Pseudo-Compressible Square

We consider a first benchmark inspired by [37], whose first purpose was to study the behavior of the incompressible Navier-Stokes solver.

The target is to reproduce the following steady analytical solution on a $(0, 1) \times (0, 1)$ square:

$$ \left\{ \begin{array}{l} \vec{u}_x = -5xy^4 + Ax^2 \\ \vec{u}_y = -0.5 + y^5 \\ p = 0.5(y^5 - y^{10}) + 5\nu y^4 \\ f_x = 5xy^8 + 10xy^3 + 60\nu xy^2 - 2\nu A - 15Ax^2y^4 + 2A^2x^3 \\ f_y = 0 \\ \varphi = 2Ax \end{array} \right. $$

To do so, the values of $f_x$, $f_y$, $\varphi$ and the boundary solutions are set, as described on Figure 13. The pressure boundary conditions are set free.

Two convergence studies are conducted with two values of $A$: 0.01 and 0.1. Five mesh sizes are considered: $2 \times 10^{-3}$, $5 \times 10^{-3}$, $1 \times 10^{-2}$, $2 \times 10^{-2}$ and $5 \times 10^{-2}$.

$\rho$ is set to 1, and $\eta$ to 0.001. The simulation starts with a null velocity, and runs until a time 20. The time step is set to 0.1. Results at time 20 are compared to the analytical solution. An example of the velocity and pressure fields are plotted in Figure 14.
<figure><img src="image_19.png" /><figcaption>Figure 13. Dirichlet velocity boundary conditions for the Square Benchmark. The other velocity components are set with homogeneous Neumann boundary conditions.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 14. Final velocity and pressure fields for A = 0.1 and h = 2 × 10⁻³.</figcaption></figure>

The error criteria based on the L2 norm is used:

$$err_x = \sqrt{\frac{\int_{\Omega}(x_{simu} - x_{an})^2 dS}{\int_{\Omega}(x_{an})^2 dS}} \quad (39)$$

where x stands for the velocity vector or the pressure. 20 is long enough to ensure the steady state, when differences between two time steps are in the order of magnitude of the machine precision (see Figure 15).

<figure><img src="image_21.png" /><figcaption>Figure 15. L2 norm error evolution for A = 0.1 and h = 5 × 10⁻³. The final time considered it high enough to guaranty the steady state regime on the final iteration.</figcaption></figure>

Results of the convergence study for both variables are plotted on Figure 16 for A = 0.1 and Figure 17 for A = 0.01.

The order of convergence for the case A = 0.01 is of the same order as the one of the incompressible formulation computed by Hachem [37], proving the good behavior of the
solver for this case. For $A = 0.1$, the order of convergence is reduced (close to 1), showing a limit of the range of applicability of the formulation.

<figure><img src="image_22.png" /><figcaption>**Figure 16.** Space convergence analysis for the Square Benchmark with $A = 0.01$. (a) Velocity L2 norm error. (b) Pressure L2 norm error.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 17.** Space convergence analysis for the Square Benchmark with A = 0.1. (a) Velocity L2 norm error. (b) Pressure L2 norm error.</figcaption></figure>

6.4. *Combination with the Level Set method*

Now that the pseudo-compressible Navier-Stokes solver has been validated, it is used
to account for a phase change process between fluids of different densities within the CSF
approach. Considering back the surface mass transfer vector $\vec{m}$, the conservation of mass
in both fluids considering a Dirac distribution of mass transfer reads:
$$ egin{cases} 
\frac{\partial H_{\alpha} \rho_{V}}{\partial t} + \vec{\nabla} \cdot (H_{\alpha} \rho_{V} \vec{u}) = (\vec{m} \cdot \vec{\nabla} \alpha) \delta_{\alpha} \\ 
\frac{\partial (1 - H_{\alpha}) \rho_{L}}{\partial t} + \vec{\nabla} \cdot ((1 - H_{\alpha}) \rho_{L} \vec{u}) = -(\vec{m} \cdot \vec{\nabla} \alpha) \delta_{\alpha} 
\end{cases} 
\quad (40) 
$$

This can be rewritten as a velocity divergence equation combined with a Level Set convection equation form.

The momentum conservation equation is enriched with a Surface Tension term in the way of [47]. It is introduced in an implicit way as a volume source term $\vec{\gamma}_{\alpha}$:

$$ \vec{\gamma}_{\alpha} = \gamma_{0} \left[ \vec{\nabla} \cdot \left( \frac{\vec{\nabla} \alpha}{|\vec{\nabla} \alpha|} \right) \right] \delta_{\alpha} \vec{\nabla} \alpha \quad (42) $$

This leads to the following system:

$$ \left\{ \begin{aligned} \vec{\nabla} \cdot \vec{u} &= \left( \frac{1}{\rho_V} - \frac{1}{\rho_L} \right) (\vec{m} \cdot \vec{\nabla} \alpha) \delta_{\alpha} \\ \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla}) \vec{u} \right) &= -\vec{\nabla} P + \vec{\nabla} \cdot 2\eta \dot{\epsilon}(\vec{u}) + \vec{\gamma}_{\alpha} + \rho \vec{g} \end{aligned} \right. \quad (43) $$

$$ \begin{aligned} \rho \left( \frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \vec{\nabla}) \vec{u} \right) &= -\vec{\nabla} P + \vec{\nabla} \cdot 2\eta \dot{\epsilon}(\vec{u}) + \vec{\gamma}_{\alpha} + \rho \vec{g} \quad (44) \\ \frac{\partial \alpha}{\partial t} + \left( \vec{u} - \frac{\rho}{\rho_V \rho_L} \vec{m} \right) \cdot \vec{\nabla} \alpha &= 0 \quad (45) \end{aligned} $$

$$ \frac{\partial \alpha}{\partial t} + \left( \vec{u} - \frac{\rho}{\rho_V \rho_L} \vec{m} \right) \cdot \vec{\nabla} \alpha = 0 \quad (45) $$

The properties $\rho$ and $\eta$ are computed thanks to the following mixing laws:

$$ \rho = H_{\alpha} \rho_{V} + (1 - H_{\alpha}) \rho_{L} \quad (46) $$

$$ \eta = H_{\alpha} \eta_{V} + (1 - H_{\alpha}) \eta_{L} \quad (47) $$

This system well describes a phase change process with an arbitrary distribution $\vec{m}$. Special attention should be paid on the velocity convecting the Level Set. The velocity field $\vec{u}$ as well as the correction factor are not constant perpendicular to the interface. This means a tendency of this formulation to deform the Level Set close to the 0 isovalue. The presence of $\vec{\nabla} \alpha$ close to $\delta_{\alpha}$ allows small deviations from the distance property. However, this does not stand above a certain point. To reduce this deformation impact, the value $\rho(\alpha = 0)$ can be taken instead of $\rho$ for the correction term, without disturbing the mass conservation property. What is really important is indeed the proper convection of the 0 isovalue. To consider this also avoid the spatial error made on $\rho$. Otherwise, the Level Set function shall be reinitialised regularly. The reinitilization technique is the one described by Shakoor et al. [48]. This is a Direct Reinitialization with Trees (DRT) that uses a geometrical approach to redefine the distance function from the 0 isovalue of the Level Set function. It is based on k-dimensional trees applied for Nearest Neighbor Searches to reduce the computational time of the reinitialization procedure. Further details are available in [48] that presents the algorithm used in this work.

This formulation is tested on two different benchmarks.

## 6.5. A 1D Benchmark—Moving Interface

The first benchmark is the convection of a straight line, meaning no theoretical influence of the surface tension term. The case is solved in 2D, in a 1 m × 1 m square domain. The vertical interface separates a vapor phase of small density on the left part from a liquid phase of high density on the right part. A mass transfer $\dot{m}$ is imposed from the vapor to the liquid phases. The fluids can only escape on the right side, leading to an interface displacement from left to right along with the liquid being pushed away due to the change
of density. No slip boundaries conditions are fixed on the left side, and free slip conditions are fixed on the top and bottom sides. Pressure conditions are set free. This is summed up in Figure 18.

<figure><img src="image_24.png" /><figcaption>**Figure 18.** Description of the Moving Interface Benchmark. Dirichlet conditions are described in this figure. Boundary conditions for the velocity components and pressure that are not mentioned are set free.</figcaption></figure>

A constant interface velocity of $u_I = 2 \times 10^{-3} \text{ m s}^{-1}$ is imposed. Different density ratios $\rho_L/\rho_V$ are considered ranging from 2 to 2000 by changing $\rho_V$. $\dot{m}$ is set accordingly to have the desired $u_I$. The interface starts at a position 0.4 from the left side. The benchmark is conducted for the fluid properties gathered in Table 2.

<table><caption>**Table 2.** Physical properties of the considered fluids of the Moving Interface benchmark.</caption><thead><tr><th></th><th>Density ρ<br>(kg m<sup>-3</sup>)</th><th>Viscosity η<br>(Pa s)</th><th>Surface Tension γ<sub>0</sub><br>(J m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Vapor</td><td>5 × 10<sup>-1</sup>/5/5 × 10<sup>1</sup>/5 × 10<sup>2</sup></td><td>1.2 × 10<sup>-5</sup></td><td rowspan="2">6 × 10<sup>-2</sup></td></tr><tr><td>Liquid</td><td>10<sup>3</sup></td><td>2.8 × 10<sup>-4</sup></td></tr></tbody></table>

No reinitialization was conducted for this benchmark, leading to a deformed Level Set. This enhances the capability of the formulation to handled non reinitialised Level Set on simple configurations.

The case was computed with four mesh sizes (in m): 5 × 10⁻³, 1 × 10⁻², 2 × 10⁻² and 5 × 10⁻². Four time steps were set accordingly to maintain a constant CFL value of 0.2 related to the interface velocity $u_I = \dot{m}/\rho_V$.

Two convergence studies were conducted. In the first one, the parameter $\epsilon$ was set constant to 0.1 m. In the second one, $\epsilon$ was set consistently with the mesh size $h$: $\epsilon = 3h$.

The considered error is also computed with (35) but with $N = 10$ sample points.

### 6.5.1. Constant $\epsilon$

Positions for every mesh sizes for the case $r = 2000$ is plotted on Figure 19, and the convergence study for every ratios is plotted in Figure 20.
<figure><img src="image_25.png" /><figcaption>Figure 19. Space-time convergence analysis for the Moving Interface benchmark for ε = 0.1 m.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 20. Interface position in the case r = 2000 for ε = 0.1 m for the Moving Interface Benchmark.</figcaption></figure>

Displacement of the interface is properly computed according to the analytical solution, validating the combination of the two solvers. This is true whatever the ratio of density. The convergence is well observed, with an order of 2 whatever the density ratio. This is coherent with the results of the pseudo-compressible square benchmark with A = 0.01.

Another point of interest is to verify if the mass conservation is properly respected. The good capturing of the interface displacement tells us that the mass of vapor is consistent with the analytical value. However, a liquid outflux exists on the right boundary, meaning a potential source of error. We therefore compare the total mass of the 2D system M:

$$M = H \int_{0}^{s} \rho_{V} dx + H \int_{s}^{L} \rho_{L} dx \quad (48)$$

where L and H are the width and the height of the domain.

The analytical mass $M_{an}$ is given by the analytical value of s. It is compared with the simulated mass $M_{simu}$ computed by (48), and with the integrated mass $M_{out}$ from the liquid outflux:

$$M_{out} = M(t = t_0) - \int_{0}^{t} \int_{0}^{H} \vec{u}_{L} \cdot \vec{e}_{x} dy dt \quad (49)$$
where $\vec{e}_x$ is the horizontal vector pointing to the right.

It is also compared with the mass $M_{div}$ computed by integrating the divergence of the velocity field:

$$M_{div} = M(t = t_0) - \int_{0}^{t} \int_{0}^{H} \int_{0}^{L} \vec{\nabla} \cdot \vec{u} \, dx \, dy \, dt \quad (50)$$

We consider the deviations of $M_{simu}$, $M_{out}$ and $M_{div}$ from $M_{an}$ thanks to the relative difference $|M_* - M_{an}|/M_{an}$. These three quantities are plotted in Figure 21 for the configuration $h_{min} = 5 \times 10^{-3}$ m and $\Delta t = 0.5$ s. As the interface is properly computed, the simulated mass is naturally well computed. What is important is to notice that the masses estimated from the outflux and from the velocity divergence is consistent with the absolute mass, as the difference remains below 1%. This validates the efficiency of the system to properly deal with mass conservation.

<figure><img src="image_27.png" /><figcaption>Figure 21. Comparison of the analytical and simulated mass of the moving interface benchmark for the configuration $h_{min} = 5 \times 10^{-3}$ m and $\Delta t = 0.5$ s. The estimations from the outflux and the velocity divergence are consistent with the instant mass value.</figcaption></figure>

### 6.5.2. Varying $\varepsilon$

The idea of the formulation is to reduce the interface size as much as possible, as the goal is to model a quasi sharp interface. The physical mixing length of a vapor liquid interface is of the order of magnitude of an Å, though this scale is rarely reached in practice. $\varepsilon$ is then set according to the mesh size value. The same law as for the Stefan problem described in Section 5 is taken: $\varepsilon = 3h$.

Doing so, two problems arise. Firstly, numerical instabilities described by [49] are not handled anymore and disturb the interface. This is not surprising, as the BKZ condition $(\Delta t_{BZK} < \sqrt{(\rho_V + \rho_L)h^3/4\pi\gamma_0})$ is exceeded by 2 orders of magnitudes, and the way the surface tension is handled does not fully control these instability, according to [49]. However, this problem can be handled by reducing the time step.

Secondly, even while canceling the surface tension term or stabilising the case by increasing viscosity values by a factor 100, no proper convergence is observed up to a certain point, and a residual error of around 1% remains. For example, interface positions for a ratio $r = 2000$ with higher viscosity are plotted Figure 22.
<figure><img src="image_28.png" /><figcaption>Figure 22. Interface position in the case $r = 2000$ for $\varepsilon = 3h$.</figcaption></figure>

This can be explained by the discretization of the Dirac function $\delta_{\alpha}$. Considering $\varepsilon$ proportional to $h$ means that the $\delta_{\alpha}$ function is discretized with the same number of point whatever the mesh size: 7 discretization points for $\varepsilon = 3h$ (see Figure 23). However most important is the capability of the function $\delta_{\alpha}$ to guarantee an integral of 1 over the interface. However, this is not guaranteed with a fixed amount of sample points in an unstructured mesh. To demonstrate this, a Python code has been implemented to assess the order of magnitude of the error made on the integral computation due to the discretization in 1D.

For the sake of an example, considering a constant mesh size $h$ of 1, and still with $\varepsilon = 3h$, seven reference sample points are considered, enriched with two extremity sample points at a distance $4h$ from the center to cover the entire function:

$$ \tilde{\alpha}_i \in [0,8] = -(3+1)h + hi \quad (51) $$

A batch of 1000 groups of sample points are created thanks to the random function `rd()` from the ` numpy` library of Python, that allocates a random real between 0 and 1. Every sample point is randomly placed around a reference sample point at a maximum distance $\zeta$. This is supposed to represent the non uniformity of the mesh. The bigger $\zeta$, the more distorted the mesh is.

$$ \{\alpha_i\}_{j \in [1,1000]} = \left\{ \tilde{\alpha}_i + \zeta \frac{rd() - 0.5}{0.5} \right\} \quad (52) $$

An example of a resulting sampling is plotted in Figure 23.

For different values of $\zeta$, a thousand $\delta_{\alpha}$ functions with randomised discretization points are created. Integrals of these functions are then computed. Mean values of the error made are plotted in Figure 24.

Tests has been carried out with groups of points bigger than 1000 without changing the results. Thus, these results are assumed to represent well the discretization error conducted with an unstructured mesh. The order of magnitude of the error is of the percent. This is the same as the one observed for the interface displacement with $\varepsilon$ changing with $h$. This explains why the formulation for a varying interface thickness does not properly converge in mesh size, as a permanent residual error exists due to the sampling of the Dirac function: if the integral is not properly computed, the integrated velocity field resulting from the velocity divergence equation deviates from the analytical solution. Then, the convection velocity of the Level Set equation also deviates from the analytical solution due to the same error, leading to the error on the interface position.
<figure><img src="image_29.png" /><figcaption>**Figure 23.** Example of discretization of the Dirac function $\delta_\alpha$ in 1D. A thousand similar discretizations are conducted in the same way using randomized positions around the reference points.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 24.** Influence of the value of $\zeta$ on the mean value of the integral absolute deviation of the discretized Dirac functions of every batch.</figcaption></figure>

### 6.6. A 2D Benchmark—Growing Bubble

For this 2D benchmark, the contribution of the surface tension term is more strongly observed. Inside a square domain of 0.008 m × 0.008 m, a circular bubble of low density grows with a constant mass transfer inside a fluid of higher density, as shown in Figure 25. This fluid is pushed away outside of the domain, whose boundary conditions are of homogeneous Neumann type. Inside the bubble, a small zone is set at 0 velocity to prevent the derivation of the Bubble.

A constant mass transfer of 0.1 kg m⁻² s⁻¹ from the vapor phase towards the liquid phase is set. The initial radius of the bubble is 0.001 m. The benchmark is conducted for the following fluid properties gathered in Table 3.

This time, the Level Set was reinitialized to avoid important distortions of the Level Set. The reinitialization is conducted every 0.001 s, thus 10 times during the simulation. This is high enough to ensure that the Level Set is not too deformed, but low enough to avoid any systematic error due to the DRT method detailed in [48].
<figure><img src="image_31.png" /><figcaption>**Figure 25.** Description of the Growing Bubble Benchmark.</figcaption></figure>

<table><caption>**Table 3.** Physical properties of the considered fluids for the Growing Bubble test case.</caption><thead><tr><th></th><th>Density ρ<br>(kg m<sup>-3</sup>)</th><th>Viscosity η<br>(Pa s)</th><th>Surface Tension γ<sub>0</sub><br>(J m<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Vapor</td><td>1</td><td>1.78 × 10<sup>-5</sup></td><td rowspan="2">7 × 10<sup>-2</sup></td></tr><tr><td>Liquid</td><td>10<sup>3</sup></td><td>10<sup>-3</sup></td></tr></tbody></table>

The case was computed with four mesh sizes: (in m): 2 × 10<sup>-5</sup>, 5 × 10<sup>-5</sup>, 1 × 10<sup>-4</sup> and 2 × 10<sup>-4</sup>. Four different time steps were set accordingly to maintain a constant CFL value of 0.025 related to the interface velocity $\dot{m}/\rho V$. With these parameters, the BKZ condition is roughly respected.

Two studies were conducted: in the first one, the parameter $\epsilon$ was set constant to 2 × 10<sup>-4</sup> m. In the second one, $\epsilon$ was set consistently with the mesh size $h$: $\epsilon = 3h$.

The considered error is based on the radius $R$, among $N = 10$ sample points:

$$err = \frac{1}{N} \sum_{i \in \{1, N\}} |R_{simu}^i - R_{an}^i| \qquad (53)$$

The study with constant $\epsilon$ leads to a residual error on the radius, even with the smallest values of time step and mesh size, as shown in Figure 26.

The reason lies inside the thick nature of the interface. When this parameter is not negligible in comparison to the curvature of radius $R$, a bias in the formulation is introduced. Globally, the amount of mass is conserved, but locally, the 0 isovalue is not properly convected. This can be shown while expressing analytically the velocity for this benchmark. The Level Set is here positive in the bubble, and negative outside, with $H_\alpha$ being 1 inside the bubble, and 0 otherwise. Considering that the Level Set is properly reinitialised, writing

$\vec{u} = u_r(r)\vec{e}_r$, and integrating $\vec{\nabla} \cdot \vec{u} = \frac{1}{r} \frac{\partial r u_r}{\partial r} = \left[\frac{1}{\rho}\right] (\vec{m} \cdot \vec{\nabla} \alpha) \delta_\alpha$ between $R - \epsilon$ and $r$ leads to:

$$u_r = \left[\frac{1}{\rho}\right] |\vec{m}| \left( (1 - H_\alpha) - \frac{1}{r} \int_{-\epsilon}^{r-\bar{R}} H_\alpha \, d\alpha \right) \qquad (54)$$

Adding the velocity correction term, the 0 isovalue convection velocity reads:

$$\vec{u}_I(\alpha = 0) = -\frac{\vec{m}}{\rho V} + \left[\frac{1}{\rho}\right] \frac{\vec{m}}{R} \int_{-\epsilon}^{0} H_{\alpha} \, d\alpha \qquad (55)$$
Which is not the analytical solution $-\frac{\vec{m}}{\rho_V}$. A second order term appears, and is not negligible when $\epsilon$ is not negligible in regards to the curvature radius $R$. This shows the need of an order 2 correction term to account for this effect. If the Level Set is not properly reinitialised, an order 1 linear approximation leads to a modification of the correction term thanks to the gradient of the Level Set.

<figure><img src="image_32.png" /><figcaption>Figure 26. Impact of the second order correction term on the LS convection velocity the Growing Bubble benchmark for $\epsilon = 2 \times 10^{-4}$ m.</figcaption></figure>

The convection equation now reads:

$$ \frac{\partial \alpha}{\partial t} + \left( \vec{u} - \frac{\rho}{\rho_V \rho_L} \vec{m} - \left[ \frac{1}{\rho} \right] \vec{m} \frac{\kappa_\alpha}{|\vec{\nabla} \alpha|} \int_{-\epsilon}^{0} H_\alpha \, d\alpha \right) \cdot \vec{\nabla} \alpha = 0 \quad (56) $$

With $\kappa_\alpha$ the curvature that is computed thanks to the divergence of the normal given by the Level Set, as for the Surface Tension term. The integral can be computed analytically knowing the analytical expression of $H_\alpha$.

Doing so, the problem is solved and the error vanishes for the smallest values of time step and mesh size as shown in Figure 26.

Results of the convergence study are shown in Figure 27.

Now, a proper convergence is observed. An order of convergence close to 2 is coherent with the former results of the moving interface.

Regarding the convergence study when varying $\epsilon$, the radius for different time steps and mesh size are plotted in Figure 28.

The error is not very significant for small $h$ and $\Delta t$, but convergence is slower with mesh size due to the Dirac discretization that requires a very small mesh size. The order of magnitude of the error is again the same as the one estimated thanks to our Python code, explaining the convergence of the scheme for varying $\epsilon$. This is however an error of relative small importance for more complex cases, and that can be reduced by raising the number of discretization points to describe $\delta_\alpha$.

The mass conservation is also verified in this configuration. We follow the same approach as in Section 6.5.1, the only difference being that this configuration is 2D: the outflux has to be computed on all the boundaries of the domain. The different deviations of the mass estimations for the configuration $h_{min} = 2 \times 10^{-5}$ m and $\Delta t = 5 \times 10^{-6}$ s are plotted in Figure 29. The masses estimated from the outflux and from the velocity divergence are still consistent with the absolute mass. This validates further the efficiency of the system to properly deal with mass conservation.
<figure><img src="image_33.png" /><figcaption>**Figure 27.** Space-time convergence analysis for the Growing Bubble benchmark for $\varepsilon = 2 \times 10^{-4}$ m.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 28.** Radius for different mesh sizes and time steps for the growing bubble benchmark for $\varepsilon = 3h$.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 29.** Comparison of the analytical and simulated mass of the growing bubble benchmark for the configuration $h_{min} = 2 \times 10^{-5}$ m and $\Delta t = 5 \times 10^{-6}$ s. The estimations from the outflux and the velocity divergence are still consistent with the instant mass value.</figcaption></figure>
## 7. The Complete Phase Change Solver

All the ingredients are now ready for the complete formulation of a phase change solver. The mass transfer and change of density are taken into account thanks to the pseudo-compressible solver along with Level Set convection. The surface mass transfer vector $\vec{m}$ is computed thanks to the HFJC method described in Section 3.

Special attention should be paid on the distance function property. The HFJC relies on a proper distance function, and can run improperly if the Level Set is deformed. To ensure this, a working distance function used for the HFJC method is created by reinitializing the Level Set at each time step without modifying $\alpha$ itself.

Two benchmarks are studied to assess the validity of the complete formulation.

### 7.1. A 1D Benchmark—Compressible Stefan Test Case

The compressible Stefan test case is very similar to the incompressible configuration, except that the liquid phase is pushed back away from the wall. There is thus a convection term inside the liquid phase that changes the formulation of the energy conservation equation, and then the temperature distribution.

In the liquid phase, the new equation reads:

$$ \frac{\partial T}{\partial t} + u_L \frac{\partial T}{\partial x} = D_L \frac{\partial^2 T}{\partial x^2} \qquad (57) $$

The velocity $u_L$ is spatially constant, and determined thanks to the velocity jump at the interface:

$$ u_L = \left[ \frac{1}{\rho} \right] |\vec{m}| \qquad (58) $$

Boundary conditions (30) do not change. Interface still initialized at 0.0001 m.

The solution of this new version of Stefan benchmark reads:

$$ T(x,t) = \begin{cases} T_{\omega} + \frac{T_{sat} - T_{\omega}}{\text{erf}(\chi)} \text{erf}\left(\chi \frac{x}{s(t)}\right) & \text{for } x \in [0, s(t)] \\ T_{\infty} + \frac{T_{sat} - T_{\infty}}{\text{erfc}\left(\chi \frac{\rho_V}{\rho_L} \sqrt{\frac{D_V}{D_L}}\right)} & \\ \times \text{erfc}\left(\chi \frac{x - \left(1 - \frac{\rho_V}{\rho_L}\right) s(t)}{s(t)} \sqrt{\frac{D_V}{D_L}}\right) & \text{for } x \in [s(t), +\infty[ \end{cases} \qquad (59) $$

The interface displacement is still proportional to the square root of the time:

$$ s(t) = 2\chi\sqrt{D_V t} \qquad (60) $$

With $\chi$ given from the energy conservation at the interface that leads to:

$$ \rho_V \mathcal{L} \chi \sqrt{D_V} + \frac{k_V (T_{sat} - T_{\omega}) e^{-\chi^2}}{\sqrt{\pi D_V} \text{erf}(\chi)} + \frac{k_L (T_{sat} - T_{\infty}) e^{-\chi^2} \rho_V^2 D_V}{\sqrt{\pi D_L} \text{erfc}\left(\chi \frac{\rho_V}{\rho_L} \sqrt{\frac{D_V}{D_L}}\right)} = 0 \qquad (61) $$

Only one configuration is tested: $T_{\omega} = T_{sat} + 10$ K, $T_{\infty} = T_{sat} - 1$ K. Then $\chi = 2.52 \times 10^{-2}$. The physical properties are given in Table 4.
<table><caption>**Table 4.** Physical properties of the considered fluids for the Compressible Stefan test case.</caption><thead><tr><th></th><th>Density ρ<br>(kg m<sup>-3</sup>)</th><th>Specific Heat Capacity c<sub>p</sub><br>(J kg<sup>-1</sup> K<sup>-1</sup>)</th><th>Thermal Conductivity k<br>(W m<sup>-1</sup> K<sup>-1</sup>)</th><th>Viscosity η<br>(Pa s)</th><th>Latent Heat L<br>(J kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Vapor</td><td>5.97 × 10<sup>-1</sup></td><td>2.030 × 10<sup>3</sup></td><td>2.48 × 10<sup>-2</sup></td><td>1.20 × 10<sup>-5</sup></td><td rowspan="2">2.26 × 10<sup>6</sup></td></tr><tr><td>Liquid</td><td>9.584 × 10<sup>2</sup></td><td>4.216 × 10<sup>3</sup></td><td>6.76 × 10<sup>-1</sup></td><td>2.8 × 10<sup>-4</sup></td></tr></tbody></table>

The test case is computed for an unstructured non uniform 2D mesh of dimension $2 \times 10^{-4}$ by $2 \times 10^{-3}$ m² with a maximum of 50.000 elements whose mesh size $h$ at the interface varies among 5 values (in m): $1 \times 10^{-6}$, $2 \times 10^{-6}$, $4 \times 10^{-6}$, $7 \times 10^{-6}$ and $1 \times 10^{-5}$. Five different time steps were set accordingly to maintain a constant CFL value of 0.2 related to the initial interface velocity $\dot{m}(t = t_0)/\rho_V$. The initial position of the interface is this time $1 \times 10^{-4}$ m. The interface thickness is set to $3h$ with an extension zone of size $4h$.

Results of interface position and mass transfer rate values versus time are plotted in Figures 30 and 31.

<figure><img src="image_36.png" /><figcaption>**Figure 30.** Position of the interface for different mesh sizes and time steps for the compressible Stefan problem benchmark.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Figure 31.** Mass transfer for different mesh sizes and time steps for the compressible Stefan problem benchmark.</figcaption></figure>

The convergence is considered regarding the error on the interface position, computed with (35) with $N = 100$ sample points. The convergence is plotted in Figure 32.
<figure><img src="image_38.png" /><figcaption>**Figure 32.** Space-time convergence analysis for the compressible Stefan problem benchmark.</figcaption></figure>

The convergence is well observed, with an order 1 that is coherent with the prediction of (16). The computed errors are above the order of magnitude of the one caused by the discretization of the Dirac, explaining why it does not disturb the convergence. This may appear for smaller time steps and mesh sizes.

The mass conservation is once again verified. We strictly follow the same approach as in Section 6.5.1. The different mass deviations for the configuration $h_{min} = 2 \times 10^{-6}$ m and $\Delta t = 2 \times 10^{-5}$ s are plotted in Figure 33. The masses estimated from the outflux and from the velocity divergence are still consistent with the absolute mass, though the error is a little higher but still below 1%.

<figure><img src="image_39.png" /><figcaption>**Figure 33.** Comparison of the analytical and simulated mass of the compressible Stefan benchmark for the configuration $h_{min} = 2 \times 10^{-6}$ m and $\Delta t = 2 \times 10^{-5}$ s. The mass deviations are higher at the beginning of the simulation due to the error in mass transfer that is higher at higher heat fluxes.</figcaption></figure>

## 7.2. A 1D Benchmark—Adiabatic Problem

We propose to tackle a complementary study with adiabatic boundary conditions inspired by the work of Santiago et al. [39]. This allows to challenge differently the respect
of the energy conservation. It consists of setting an arbitrary temperature field on a domain with both liquid and vapor separated by a planar interface, and to let the system evolve to the steady state – either one of the phase vanishes, or the system temperature tends to $T_{sat}$. In our consideration, the vapor is blocked by the wall at $x = 0$, and the liquid is allowed to move with an initial right boundary at $x = L_0 = L(t_0)$, as shown in Figure 34.

<figure><img src="image_40.png" /><figcaption>Figure 34. Schematic description of the adiabatic problem.</figcaption></figure>

If the initial temperature field $T(0 < x < L, t_0)$ and the initial interface position $s_0 = s(t_0)$ are set so that the steady state is a mixture of liquid and vapor at $T_{sat}$, then the available thermal energy $\mathcal{E}$ over a section S reads:

$$ \mathcal{E} = S\rho_V c_{pV} \int_0^{s_0} (T(x, t_0) - T_{sat}) dx - S\rho_L c_{pL} \int_{s_0}^{L_0} (T_{sat} - T(x, t_0)) dx \quad (62) $$

Over the conditions that $\mathcal{E}$ is positive, meaning that there is an excess of energy regarding the steady state at $T_{sat}$ for the same proportion of liquid and vapor, all this energy is transformed to vaporise a portion of liquid:

$$ \mathcal{E} = \rho_V S \mathcal{L}(s_\infty - s_0) = \rho_V S \mathcal{L} \Delta s \quad (63) $$

with $s_\infty = s(t \to \infty)$ the final interface position.

Thus, for a given initial temperature field that fulfills these conditions, the final displacement of the liquid vapor interface is:

$$ \Delta s = \left( \rho_V c_{pV} \int_0^{s_0} (T(x, t_0) - T_{sat}) dx - \rho_L c_{pL} \int_{s_0}^{L_0} (T_{sat} - T(x, t_0)) dx \right) / (\rho_V \mathcal{L}) \quad (64) $$

This is related to a certain amount of liquid that vaporizes, leading to the liquid interface displacement:

$$ \Delta L = \left(1 - \frac{\rho_V}{\rho_L}\right) \Delta s \quad (65) $$

In our study, the initial temperature profile considered reads:

$$ T(x, t_0) = \begin{cases} (T_{sat} - T_\omega) \left(\frac{x}{s_0}\right)^2 + T_\omega & \text{for } x \in [0, s_0] \\ (T_{sat} - T_\infty) \left(\frac{L_0 - x}{L_0 - s_0}\right)^2 + T_\infty & \text{for } x \in [s_0, L_0] \end{cases} \quad (66) $$

with $T_\omega = 1100$ °C and $T_\infty = 99.9$ °C. The liquid vapor interface is set at $s_0 = 0.5$ mm from the wall, and the liquid boundary is initially at $L_0 = 1$ mm from the wall. The material properties of liquid and water are taken from Table 4.
The test case is computed for an unstructured non uniform 2D mesh of dimension
$2 \times 10^{-4}$ by $2 \times 10^{-3} \text{ m}^2$ with a maximum of 25.000 elements whose mesh size $h_{min}$ at the
interface is set to $2 \times 10^{-6} \text{ m}$. The end of the liquid interface is represented by another Level
Set to enable its movement. To guaranty adiabatic conditions, the rest of the domain ($x > L$)
is set with a reduced conductivity ($1 \times 10^{-6} \text{ W m}^{-1} \text{ K}^{-1}$). This way the temperature gradient
is approximately null at the liquid boundary. A time step of $1 \times 10^{-5} \text{ s}$ is set. The simulation
runs until a time $t_{\infty} = 5 \text{ s}$ to ensure that it reaches the steady state.

The liquid vapor interface and the liquid boundary position are plotted in Figure 35.
First of all, the steady state is consistent with the analytical value. The absolute final error
on both interface positions is around $5 \times 10^{-6}$ m, which is of the order of magnitude of the
minimum mesh size, and three orders of magnitude lower than the system lengthscale,
proving a proper respect of the energy conservation. It is also interesting to notice that the
system vaporises too much before condensing to reach the steady state, which is consistent
with the difference of diffusivity between the liquid and the vapor phases. The vapor phase
diffusivity being higher, conduction inside the vapor is stronger and the steady state is
reached faster. Thus, all the thermal potential energy of vapor is given to the interface
before the liquid recovers part of it.

<figure><img src="image_41.png" /><figcaption>**Figure 35.** Comparison of the simulated liquid vapor interfaces and liquid boundary positions compared with the steady state analytical solution of the adiabatic problem.</figcaption></figure>

**8. Conclusions and Perspectives**

In this article, an extension of an existing Finite Element solver based on the Continuous Surface Force approach and adapted to phase change systems has been developed. In comparison to other phase change approaches, it has the advantage of an easy and light numerical implementation. The main conclusions are gathered hereafter:

* This methods has been tested and validated on several analytic tests and isochoric Stefan Problem without and with subcooling.

* Then a remeshing algorithm has been added and good results are showed for the combination of these two methods on a stiffer Stefan problem.

* A Navier-Stokes solver embedding a source term in the mass conservation equation has been developed and coupled with a new formulation of the Level Set convection equation to account for a density change during phase change. This framework has been tested and validated on simple 1D and 2D benchmarks. The discretization of the interface has been shown to be a critical feature to control for a proper behavior of the method.
* Finally, these two frameworks have been merged to create a complete phase change solver. It has been tested on a compressible Stefan problem. The conclusion is that this solver is valid for stable configurations where no phase is overheated or undercooled.

* The model is based on several simplifications: the fluids are considered incompressible, only thermal energy variations are prevailing, the interface is in thermal equilibrium, and the viscous dissipation is negligible. This framework is then applicable only in those conditions. Thus, violent thermal shock waves are not modeled, which can be troublesome when dealing with boiling modeling. Furthermore, second order effects due to the compressibility of vapor are not taken into account, and the simulation can under-estimate hydrodynamical effects.

Future works are expected to deal with 3D tests and add a suitable wetting model.

**Author Contributions:** Conceptualization, C.B. and R.V.; software, C.B. and E.H.; validation, C.B. and L.C.-B.; formal analysis, C.B.; investigation, C.B. and L.C.-B.; writing—original draft preparation, C.B. and L.C.-B.; writing—review and editing, R.V.; supervision, E.H. and R.V.; project administration, E.H.; funding acquisition, E.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The author gratefully acknowledges support from the companies involved in the INFINITY industrial ANR chair: Aubert & Duval, Cevifal, CMI Group, Faurecia, Framatome, Industeel CRMC (Arcelor Mittal), Lisi Aerospace, Montupet, Safran, Science Computer and Consultants (SCC) and Transvalor.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Abbreviations**

The following abbreviations are used in this manuscript:

<table><tr><td>LS</td><td>Level Set</td></tr><tr><td>CSF</td><td>Continuous Surface Force</td></tr><tr><td>t</td><td>Time (s)</td></tr><tr><td>ū</td><td>Velocity (m·s<sup>-1</sup>)</td></tr><tr><td>P</td><td>Pressure (Pa)</td></tr><tr><td>T</td><td>Temperature (K)</td></tr><tr><td>q̇</td><td>Heat flux (W·m<sup>-2</sup>)</td></tr><tr><td>ṁ</td><td>Surface mass transfer rate vector (kg·s<sup>-1</sup>·m<sup>-2</sup>)</td></tr><tr><td>s</td><td>Interface position (m)</td></tr><tr><td>R</td><td>Radius (m)</td></tr><tr><td>M</td><td>Mass (kg)</td></tr><tr><td>dṀ</td><td>Local mass transfer rate (kg·s<sup>-1</sup>)</td></tr><tr><td>ū<sub>I</sub></td><td>Interface velocity (m·s<sup>-1</sup>)</td></tr><tr><td>dS</td><td>Infinitesimal surface (m<sup>2</sup>)</td></tr><tr><td>dV</td><td>Infinitesimal volume (m<sup>3</sup>)</td></tr><tr><td>ρ</td><td>Density (kg·m<sup>-3</sup>)</td></tr><tr><td>η</td><td>Dynamic viscosity (Pa·s)</td></tr><tr><td>c<sub>p</sub></td><td>Specific heat capacity at constant pressure (J·K<sup>-1</sup>·kg<sup>-1</sup>)</td></tr><tr><td>k</td><td>Thermal conductivity (W·m<sup>-1</sup>·K<sup>-1</sup>)</td></tr><tr><td>T<sub>sat</sub></td><td>Saturation temperature of the fluid at atmospheric pressure (K)</td></tr><tr><td>L</td><td>Specific enthalpy of vaporization at saturation temperature (J·kg<sup>-1</sup>)</td></tr><tr><td>D</td><td>Thermal diffusivity (m<sup>2</sup>·s<sup>-1</sup>)</td></tr><tr><td>γ<sub>0</sub></td><td>Surface tension coefficient (J·m<sup>-2</sup>)</td></tr><tr><td>ġ</td><td>Gravity (m·s<sup>-2</sup>)</td></tr><tr><td>ε̇</td><td>Strain rate tensor (s<sup>-1</sup>)</td></tr><tr><td>χ</td><td>Stefan constant (Ø)</td></tr></table>
L Liquid

V Vapor

ω Wall

∞ Far from the interface

n̄ Normal vector of the liquid-vapor interface

d Distance function

α Level Set function

Hα Heaviside function

δα Dirac function

ε Characteristic length of the smoothed interface

∇⃗ Gradient operator

∇· Divergence operator

λ Diffusivity of the extension solver

## References

1. Nukiyama, S. The maximum and minimum values of the heat Q from metal to boiling water under atmospher pressure. *J. Jpn. Soc. Mech. Eng.* **1934**, *37*, 367–374.

2. Zuber, N. *Hydrodynamic Aspects Of Boiling Heat Transfer (Thesis)*; Technical Report; University of California: Los Angeles, CA, USA, 1959. [CrossRef]

3. Berenson, P.J. *Transition Boiling Heat Transfer from a Horizontal Surface*; Technical Report 17; The National Science Foundation, Cambridge, MA, USA 1960.

4. Chen, J.C. Correlation for Boiling Heat Transfer to Saturated Fluids in Convective Flow. *Industrial & Engineering Chemistry Process Design and Development* **1966**, *5*, 322–329. [CrossRef]

5. Judd, R.L.; Hwang, K.S. A Comprehensive Model for Nucleate Pool Boiling Heat Transfer Including Microlayer Evaporation. *J. Heat Transfer* **1976**, *98*, 623–629. [CrossRef]

6. Cooper, M. Saturation Nucleate Pool Boiling—A Simple Correlation. In *First U.K. National Conference on Heat Transfer*; Elsevier: Amsterdam, The Netherlands, 1984; pp. 785–793. [CrossRef]

7. Ramilison, J.M.; Lienhard, J.H. Transition Boiling Heat Transfer and the Film Transition Regime. *J. Heat Transfer* **1987**, *109*, 746–752. [CrossRef]

8. Steiner, H.; Kobor, A.; Gebhard, L. A wall heat transfer model for subcooled boiling flow. *Int. J. Heat Mass Transfer* **2005**, *48*, 4161–4173. [CrossRef]

9. Meduri, P.K.; Warrier, G.R.; Dhir, V.K. Wall heat flux partitioning during subcooled forced flow film boiling of water on a vertical surface. *Int. J. Heat Mass Transfer* **2009**, *52*, 3534–3546. [CrossRef]

10. Brennen, C.E. *Fundamentals of Multiphase Flow*, reprint ed.; Cambridge University Press: Cambridge, UK, 2006.

11. Lee, W.H. A Pressure Iteration Scheme for Two-Phase Flow Modeling. In *Computational Methods for Two-Phase Flow and Particle Transport*; World Scientific: Singapore, 1980; pp. 61–82. [CrossRef]

12. Rattner, A.S.; Garimella, S. Simple Mechanistically Consistent Formulation for Volume-of-Fluid Based Computations of Condensing Flows. *J. Heat Transfer* **2014**, *136*, 071501. [CrossRef]

13. Giustini, G.; Issa, R.I. A method for simulating interfacial mass transfer on arbitrary meshes. *Phys. Fluids* **2021**, *33*, 087102. [CrossRef]

14. Son, G.; Dhir, V.K. Numerical Simulation of Film Boiling Near Critical Pressures With a Level Set Method. *J. Heat Transfer* **1998**, *120*, 183–192. [CrossRef]

15. Juric, D. Computations of Boiling Flows. *Int. J. Multiph. Flow* **1998**, *24*, 387–410. [CrossRef]

16. Welch, S.W.; Wilson, J. A Volume of Fluid Based Method for Fluid Flows with Phase Change. *J. Comput. Phys.* **2000**, *160*, 662–682. [CrossRef]

17. Kang, M.; Fedkiw, R.P.; Liu, X.D. A Boundary Condition Capturing Method for Multiphase Incompressible Flow. *J. Sci. Comput.* **2000**, *15*, 323–360. [CrossRef]

18. Tanguy, S.; Ménard, T.; Berlemont, A. A Level Set Method for vaporizing two-phase flows. *J. Comput. Phys.* **2007**, *221*, 837–853. [CrossRef]

19. Gibou, F.; Chen, L.; Nguyen, D.; Banerjee, S. A level set based sharp interface method for the multiphase incompressible Navier–Stokes equations with phase change. *J. Comput. Phys.* **2007**, *222*, 536–555. [CrossRef]

20. Esmaeeli, A.; Tryggvason, G. Computations of film boiling. Part I: numerical method. *Int. J. Heat Mass Transf.* **2004**, *47*, 5451–5461. [CrossRef]

21. Sato, Y.; Ničeno, B. A sharp-interface phase change model for a mass-conservative interface tracking method. *J. Comput. Phys.* **2013**, *249*, 127–161. [CrossRef]

22. Khalloufi, M.; Valette, R.; Hachem, E. Adaptive Eulerian framework for boiling and evaporation. *J. Comput. Phys.* **2020**, *401*, 109030. [CrossRef]

23. Kumar Singh, N.; Premachandran, B. A coupled level set and volume of fluid method on unstructured grids for the direct numerical simulations of two-phase flows including phase change. *Int. J. Heat Mass Transfer* **2018**, *122*, 182–203. [CrossRef]
24. Hachem, E.; Rivaux, B.; Kloczko, T.; Digonnet, H.; Coupez, T. Stabilized finite element method for incompressible flows with high Reynolds number. *J. Comput. Phys.* **2010**, *229*, 8643–8665. [CrossRef]

25. Unverdi, S.O.; Tryggvason, G. A Front-Tracking Method for Viscous, Incompressible, Multi-fluid Flows. *J. Comput. Phys.* **1992**, *100*, 25–37. [CrossRef]

26. Sussman, M.; Smereka, P.; Osher, S. A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flows. *J. Comput. Phys.* **1994**, *114*, 146–159. [CrossRef]

27. Osher, S.; Fedkiw, R. Level Set Methods and Dynamic Implicit Surfaces; Applied Mathematical Sciences; Springer: Berlin/Heidelberg, Germany, 2004; Volume 153.

28. Feng, Q.; Johannsen, K. Experimental results of maximum transition boiling temperature during upflow in a circular tube at medium pressure. *Exp. Therm. Fluid Sci.* **1991**, *4*, 90–102. [CrossRef]

29. Ebrahim, S.A.; Chang, S.; Cheung, F.B.; Bajorek, S.M. Parametric investigation of film boiling heat transfer on the quenching of vertical rods in water pool. *Appl. Therm. Eng.* **2018**, *140*, 139–146. [CrossRef]

30. Liscic, B.; Tensi, H.M.; Luty, W. *Theory and Technology of Quenching*; Springer: Berlin/Heidelberg, Germany, 1992.

31. Kim, J.; Benton, J.F.; Wisniewski, D. Pool boiling heat transfer on small heaters: effect of gravity and subcooling. *Int. J. Heat Mass Transfer* **2002**, *45*, 3919–3932. [CrossRef]

32. Ikkene, R.; Koudil, Z.; Mouzali, M. Pouvoir de refroidissement des solutions de trempe à base de polymères hydrosolubles. *Comptes Rendus Chim.* **2008**, *11*, 297–306. [CrossRef]

33. Kopun, R.; Škerget, L.; Hriberšek, M.; Zhang, D.; Stauder, B.; Greif, D. Numerical simulation of immersion quenching process for cast aluminium part at different pool temperatures. *Appl. Therm. Eng.* **2014**, *65*, 74–84. [CrossRef]

34. Hachem, E.; Feghali, S.; Codina, R.; Coupez, T. Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation: a monolithic approach to fluid-structure interaction. *Int. J. Numer. Meth. Engng* **2013**, *94*, 805–825. [CrossRef]

35. Aslam, T.D. A partial differential equation approach to multidimensional extrapolation. *J. Comput. Phys.* **2004**, *193*, 349–355. [CrossRef]

36. Patankar, S.V. *Numerical Heat Transfer and Fluid Flow*; Series in Computational Methods in Mechanics and Thermal Sciences; Hemisphere Publ. Co.: New York, NY, USA, 1980.

37. Hachem, E. Stabilized Finite Element Method for Heat Transfer and Turbulent Flows Inside Industrial Furnaces. Ph.D. Thesis, Ecole Nationale Supérieure des Mines de Paris, Paris, France, 2009.

38. Hachem, E.; Kloczko, T.; Digonnet, H.; Coupez, T. Stabilized finite element solution to handle complex heat and fluid flows in industrial furnaces using the immersed volume method. *Int. J. Numer. Meth. Fluids* **2012**, *68*, 99–121. [CrossRef]

39. Santiago, R.D.; Hernández, E.M.; Otero, J.A. Constant mass model for the liquid-solid phase transition on a one-dimensional Stefan problem: Transient and steady state regimes. *Int. J. Therm. Sci.* **2017**, *118*, 40–52. [CrossRef]

40. Chaurasiya, V.; Singh, J. An analytical study of coupled convective heat and mass transfer with volumetric heating describing sublimation of a porous body under most sensitive temperature inputs: Application of freeze-drying. *Int. J. Heat Mass Transf.* **2023**, *214*, 124294. [CrossRef]

41. Chaurasiya, V.; Singh, J. Numerical investigation of a non-linear moving boundary problem describing solidification of a phase change material with temperature dependent thermal conductivity and convection. *J. Therm. Stress.* **2023**, *46*, 799–822. [CrossRef]

42. Chaurasiya, V.; Jain, A.; Singh, J. Numerical Study of a Non-Linear Porous Sublimation Problem With Temperature-Dependent Thermal Conductivity and Concentration-Dependent Mass Diffusivity. *ASME J. Heat Mass Transf.* **2023**, *145*, 072701. [CrossRef]

43. Chaurasiya, V.; Upadhyay, S.; Rai, K.; Singh, J. A temperature-dependent numerical study of a moving boundary problem with variable thermal conductivity and convection. *Waves Random Complex Media* **2023**, 1–25. [CrossRef]

44. Kim, D.G.; Jeon, C.H.; Park, I.S. Comparison of numerical phase-change models through Stefan vaporizing problem. *Int. Commun. Heat Mass* **2017**, *87*, 228–236. [CrossRef]

45. Denis, R. Modélisation et Simulation de l'effet Leidenfrost dans les Micro-Gouttes. Ph.D. Thesis, Université de Grenoble, Grenoble, France, 2014.

46. Ramezanzadeh, H.; Ramiar, A.; Yousefifard, M. Numerical investigation into coolant liquid velocity effect on forced convection quenching process. *Appl. Therm. Eng.* **2017**, *122*, 253–267. [CrossRef]

47. Khalloufi, M.; Mesri, Y.; Valette, R.; Massoni, E.; Hachem, E. High fidelity anisotropic adaptive variational multiscale method for multiphase flows with surface tension. *Comput. Methods Appl. Mech. Eng.* **2016**, *307*, 44–67. [CrossRef]

48. Shakoor, M.; Scholtes, B.; Bouchard, P.O.; Bernacki, M. An efficient and parallel level set reinitialization method—Application to micromechanics and microstructural evolutions. *Appl. Math. Model.* **2014**, *39*, 7291–7302. [CrossRef]

49. Denner, F.; van Wachem, B.G. Numerical time-step restrictions as a result of capillary waves. *J. Comput. Phys.* **2015**, *285*, 24–40. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.