# A SOLUTION TO STEFAN PROBLEM USING EULERIAN TWO FLUID VOF MODEL

Alen Cukrov, Yohei Sato, Ivanka Boras, Bojan Ničeno

► To cite this version:

Alen Cukrov, Yohei Sato, Ivanka Boras, Bojan Ničeno. A SOLUTION TO STEFAN PROBLEM USING EULERIAN TWO FLUID VOF MODEL. Brodogradnja, 2021, 72 (4), pp.141-164. 10.21278/brod72408. hal-04515158

HAL Id: hal-04515158

https://hal.science/hal-04515158

Submitted on 21 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Alen Cukrov
Yohei Sato
Ivanka Boras
Bojan Ničeno

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

http://dx.doi.org/10.21278/brod72408

ISSN 0007-215X
eISSN 1845-5859

# A SOLUTION TO STEFAN PROBLEM USING EULERIAN TWO-FLUID VOF MODEL

UDC 641.522.9:536.2:629.5.015.2

Original scientific paper

## Summary

A novel approach for the solution of Stefan problem within the framework of the multi-fluid model supplemented with Volume of Fluid (VOF) method, i.e. two-fluid VOF, is presented in this paper. The governing equation set is comprised of mass, momentum and energy conservation equations, written on a per-phase basis and supplemented with closure models via the source terms. In our method, the heat and mass transfer is calculated from the heat transfer coefficient, which has a fictitious function and depends on the local cell-size and the thermal conductivity, and the implementation is straightforward because of the usage of the local value instead of a global parameter. The interface sharpness is ensured by the application of the geometrical reconstruction scheme implemented in VOF. The model is verified for three types of computational meshes including triangular cells, and good agreement was obtained for the interface position and the temperature field. Although the developed method was validated only for Stefan problem, the application of the method to engineering problems is considered to be straightforward since it is implemented to a commercial CFD code only using a local value; especially in the field of naval hydrodynamics wherein the reduction of ship resistance using boiling flow can be computed efficiently since the method handles phase change processes using low-resolution meshes.

Key words: Stefan problem; two-fluid VOF; locally imposed closures; boiling flow; mesh versatility

## 1. Introduction

The sustainable development and design in industrial sectors concerned with high heat outputs, e.g., in power engineering and metallurgy, are strongly dependent on well-optimized processes that involve boiling phenomena. Thus, for example, in power engineering nucleate boiling is a favorable boiling regime because of high heat flow rates that can be achieved during cooling inside pressure vessels like boilers, thereby avoiding dangerous overheating of the heat transfer surface that may occur if the liquid coolant cannot reach that surface. This unfavorable boiling heat transfer regime is known as film boiling and, in contrary, is very common in metallurgy, where achievement of desired mechanical properties of metal material via
quenching process is followed by the occurrence of vapor blanked around the hot object due to
high heat flow rate released in the quenchant liquid.

In addition to these classical industrial appearances of boiling processes is the application
of the process in the field of naval architecture, where it can be addressed to, already
well-established, air lubrication systems used in reduction of ship resistance; and whose
theoretical foundations are given in Pavlov et al. [1]. To reduce the ship resistance, the bubble
injection technique under a vessel has been actively studied in recent years, and it was applied
to full-scale ships. For example, Kodama et al.[2] reported approximately 10-15% of saving
of the total energy consumption for an experimental ship, and Mizokami et al.[3] achieved at
most 12% fuel saving. In these systems, blowers are used to inject air below a hull, which
requires energy input. From the numerical point of view, two distinct approaches may be found
in the literature in this regard: the more fundamental ones, dealing with micro-bubble drag
reduction (MBDR); and the macroscopic studies aimed to study the flow features around the
bodies. An example of the former is the numerical simulation of an MBDR on a flat plate
conducted using the two-fluid model in Mohanarangam et al. [4]. In the latter case, a scaled
bulk carrier is a subject of the study by Sindagi et al. [5]. Using the single-fluid VOF method
the ship resistance is obtained for two different ship speeds; different air injection rates and
different diameters of the air injection holes. The references on previous work carried out in
both assessments, micro- and macro-scale can also be found therein.

If a heating system is implemented at the fore-bottom of the vessel and boiling takes place
there as illustrated in Fig. 1, the reduction of ship resistance due to the bubbles and the low
viscosity of water with higher temperature may be achieved. The vapor bubbles condensate in
the downstream because of lower temperature, but the non-condensable gas bubble may remain.
As the heat source, the heat energy from condenser in a steam power plant on board could be
used in case that the plant is equipped. Recent research shows that a few volume percent of
vapor bubble can reduce up to 45% of drag [6].

<figure><img src="image_3.png" /><figcaption>Fig. 1 Conceptual sketch of the ship resistance reduction due to boiling flow.</figcaption></figure>

However, the publications that study the latter effect, that is, the reduction of the ship
resistance involving boiling flow using the numerical simulation were not found in the
literature.

To optimize such processes efficiently, reliable numerical simulations are required to
estimate boiling heat transfer using affordable computational resources and time. Furthermore,
the numerical simulations have to be capable to compute various flow situations and to reflect
properly the physics underlying the boiling phenomena in the sense that the dominant mass
transfer rate occurs at the interface between the phases. The development of computational
methods for application in the former case (nucleate boiling) has been subject of extensive
studies in the past, simulating thereby the boiling curve from the side of nucleate boiling on.
This approach is utilized, for example, in Sato and Ničeno [7]. In the study presented by the
authors, the boiling curve is being modeled from nucleate to film boiling, with consideration of the Critical Heat Flux (CHF), conjugate heat transfer and bubble micro-layer. However, the reliable computational models are still needed for the latter case (film boiling), where the transition from the initial film boiling to succeeding nucleate boiling regime needs to be modeled. Therefore, the long-term goal of our research is to develop the computational method which can simulate the transition from film boiling to nucleate boiling.

In development of such models, the first step is the verification of Stefan problem: the one-dimensional phase-change problem. This moving boundary problem provides the information about unsteady evolution of liquid-vapor interface and the corresponding temperature field throughout the domain composed of those phases, thus being suitable for verification of the applied interface tracking algorithm and accompanied mass transfer model. Stefan problem has been typically used as a benchmark in the boiling process, which consists of superheated vapor and saturated liquid. The examined numerical methods are the interface tracking methods, such as, for example, single-fluid Volume-of-Fluid (VOF), coupled with Direct Numerical Simulation (DNS) of turbulence, in which the volume fraction of vapor/liquid is tracked by solving the advection equation. Sun *et al.* [8] proposed the energy jump model, and they obtained solutions for different density ratios, assuming zero thermal conductivity of phase at saturation temperature and equal specific heat capacities of both phases. Stefan problem is used in verification of VOF-based phase change model proposed by Pan *et al.* [9]. The model is based on fixing the interface temperature to saturation temperature and was verified for three different working fluids: HFE7100, R113 and water. A systematic introduction in the solution of Stefan problem using both, the analytical and numerical approach, was given in Perez-Raya and Kandlikar [10]. The authors presented the necessary steps for obtaining the solution using VOF and presented the theoretical framework that underlies final expressions used for verification purposes. The accuracy of three different mass transfer models, namely: Lee [11], Rattner [12] and Sun model [13], in solving Stefan problem was investigated by Kim *et al.* [14]. The transient interface evolution was examined therein, using different mesh resolutions, time increment and, in Lee model, model constants. Perez-Raya and Kandlikar [15] used the sharp interface approach for the solution of Stefan problem. In their approach, the saturation temperature is imposed in the interface cells and the estimation of temperature gradient is calculated from the distance between the interface and neighboring cell centers. More recently, the idea of incorporation of a gradient of volume fraction field into a model that relies on empirical constant is proposed in Chen *et al.* [16]. By considering the cell size, the Chen's approach can be applied to non-uniform quadrilateral mesh. The other studies that utilize interface tracking approach, and some of them apply it in solution of Stefan problem, can be found in the review of Kharangate and Mudawar [17].

Compared to the abovementioned single-fluid models, a two-fluid model is considered to be able to predict more accurate results on coarser mesh owing to the unnecessity of resolving the interface geometry/shape, though two sets of momentum and energy equations must be solved. However, due to the averaging process in the formulation of the governing equations, the shape and area of the liquid-vapor interface are lost, and the dynamics related to the interface, e.g. surface tension and interphase forces, cannot be directly computed. To mitigate this, several approaches to handle interfaces inside the two-fluid model were proposed; a review can be found in Mer *et al.* [18]. Among the modeling approaches presented in [18], further exploration in the paper is made on the following models, namely: Large Interface Model (LIM), Generalized Large Interface Model (GLIM) and Large Bubble Model (LBM). LIM and LBM were applied in the solution of Stefan problem by Fleau [19]. The description of LIM is given in Coste [20]. LIM identifies the large interface, but does not reconstruct it. In order to determine the distance from interface to the adjacent cell, the method utilizes a three-cell stencil
concept: one cell contains the interface, while other two are occupied by single representatives of the involved phases. The general idea behind LBM can be found in Dènefle et al. [21]. LBM is three- or four-field approach in the framework of two-fluid model, that considers two continuous and one or two dispersed phases. The distinction between the continuous phases, that are simulated, and dispersed phases, that are modeled, is made with a spatial cutting length, while the mass exchange between the continuous and dispersed field of the same phase is established via the mass transfer terms. However, to the authors' knowledge, there appears not to exist evaluation of the two-fluid VOF for Stefan problem in the literature.

In this paper, the two-fluid VOF model available in CFD package Ansys® Fluent, Release 20.1 as Multi-Fluid VOF model is adopted for the solution of Stefan problem. The advantage of the two-fluid VOF over the single-fluid VOF is in avoidance of averaging of the material properties in the energy equation, which reduces the accuracy in the heat-flux calculation [15].

The novelty of this paper is the application of two-fluid VOF method to Stefan problem and the validation of the model to not only quadrilateral meshes but also triangular meshes. The originality of our model is the usage of local mesh size for the heat transfer coefficient, which is explained in the section of the numerical model. The paper is organized as follows. The theoretical solution of Stefan problem is presented in Chapter 2. The proposed computational method is the content of Chapter 3. Chapter 4 outlines the details related to performed numerical simulations. The obtained results together with the applied error measures are discussed in Chapter 5.

## 2. Theoretical solution of Stefan problem

The analytical solution of Stefan problem, presented here, refers to expressions for computation of interface position and temperature distribution during planar interface evaporation in the system composed of superheated vapor and saturated liquid, reinterpreting thereby the work of Perez-Raya and Kandlikar [10]. The equations involved in derivation process are the energy balance at the interface and one-dimensional heat conduction equation in a continuum. In what follows, a concise exposition of the theoretical framework underlying the applied analytical expressions is given.

### 2.1 Conditions at the interface

In the analytical solution used in this study, the interface condition that accounts only for occurrence of heat conduction in the system is incorporated. This is because in the case with superheated vapor and saturated liquid the heat transfer takes place only in vapor phase,

$$ \frac{\partial T_l}{\partial x} = 0 \qquad (1) $$

Hence, the interface condition is given by:

$$ X'(t)\rho_v h_{fg} = -k_v \frac{\partial T}{\partial x}\bigg|_{X(t)} \qquad (2) $$

where $X'(t)$ is the interface velocity, $\rho_v$ is the vapor density, $h_{fg}$ is the latent heat, $k_v$ is the thermal conductivity of the vapor, while the negative signed gradient denotes heat flux in normal direction to the interface. Therefore, the mass generated in phase change process is proportional to the balance of the heat fluxes applied from each phase in normal direction to the interface.
## 2.2 Interface position

In order to derive analytical model for this moving boundary problem, introduction of similarity variable, $\eta$, is required to relate the spatial coordinate, $x$, and time, $t$, as:

$$ \eta = x/\sqrt{t} \qquad (3) $$

Written in terms of the interface position, that is, $x = X(t)$ and $\eta = \xi$, this relation yields the equation for determination of interface displacement:

$$ X(t) = \xi\sqrt{t} \qquad (4) $$

where $\xi$ is the constant that has to be determined using iterative procedure. The derivation of Eq. (4) with respect to time gives the interface velocity:

$$ X'(t) = \frac{\xi}{2\sqrt{t}} \qquad (5) $$

The substitution of Eq. (5) into Eq. (2) and involvement of solution that is obtained from incorporation of similarity variable in one-dimensional unsteady heat conduction equation gives:

$$ \lambda_v h_{fg} = \frac{c_v(T_w - T')}{\sqrt{\pi} \exp(\lambda_v^2) \operatorname{erf}(\lambda_v)} \qquad (6) $$

where $\lambda_v$ is constant that has to be determined using an iterative procedure, $c_v$ is the specific heat capacity of the vapor, $T_w$ is the temperature of the superheated wall and $T'$ is the saturation temperature. The constant $\lambda_v$ embeds the value of similarity variable at the interface, that is, the constant $\xi$, and is defined as:

$$ \lambda_v = \frac{\xi}{2\sqrt{\alpha_v}} \qquad (7) $$

where $\alpha_v$ is the thermal diffusivity of the vapor phase. Substitution of Eq. (7) into Eq. (4) yields the instantaneous interface position in the case with superheated vapor and saturated liquid [10]:

$$ X(t) = 2\lambda_v\sqrt{\alpha_v t} \qquad (8) $$

In the cases where initial vapor layer is needed for computational purposes, the necessary time shift of the simulation time is accomplished using Eq. (8).

## 2.3 Temperature field

The analytical solution of the temperature field inside the system composed of two phases, vapor and liquid, distinguished by a sharp interface, is derived from the partial differential equation (PDE) for one-dimensional unsteady heat conduction in a single continuum. The time and space coordinates present therein are transformed using the similarity variable introduced in the previous section. Two boundary conditions are required to obtain the
solution of the PDE: at the interface and at the wall. The temperature distribution at certain time instance, t, is computed as:

$$T_a(x, t) = T_w - (T_w - T') \frac{\operatorname{erf}\left(\frac{x}{2\sqrt{\alpha_v t}}\right)}{\operatorname{erf}(\lambda_v)} \quad (9)$$

### 3. Numerical method

Within the framework of two-fluid model adopted here for solution of laminar two-phase flow with phase change, two separate sets of conservation equations are solved in each computational cell; one for each phase present in the domain. The theoretical foundations underlying the two-fluid model are given in Drew and Passmann [22], while the model equations available in Fluent and applied in this study are taken from [23]. The communication between the two phases is established via source terms in the equations, referred to as phase interaction terms. In addition, a single pressure field is shared between the phases. The temperature in the liquid phase is assumed to be constant at the saturation temperature.

#### 3.1 Conservation equations

The evolution of planar liquid-vapor interface due to evaporation, studied here, is described with three laws of classical physics, namely: the conservations of mass, momentum and energy:

$$\frac{\partial}{\partial t}(\alpha_q \rho_q) + \nabla \cdot (\alpha_q \rho_q \vec{v}_q) = \dot{m}_{pq} \quad (10)$$

$$\frac{\partial}{\partial t}(\alpha_q \rho_q \vec{v}_q) + \nabla \cdot (\alpha_q \rho_q \vec{v}_q \vec{v}_q) = -\alpha_q \nabla p + \nabla \cdot \bar{\tau}_q + \vec{R}_{pq} + \dot{m}_{pq} \vec{v}_{pq} \quad (11)$$

$$\frac{\partial}{\partial t}(\alpha_q \rho_q h_q) + \nabla \cdot (\alpha_q \rho_q \vec{u}_q h_q) = \alpha_q \frac{dp_q}{dt} + \bar{\tau}_q \cdot \nabla \vec{u}_q + Q_{pq} + \dot{m}_{pq} h_{pq} \quad (12)$$

where subscripts q and pq refer to particular phase and the phase pair involved in certain interfacial transfer process, respectively. Appropriate closure models are required for the interphase transfer terms that appear at the right-hand side in these equations, those are, interfacial mass transfer, $\dot{m}_{pq}$, the momentum transfer, $\vec{R}_{pq}$, and the heat transfer, $Q_{pq}$, terms. These terms stem from averaging process used in derivation of the set of conservation equations and more details regarding the closures applied therein is given in the forthcoming section.

#### 3.2 Closure models

##### 3.2.1 Interfacial area

The interphase transfer terms in the presented set of conservation equations for a continuum composed of interpenetrating phases have to be supplemented with suitable correlations. A common issue, however, in all the terms is related to modeling of the interfacial area, a flow regime dependent ratio of dispersed phase's projected area and cell volume. When free surface flow modeling is considered, the interfacial area per unit volume is computed as:

$$A_i = |\nabla \alpha_q| \quad (13)$$

This idea is also inherited for free surface modeling within Algebraic Interfacial Area Density (AIAD) modeling approach for handling different interface scales inside two-fluid
model, presented by Höhne and Vallée [24]. Since the application of Eq. (13) yields the application of closure models only within the interface zone, the closures used in this study can be referred as “locally imposed closures”.

### 3.2.2 Interphase momentum transfer

In general, the interphase momentum transfer consists of the drag force, lift force, wall lubrication force, turbulent dispersion force and virtual mass force. Since in the solution of Stefan problem there is a distinction between the phases with a free surface, the incorporation of drag force is necessary according to the study of stratified two-fluid modeling by Štrubelj et al. [25]. The role of drag force, as noted by the authors, is to ensure equal velocities of the phases within the interface region, and the physical basis of the correlation applied for computation of the drag force coefficient is not mandatory since the drag force term has no real physical meaning in a stratified free-surface flow as it has in a dispersed flow – its presence in the momentum equation is due to averaging process used in derivation of two-fluid model. Thus, the drag force is only taken into account in this paper. The general formulation of interphase drag is given by:

$$ \vec{R}_{pq} = K_{pq}(\vec{u}_p - \vec{u}_q) \qquad (14) $$

where $K_{pq}$ is drag force coefficient and $\vec{u}_p - \vec{u}_q$ is the relative velocity between the phases. In this study, the drag force coefficient is modeled using the default anisotropic drag model available in Fluent. This model allows for higher drag in the normal direction to the interface and lower drag in the direction tangential to the interface, with the default values of normal and tangential interfacial drag friction factors equal to 1e+6 and 1e+3, respectively. In addition, application of this model contributes to the stability of the computation [26].

### 3.2.3 Interphase heat transfer

The relation for computation of heat transfer across the interface between the phases reads:

$$ Q_{pq} = h_{pq}A_i(T_p - T_q) \qquad (15) $$

where $h_{pq}$ is the heat transfer coefficient, $A_i$ is the interfacial area and $T_p - T_q$ is the temperature difference between the interacting phases. The heat transfer coefficient at the interface is given by:

$$ h_{pq} = h_{qp} = \frac{\lambda_q N u_p}{d_p} \qquad (16) $$

where $\lambda_q$ is the thermal conductivity of the liquid, while $Nu_p$ and $d_p$ are Nusselt number in the vapor phase and the user-prescribed characteristic dimension, respectively. In the present study, the heat transfer coefficient, $h_{pq}$, is modeled using two-resistance model; wherein the different correlation or expression of the heat transfer coefficient is applied at each side of the vapor-liquid interface.

At the vapor side, the heat transfer coefficient is computed using cell-size dependent Nusselt number correlation, which we proposed in this paper:
$$
Nu_p = \frac{2 d_p}{\Delta x} \tag{17}
$$

where $d_p$ is the user-prescribed vapor phase characteristic dimension and $\Delta x$ is the approximated cell size. The derivation of Eq. (17) is given in Appendix. Note that the vapor phase diameter, $d_p$, defined in Eq. (17) is cancelled, when it is substituted into Eq. (16). Thus, an arbitral value can be given as $d_p$. The cell size is calculated as:

$$
\Delta x = \sqrt[3]{V_c} \tag{18}
$$

where $V_c$ is the cell volume.

At the liquid side, the saturation temperature is maintained at the interface using zero-
resistance model. This model refers to application of a large heat transfer coefficient, $h_l \to +\infty$, to ensure the saturation temperature at the interface throughout the computation process. This idea corresponds to application of large coefficient method to fix the saturation temperature at the interface, as discussed in the work from Perez-Raya and Kandlikar [15].

3.2.4 Interphase mass transfer

The mass transfer at liquid-vapor interface is modeled using the thermal phase change
model implemented in Fluent that reads:

$$
\dot{m}_{lv} = - \frac{c_v h_v A_i (T_v - T') + c_l h_l A_i (T_l - T')}{h_{fg}} \quad (19)
$$

where $A_i$ is the interfacial area computed with Eq. (13), $T'$ is saturation temperature, $h_{fg}$ is the latent heat of vaporization, $C_v$ and $C_l$ are, respectively, vapor and liquid phase's scaling factors, $h_v$ and $h_l$ are, respectively, the heat transfer coefficients at vapor and liquid side of the interface, while $T_v$ and $T_l$ are temperatures of vapor and liquid phase, respectively.

In Eq. (19), $C_v$ and $C_l$ equal to unity theoretically, although they are treated as the tuning parameters in Fluent. Since we intend to obtain a solution of the governing equations without introducing any tuning parameters except computational mesh, they are set to be unity in this paper.

3.3 Discretization of equations

The interpolation near the interface is accomplished using Geo-Reconstruct scheme.
Furthermore, the momentum and energy equations are discretized using second-order upwind
interpolation scheme. The pressure equation is solved with PRESTO! method [23], while the
gradient discretization is carried out using least square cell based method. The coupling between
pressure and velocity fields is realized using phase coupled SIMPLE algorithm. The governing
equations are discretized in time using the implicit first-order Euler method.
# 4. Conditions of simulations

## 4.1 Material properties

Since our final objective of the simulations is an application to engineering problems, we select the water as the working fluid instead of an artificial material. The properties of the water at 1.013 *bar* are listed in Table 1. The saturation temperature is 373.15 *K* at this system pressure.

## 4.2 Computational domain and boundary condition

The computational domain is designed as a beam divided into finite number of cells in the *x*-direction of the two-dimensional (2-D) Cartesian coordinate system and with one cell layer in the *y*-direction, i.e. a special case of 2-D domain that is also known as 1.5-D mesh. A schematic view of the computational domain together with the boundary condition and initial phase distribution is shown in Fig. 2.

<table><caption>**Table 1** The properties of vapor and liquid at 1.013 bar [27].</caption><thead><tr><th></th><th>Density<br>[kg/m³]</th><th>Specific heat<br>capacity<br>[kJ/(kg K)]</th><th>Dynamic<br>viscosity<br>[µNs/m²]</th><th>Thermal<br>conductivity<br>[W/(m K)]</th><th>Latent heat of<br>vaporization<br>[kJ/kg]</th></tr></thead><tbody><tr><td>Vapor</td><td>0.597</td><td>2.03</td><td>12.55</td><td>0.025</td><td rowspan="2">2256.7</td></tr><tr><td>Liquid</td><td>958.3</td><td>4.22</td><td>277.53</td><td>0.679</td></tr></tbody></table>

The left boundary of the computational domain is the wall, where the constant superheat of 10 *K*, $\Delta T$, is prescribed for all the simulation cases. The symmetry condition is imposed on the top and bottom boundaries to ensure the one-dimensionality. The right boundary is defined as the pressure outlet, with prescribed value of volume fraction and the bulk liquid's temperature.

<figure><img src="image_4.png" /><figcaption>**Fig. 2** A schematic representation of the domain accompanied with the applied boundary conditions and initial phase distribution.</figcaption></figure>

## 4.3 Initial condition

The vapor volume fraction is set to unity in the first cell adjacent to the superheated wall, while in the remainder of the domain is set to zero; a vice versa approach is applied for liquid volume fraction as depicted in Fig. 3a. Here $\Delta x_c$ is the cell length of the wall adjacent cell. The velocity fields of both phases are initially set to zero in all the cells throughout the domain. The initial vapor temperature inside the vapor phase is defined by linear profile changing from $T'+\Delta T$ to $T'$ as illustrated in Fig. 3b. The initial liquid temperature is, on the other hand, set to $T'$ in the whole domain.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** The initial condition of the vapor volume fraction (a) and the vapor temperature (b).</figcaption></figure>

4.4 Computational mesh

The presented numerical model for solution of Stefan problem using two-fluid VOF model has been assessed to three different types of computational meshes: (a) uniform quadrilateral, (b) stretched quadrilateral and (c) the hybrid of stretched quadrilateral and triangular cells (hereinafter also named hybrid mesh), as shown in Fig. 4, respectively. The mesh is generated by using the GMSH code [28].

<figure><img src="image_6.png" /><figcaption>**Fig. 4** Types of meshes used in computations: (a) uniform quadrilateral, (b) stretched quadrilateral and (c) hybrid of stretched quadrilateral and triangular cells.</figcaption></figure>

The first mesh type is referred as the standard mesh, while the last two types serve to
prove the capability of the proposed model under more realistic conditions regarding
computational domain, that is, when the gradients of the dependent variables are computed
using graded mesh in wall adjacent region. With respect to mesh resolution, the model
verification is carried out using meshes with 50, 100 and 200 cells, hereinafter also named
coarse, medium and fine, respectively. The graded mesh, where included, is composed of ten
quadrilateral cells with expansion ratio of 1.2, while the number of cells with the dominant cell-
shape has been varied.

The time steps applied in the simulations are mesh resolution dependent and written in Table 2. The total simulation time is divided in two parts. Firstly, the simulation is run for 0.1 ms to obtain more accurate initial distribution of interface displacements; otherwise, an initial discrepancy from the analytical solution is observed that shifts the normal interface evolution to later time. Thus, two-orders of magnitude lower time step size than in the
remainder of simulation is used together with the number of time steps that is 10 % of the number of time steps in the second part of the simulation. Then, the simulation is run for 0.1 s with increased time step size. Furthermore, due to presence of initial vapor cell in the simulation, the simulation time has to be shifted for a constant value obtained from Eq. (8). Therefore, the results at specific time instances are obtained using linear interpolation from neighboring values.

<table><caption>**Table 2** A general, mesh size dependent, time stepping procedures applied for the uniform and stretched meshes.</caption><thead><tr><th rowspan="3">Number of cells</th><th colspan="4">Simulation time</th><th rowspan="3">Total simulation time [s]</th></tr><tr><th colspan="2">First part</th><th colspan="2">Second part</th></tr><tr><th>Number of time steps</th><th>Time step size [s]</th><th>Number of time steps</th><th>Time step size [s]</th></tr></thead><tbody><tr><td>50</td><td>20</td><td>5.00e-6</td><td>200</td><td>5.00e-4</td><td rowspan="3">0.1</td></tr><tr><td>100</td><td>40</td><td>2.50e-6</td><td>400</td><td>2.50e-4</td></tr><tr><td>200</td><td>80</td><td>1.25e-6</td><td>800</td><td>1.25e-4</td></tr></tbody></table>

The simulations carried out using medium and fine hybrid meshes required, however, finer time step sizes, listed in Table 3.

## 4.5 Implementation into Fluent

Since the two-fluid modeling approach utilized the concept of interpenetrating continua, in each cell of the computational mesh the relevant dependent variable's has to be assigned for both phases. Thus, three field variables are initialized on a per phase basis using Define-Init User-Defined Function (UDF) macro: volume fraction, velocity and temperature. Furthermore, the proposed correlation for computation of the heat transfer coefficient, Eq. (16) is incorporated into the Fluent code using Define-Exchange-Property UDF macro. Finally, in the cases where the hybrid mesh is used, the interface position is calculated using Define-Report-Definition-FN UDF macro.

<table><caption>**Table 3** Time stepping procedures for the cases of certain hybrid meshes.</caption><thead><tr><th rowspan="3">Number of cells</th><th colspan="4">Simulation time</th><th rowspan="3">Total simulation time [s]</th></tr><tr><th colspan="2">First part</th><th colspan="2">Second part</th></tr><tr><th>Number of time steps</th><th>Time step size [s]</th><th>Number of time steps</th><th>Time step size [s]</th></tr></thead><tbody><tr><td>100</td><td>80</td><td>1.25e-6</td><td>800</td><td>1.25e-4</td><td rowspan="2">0.1</td></tr><tr><td>200</td><td>100</td><td>1.00e-6</td><td>1000</td><td>1.00e-4</td></tr></tbody></table>

In all the computations, the default relaxation factors for the iterative calculations implemented in Fluent were used.

# 5. Results and discussions

## 5.1 Definition of the interface position

In the computations performed with the meshes composed of quadrilateral cells, the interface is considered as an iso-surface of vapor volume fraction equal to 0.5, that is, $\alpha_v = 0.5$. Thus, the information on interface position is obtained by tracking the x-coordinate of the iso-surface. In the case of hybrid mesh, inside the mesh region composed of triangular cells, this definition of interface, however, leads to interface distributed over two or more cells in a
piecewise-linear segments. Therefore, a cell-shape independent approach is used, with the interface position given by:

$$X(t) = \frac{\sum_{i=1}^{n} \alpha_{v,i} V_{c,i}}{\sum_{i=1}^{n} V_{c,i}} L \qquad (20)$$

where $\alpha_{v,i}$ is the vapor volume fraction in i-th cell, $V_{c,i}$ is i-th cell's volume and L is the length of the domain. The sum in denominator of Eq. (20) is the volume of the domain.

## 5.2 Definition of computational errors

We define the computational errors of the interface position and the temperature field here. The relative error of the interface position at time t is given by:

$$E_{X_{int}}(t) = \frac{|X_{int,a}(t) - X_{int,n}(t)|}{X_{int,a}(t)} \qquad (21)$$

where $X_{int,a}(t)$ and $X_{int,n}(t)$ refer to the x-coordinate of the interface position obtained from the analytical solution and the numerical simulation, respectively. In this study, the relative error of the interface position is calculated for two time instances: in the middle of the simulation, $t = t_{end}/2$, and at the end of the simulation, $t = t_{end}$.

We evaluate the relative error in the temperature field at the end of the simulation time which is defined as:

$$E_T = \max(E_T(x_{c1}, t_{end}), E_T(x_{c2}, t_{end}), \dots, E_T(x_{cn}, t_{end})) \qquad (22)$$

where $E_T(x_{ci}, t_{end})$ is the error of the temperature at the cell center of ci-th cell calculated by:

$$E_T(x_{ci}, t_{end}) = \frac{|T_a(x, t_{end}) - T_n(x_{ci}, t_{end})|}{\Delta T_w} \qquad (23)$$

where $T_a(x, t_{end})$ and $T_n(x, t_{end})$ are the temperatures at the location of the cell center, x, for the analytical solution and the numerical simulation, respectively, and $\Delta T_w$ is wall superheat. The temperature of the vapor is used for $T_n$, although in the two-fluid modeling approach two energy equations are solved and two temperature fields are present in the domain. This is because the temperature field calculated from the energy equation of the liquid phase is constant at the saturation temperature.

## 5.3 Uniform mesh

The advancement of vapor-liquid interface in the performed numerical simulation reproduces the analytical solution with a commensurate level of accuracy. This is shown in Fig. 5a. Additionally, the obtained result is supplemented with ± 2 % error band in order to make a stringent estimation of solution accuracy, as indicated in Figs. 5b and 5c.
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Evolution of the interface position using the uniform quadrilateral mesh: (a) for the whole simulation time, (b) the detailed view at the beginning, and (c) near the end of the simulation.</figcaption></figure>

The comparison between Fig. 5b and Fig. 5c, together with exact values given in Table 4, points out that increased mesh resolution does not necessarily provide the result that is more accurate. This is in accordance with the study in Gauss *et al.* [29], where is shown that, when two-fluid model is used in typical interface tracking case, the solution is not fully dependent on the mesh size. Here, close to the start of the computation, within the period from 10 % to 15 % of the total simulation time (Fig. 5b), the solution obtained on fine mesh exhibits the highest accuracy. However, when the end of the simulation is being approached, that is, 85 % to 90 % of the prescribed total time (Fig. 5c), the result obtained using medium mesh is found as the most accurate.

<table><caption>**Table 4** The relative computational error in estimation of interface position at selected two time instances for different resolutions of uniform quadrilateral mesh.</caption><thead><tr><th rowspan="2">Mesh type</th><th rowspan="2">Number of cells</th><th colspan="2">Relative error in interface position [%]</th></tr><tr><th>t = 0.05 s</th><th>t = 0.1 s</th></tr></thead><tbody><tr><td rowspan="3">Uniform mesh</td><td>50</td><td>0.54</td><td>1.38</td></tr><tr><td>100</td><td>0.28</td><td>1.24</td></tr><tr><td>200</td><td>0.91</td><td>1.70</td></tr></tbody></table>

The temperature field inside the domain at $t = 0.1$ s exhibit distributions shown in Fig. 6. Near the interface a discrepancy from exact, linear, temperature profile is observed. In this case the increased mesh resolution leads to more accurate result, as quantified in Table 5.
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The comparison of temperature at $t = 0.1$ s between the analytical solution and the simulation: whole domain (left) and the closer view near the interface (right).</figcaption></figure>

<table><caption>**Table 5** The maximum relative computational error in estimation of vapor phase's temperature field inside the domain at $t = 0.1$ s for uniform quadrilateral mesh.</caption><thead><tr><th>Mesh type</th><th>Number of cells</th><th>Maximum relative error in temperature field [%]</th></tr></thead><tbody><tr><td rowspan="3">Uniform mesh</td><td>50</td><td>7.49</td></tr><tr><td>100</td><td>4.38</td></tr><tr><td>200</td><td>3.44</td></tr></tbody></table>

Since the solution of Stefan problem with the approach presented in this paper differs from standard solutions obtained with single-fluid approach, it is noteworthy to outline some additional features of the model.

The consideration of liquid phase at saturation temperature imposes as the primary goal to study the vapor phase's temperature field. However, due to application of two-fluid model, there is also a temperature field associated with liquid phase that is distributed throughout the domain, as shown in Fig. 7. Fixing the liquid phase's temperature at the interface to saturation temperature by application of zero-resistance heat transfer model ensures stagnant temperature profile of the liquid, present also in the cell adjacent to interface zone with $\alpha_l = 0$ that is preceded by a jump in temperature decrease.
<figure><img src="image_9.png" /><figcaption>**Fig. 7** The vapor and liquid phase's temperature distributions in the domain composed of 100 quadrilateral, uniformly distributed, cells compared with analytical solution and accompanied with distribution of vapor volume fraction: whole domain (left) and the closer view near the interface (right); the markers denote cell center values.</figcaption></figure>

The tendency to equalize the velocities of the phases by application of the appropriate drag model, mentioned before, is shown in Fig. 8. The close agreement of the maximum velocity magnitude of the vapor (square marked) and liquid phase's velocity magnitude is outlined in the interface zone.

The liquid phase's temperature distribution and the velocity magnitude field of both phases are general features of the model and, therefore, are not considered further in this study.

## 5.4 Non-uniform meshes

### 5.4.1 Stretched quadrilateral mesh

The incorporation of stretched mesh in near wall region of a quadrilateral mesh also results in accurate estimation of interface displacement, as shown in Fig. 9a. Although shortly after the simulation start only the result obtained on the fine mesh enters the prescribed error band, as displayed in Fig. 9b, finally, as shown in Fig. 9c, it is accompanied with the coarse mesh solution.

According to exact measures, at two selected time instances, as indicated in Table 6, the interface displacement solution obtained with fine stretched quadrilateral mesh is more accurate than the solution obtained with quadrilateral uniform mesh of same number of cells. The medium mesh that yields most accurate interface displacement in the selected time instances in the case of uniform cell distribution, however, in this case has highest discrepancy, but still in reasonable limit, that is, lower than 5 %.
<figure><img src="image_10.png" /><figcaption>**Fig. 8** The distribution of vapor and liquid phase's velocity magnitude field in the domain composed of 100 quadrilateral, uniformly distributed, cells together with the profile of vapor volume fraction: whole domain (left) and the closer view near the interface (right); the markers denote cell center values.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9** Evolution of the interface position using the stretched quadrilateral mesh: (a) for the whole simulation time, (b) the detailed view at the beginning, and (c) near the end of the simulation.</figcaption></figure>
<table><caption>**Table 6** The relative computational error in estimation of interface position at selected two time instances for different resolutions of stretched quadrilateral mesh.</caption><thead><tr><th rowspan="2">Mesh type</th><th rowspan="2">Number of cells</th><th colspan="2">Relative error in interface position [%]</th></tr><tr><th>t = 0.05 s</th><th>t = 0.1 s</th></tr></thead><tbody><tr><td rowspan="3">Stretched quad.</td><td>50</td><td>0.97</td><td>0.69</td></tr><tr><td>100</td><td>2.25</td><td>4.22</td></tr><tr><td>200</td><td>0.80</td><td>1.65</td></tr></tbody></table>

The temperature field obtained using stretched quadrilateral meshes, shown in Fig. 10, has the maximum relative error in computation of temperature field lower than 10 %, when all the applied mesh resolutions are considered (Table 7). In the case of medium mesh, however, the better results are obtained in the case without stretched cells in wall adjacent zone.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** The comparison of temperature at $t = 0.1$ s between the analytical solution and the simulation: whole domain (left) and the closer view near the interface (right).</figcaption></figure>

<table><caption>**Table 7** The maximum relative computational error in estimation of vapor phase's temperature field inside the domain at $t = 0.1$ s for different resolutions of stretched quadrilateral mesh.</caption><thead><tr><th>Mesh type</th><th>Number of cells</th><th>Maximum relative error in temperature field [%]</th></tr></thead><tbody><tr><td rowspan="3">Stretched quad.</td><td>50</td><td>5.43</td></tr><tr><td>100</td><td>6.19</td></tr><tr><td>200</td><td>3.45</td></tr></tbody></table>

### 5.4.2 Hybrid mesh

The interface displacement in the case of hybrid mesh is calculated using Eq. (20) and exhibits the distribution shown in Fig. 11. Although detailed view in Fig. 11b suggests the solution obtained using the fine mesh as the most promising one, the result in Fig. 11c together
with the exact values of the interface displacement error, given in Table 8, indicate the reliability of medium size mesh.

<figure><img src="image_13.png" /><figcaption>**Fig. 11** Evolution of the interface position using the hybrid mesh: (a) for the whole simulation time, (b) the detailed view at the beginning, and (c) near the end of the simulation.</figcaption></figure>

<table><caption>**Table 8** The maximum relative computational error in estimation of vapor phase's temperature field inside the domain at $t = 0.1$ s for different resolutions of hybrid mesh.</caption><thead><tr><th rowspan="2">Mesh type</th><th rowspan="2">Number of cells</th><th colspan="2">Relative error in interface position [%]</th></tr><tr><th>t = 0.05 s</th><th>t = 0.1 s</th></tr></thead><tbody><tr><td rowspan="3">Hybrid mesh</td><td>50</td><td>0.31</td><td>2.36</td></tr><tr><td>100</td><td>2.09</td><td>3.06</td></tr><tr><td>200</td><td>2.88</td><td>3.79</td></tr></tbody></table>

Sampling the temperature field data from cell's centers yields a “zig-zag” temperature field distribution inside the zone composed of triangular cells in the hybrid mesh, as shown in Fig. 12. The error in computation of temperature field for coarse and medium resolution mesh is in the range of the errors already reported for other types of meshes with these resolutions, that is, less than 10 %. However, the maximum relative temperature field error obtained using fine mesh at the simulation end time is, compared to fine mesh calculations with other mesh types, the largest (Table 9).

The distribution of the volume fraction of the vapor phase and the temperature for the hybrid mesh at $t = 0.1$ s are shown in Fig. 13. The volume fraction sharply changed from zero to one, and the temperature field features linear change from the wall temperature to the saturation temperature in the vapor phase.
<figure><img src="image_14.png" /><figcaption>**Fig. 12** The comparison of temperature at $t = 0.1$ s between the analytical solution and the simulation: whole domain (left) and the closer view near the interface (right).</figcaption></figure>

<table><caption>**Table 9** The maximum relative computational error in estimation of vapor phase's temperature field inside the domain at $t = 0.1$ s for different resolutions of hybrid mesh.</caption><thead><tr><th>Mesh type</th><th>Number of cells</th><th>Maximum relative error in temperature field [%]</th></tr></thead><tbody><tr><td rowspan="3">Hybrid mesh</td><td>50</td><td>5.49</td></tr><tr><td>100</td><td>5.32</td></tr><tr><td>200</td><td>4.96</td></tr></tbody></table>

## 6. Conclusions

The reduction of ship resistance via air injection beneath the ship has been extensively studied in the past, but the corresponding idea of ship resistance reduction using boiling flow is relatively new in the literature, and the numerical assessments in this regard are still missing. To overcome this, a novel method for computation of mass transfer is proposed herein. In this paper, a novel method for computation of Stefan problem was proposed and the following conclusions can be drawn:

1. Stefan problem was accurately solved using the two-fluid VOF model. To this end, the local-parameter based, heat and mass closure was derived, under the assumption of the constant liquid temperature at the saturation.

2. The sharp distribution of the volume fraction of liquid/vapor was achieved by using the geometrical reconstruction scheme implemented in the VOF model.

3. The velocity jump was captured at the interface owing to the sharp distribution of the interfacial area density calculated as the volume fraction gradient.

4. The method is capable to solve Stefan problem on the hybrid mesh predominantly composed of triangular cells with a commensurate level of accuracy. The accurate solutions are also obtained for the stretched quadrilateral mesh.
5. The evaluated relative computational error in calculation of interface position and temperature field, for all the studied cases, agreed with the analytical solution in 5 % and 10 % error band, respectively.

<figure><img src="image_15.png" /><figcaption>Fig. 13 Distribution of (a) volume fraction of vapor and (b) vapor temperature for the hybrid mesh.</figcaption></figure>

Although the proposed method was validated only for Stefan problem, the application of the method to engineering problems is considered to be straightforward because the method only uses a local value: the cell volume. Hence, the problem of ship resistance reduction using boiling flow can now be studied in the field of naval hydrodynamics CFD modeling.

## ACKNOWLEDGEMENT

The help received from Mr. Nithin Moran Narayan (TU Bremen) and Dr. Amine Ben Hadj Ali (ANSYS, via Ansys Student Community forum) is gratefully acknowledged. The valuable suggestions from Dr. Mijail Febres Soria and Mr. Lubomir Bureš (both Paul Scherrer Institut) are also sincerely appreciated

## APPENDIX: The derivation of correlation for calculation of local Nusselt number

The expression for local, cell-size dependent, Nusselt number stems from the balance between thermal phase change and energy jump mass transfer models. It is thereby assumed that liquid is at saturation temperature, $T_l = T'$. Other, model related, assumptions are also incorporated.

Thus, in thermal phase change model it is also assumed that vapor and liquid phase's scaling factors are equal to unity, that is, $C_v = C_l = 1.0$ and that interfacial area is computed as magnitude of volume fraction field, Eq. (13). Hence, Eq. (19), reduces to:

$$ \dot{m}_{lv} = - \frac{h_v |\nabla \alpha_v| (T' - T_v)}{h_{fg}} \qquad (A 1) $$

where $h_{fg}$ denotes specific heat of vaporization.

As stated in the introduction, one of the goals of the proposed model is to ensure the occurrence of dominant mass transfer rate in the interface zone. Therefore, in order to define energy jump model, consider the phase change process in a three-cell stencil, as shown in Fig. 14. Let the interface is thereby located between positions $i$ and $i + 1$ at $x$-coordinate of right sided Cartesian coordinate system. Furthermore, it is assumed that the superheated vapor phase occupies the zone on the left side of the interface, while the remainder of the domain is filled with liquid at saturation temperature.

The temperature gradient at the right side of the central cell's center is calculated as [30]:
$$ \frac{dT}{dx} = \frac{T' - T_i}{h} \qquad (A 2) $$

where $h$ denotes the distance between the interface and the center of adjacent vapor cell, see Fig. 14.

<figure><img src="image_16.png" /><figcaption>Fig. 14 The mass transfer process in a three-cell stencil.</figcaption></figure>

The mass transfer rate calculated using energy jump model is defined as [30]:

$$ \dot{m}_{lv} = \frac{\lambda_l(\nabla T_l) \cdot \vec{n} - \lambda_v(\nabla T_v) \cdot \vec{n}}{h_{fg}} \cdot \frac{A_i}{V_c} \qquad (A 3) $$

where indexes $l$ and $v$ denote, respectively, liquid and vapor phase, $\nabla T$ is the temperature gradient, $\vec{n}$ is the interface normal directed from vapor to liquid phase, $A_i$ is the interface area and $V_c$ is the cell volume.

Since the liquid is at saturation temperature, $\nabla T_l = 0$, and Eq. (A 3) is reduced to:

$$ \dot{m}_{lv} = - \frac{\lambda_v(\nabla T_v) \cdot \vec{n}}{h_{fg}} \cdot \frac{A_i}{V_c} \qquad (A 4) $$

The consideration of the problem as one-dimensional allows for the following simplification:

$$ |\nabla \alpha_v| \approx \frac{A_i}{V_c} $$

Hence, the energy jump model reads:

$$ \dot{m}_{lv} = - \frac{\lambda_v \nabla T_v |\nabla \alpha_v|}{h_{fg}} \qquad (A 5) $$

Approximation of vapor phase's temperature gradient using Eq. (A 2) gives:
$$ \dot{m}_{lv} = -\frac{1}{h_{fg}} \lambda_v \frac{T'-T_i}{h} |\nabla \alpha_v| \qquad (A 6) $$

Equalization of the right-hand sides of Eqs. (A 1) and (A 6) gives:

$$ h_v(T' - T_v) = \lambda_v \frac{T'-T_i}{h} \qquad (A 7) $$

The heat transfer is present only in vapor phase and, therefore, the heat transfer coefficient is defined as:

$$ h_v = \frac{\lambda_v N u_v}{d_v} \qquad (A 8) $$

Further approximations refer to *i*-th cell temperature and distance *h*, and are given by:

$$ T_i \approx T_v $$

$$ h \approx \frac{\Delta x}{2} $$

Involvement of these approximations together with Eq. (A 8) in Eq. (A 7) yields the local, cell-size dependent, Nusselt number

$$ Nu_v = \frac{2 \cdot d_v}{\Delta x} \qquad (A 9) $$

By rearranging Eq. (A 9) we obtain:

$$ Nu_v = \frac{d_p}{\frac{\Delta x}{2}} $$

which involved in Eq. (A 8) gives:

$$ h_v = \frac{\lambda_v}{\frac{\Delta x}{2}} $$

where the diameter of dispersed phase vanished. By involving it in the mass transfer equation, Eq. (A 1), and after subsequent rearranging we have:

$$ \dot{m} h_{fg} = -\lambda_v \frac{T_v - T'}{\frac{\Delta x}{2}} |\nabla \alpha_v| \qquad (A 10) $$

Thus, the heat transfer coefficient is a fictitious quantity since it is used to express the approximate distance from the cell center of the control volume that contains the interface to
the interface itself and to involve the thermal conductivity of the vapor phase in the mass transfer model. Furthermore, its inclusion in the mass transfer model yields an approximation of the temperature gradient at the interface, that appears from energy balance at the interface, i.e., in the heat flux in the normal direction to the interface.

REFERENCES

[1] G.A. Pavlov, L. Yun, A. Bliault, S.-L. He, Air Lubrication and Air Cavity Analysis. In: Air Lubricated and Air Cavity Ships. Springer, New York (2020). https://doi.org/10.1007/978-1-0716-0425-0

[2] Y. Kodama, M. Hinatsu, T. Hori, H. Kawashima, H. Takeshi, M. Makino, M. Otsuna, Y. Sanada, Y. Murai, S. Ota, A full-scale air lubrication experiment using a large cement carrier for energy saving (result and analysis), in: Proc. Japan Soc. Naval Architects and Ocean Engineers Conference, 2008, pp. 163-166.

[3] S. Mizokami, C. Kawakita, Y. Kodan, S. Takano, S. Higasa, R. Shigenaga, Development of Air lubrication system and verification by the full scale ship test, Journal of the Japan Society of Naval Architects and Ocean Engineers, 12 (2010) 69-77. https://doi.org/10.2534/jjasnaoe.12.69

[4] K. Mohanarangam, S.C.P. Cheung, J.Y. Tu, L. Chen, Numerical simulation of micro-bubble drag reduction using population balance model, Ocean Engineering, 36 (11) (2009) 863-872.
https://doi.org/10.1016/j.oceaneng.2009.05.001

[5] S. Sindagi, R. Vijayakumar, B.K. Saxena, Parametric CFD investigation of ALS technique on reduction in drag of bulk carrier, Ships and Offshore Structures, 15:4 (2020) 417-430.
https://doi.org/10.1080/17445302.2019.1661617

[6] R. Ezeta, D. Bakhuis, S.G. Huisman, C. Sun, D. Lohse, Drag reduction in boiling Taylor-Couette turbulence, Journal of Fluid Mechanics, 881 (2019) 104-118. https://doi.org/10.1017/jfm.2019.758

[7] Y. Sato, B. Niceno, Pool boiling simulation using an interface tracking method: From nucleate boiling to film boiling regime through critical heat flux, International Journal of Heat and Mass Transfer, 125 (2018) 876-890. https://doi.org/10.1016/j.ijheatmasstransfer.2018.04.131

[8] D.-L. Sun, J.-L. Xu, L. Wang, Development of a vapor-liquid phase change model for volume-of-fluid method in FLUENT, International Communications in Heat and Mass Transfer, 39 (8) (2012) 1101-1106.
https://doi.org/10.1016/j.icheatmasstransfer.2012.07.020

[9] Z. Pan, J.A. Weibel, S.V. Garimella, A saturated-interface-volume phase change model for simulating flow boiling, International Journal of Heat and Mass Transfer, 93 (2016) 945-956.
https://doi.org/10.1016/j.ijheatmasstransfer.2015.10.044

[10] I. Perez-Raya, S.G. Kandlikar, Chapter three - Evaporation on a planar interface – Numerical simulation and theoretical analysis of heat and mass transport processes, in: E.M. Sparrow, J.P. Abraham, J.M. Gorman, T.F. Irvine, J.P. Hartnett (Eds.) Advances in Heat Transfer, Elsevier, 2016, pp. 125-190.
https://doi.org/10.1016/bs.aiht.2016.08.005

[11] W.H. Lee, Pressure iteration scheme for two-phase flow modeling, in: "Multiphase transport: Fundamentals, Reactor safety, Applications". (1980) 407-432.

[12] A.S. Rattner, S. Garimella, Simple mechanistically consistent formulation for volume-of-fluid based computations of condensing flows, Journal of heat transfer, 136 (7) (2014).
https://doi.org/10.1115/1.4026808

[13] D. Sun, J. Xu, Q. Chen, Modeling of the evaporation and condensation phase-change problems with FLUENT, Numerical Heat Transfer, Part B: Fundamentals, 66 (4) (2014) 326-342.
https://doi.org/10.1080/10407790.2014.915681

[14] D.G. Kim, C.H. Jeon, I.S. Park, Comparison of numerical phase-change models through Stefan vaporizing problem, International Communications in Heat and Mass Transfer, 87 (2017) 228-236.
https://doi.org/10.1016/j.icheatmasstransfer.2017.07.013

[15] I. Perez-Raya, S.G. Kandlikar, Discretization and implementation of a sharp interface model for interfacial heat and mass transfer during bubble growth, International Journal of Heat and Mass Transfer, 116 (2018) 30-49. https://doi.org/10.1016/j.ijheatmasstransfer.2017.08.106

[16] G. Chen, T. Nie, X. Yan, An explicit expression of the empirical factor in a widely used phase change model, International Journal of Heat and Mass Transfer, 150 (2020) 119279.
https://doi.org/10.1016/j.ijheatmasstransfer.2019.119279
[17] C.R. Kharangate, I. Mudawar, Review of computational studies on boiling and condensation, International Journal of Heat and Mass Transfer, 108 (2017) 1164-1196.
https://doi.org/10.1016/j.ijheatmasstransfer.2016.12.065

[18] S. Mer, O. Praud, H. Neau, N. Merigoux, J. Magnaudet, V. Roig, The emptying of a bottle as a test case for assessing interfacial momentum exchange models for Euler–Euler simulations of multi-scale gas-liquid flows, International Journal of Multiphase Flow, 106 (2018) 109-124.
https://doi.org/10.1016/j.ijmultiphaseflow.2018.05.002

[19] S.n. Fleau, Modélisation des écoulements eau-vapeur « tous régimes d'écoulements » par une approche multi-champ, 2017.

[20] P. Coste, A Large Interface Model for two-phase CFD, Nuclear Engineering and Design, 255 (2013) 38-50.
https://doi.org/10.1016/j.nucengdes.2012.10.008

[21] R. Denèfle, S. Mimouni, J.-P. Caltagirone, S. Vincent, Multifield hybrid approach for two-phase flow modeling – Part 1: Adiabatic flows, Computers & Fluids, 113 (2015) 106-111.
https://doi.org/10.1016/j.compfluid.2014.07.018

[22] D.A. Drew, S.L. Passman, Theory of Multicomponent Fluids, Springer, 1999.
https://doi.org/10.1007/b97678

[23] ANSYS Inc., ANSYS Fluent Theory Guide, 2019.

[24] T. Höhne, C. Vallée, Experiments and numerical simulations of horizontal two-phase flow regimes using an interfacial area density model, Journal of Computational Multiphase Flows, 2 (3) (2010) 131-143.
https://doi.org/10.1260/1757-482X.2.3.131

[25] L. Štrubelj, I. Tiselj, B. Mavko, Simulations of free surface flows with implementation of surface tension and interface sharpening in the two-fluid model, International Journal of Heat and Fluid Flow, 30 (4) (2009) 741-750.
https://doi.org/10.1016/j.ijheatfluidflow.2009.02.009

[26] ANSYS Inc., VOF Model, in, 2017.

[27] V.P. Carey, Liquid-vapor phase-change phenomena, Second edition, Taylor & Francis, 2007.

[28] C. Geuzaine, J.F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, 79 (11) (2009) 1309-1331.
https://doi.org/10.1002/nme.2579

[29] F. Gauss, D. Lucas, E. Krepper, Grid studies for the simulation of resolved structures in an Eulerian two-fluid framework, Nuclear Engineering and Design, 305 (2016) 371-377.
https://doi.org/10.1016/j.nucengdes.2016.06.009

[30] Y. Sato, B. Ničeno, A sharp-interface phase change model for a mass-conservative interface tracking method, Journal of Computational Physics, 249 (2013) 127-161.
https://doi.org/10.1016/j.jcp.2013.04.035