Received: 28 November 2016 Revised: 15 January 2017 Accepted: 15 January 2017

DOI: 10.1002/htj.21268

RESEARCH ARTICLE

The evaluation of the diffuse interface method

for phase change simulations using OpenFOAM
Nima Samkhaniani Mohamad Reza Ansari

Faculty of Mechanical Engineering, Tarbiat

Modares University, Tehran, I.R. of Iran Abstract
Correspondence In the present study, a new solver named phaseChangeHeat-
Nima Samkhaniani, Faculty of Mechanical
Engineering, Tarbiat Modares University,
Foam is implemented on the OpenFOAM cfd package to
P.O. Box 14117 13116, Tehran, I.R. of Iran simulate boiling and condensation. The solver is capturing
Email: [Link]@[Link], the interface between two immiscible phases with a color
[Link]@[Link]
function volume of fluid (CF-VOF) method. The two flu-
ids (vapor and liquid) are assumed Newtonian and incom-
pressible. The surface tension is modeled with continuous
surface force (CSF) which is improved with a Lafaurie fil-
ter to suppress the spurious current. The mass flux across
the interface in the phase change process is determined
by either Lee or Tanasawa mass transfer models. Addi-
tionally, the slight variation of saturation temperature with
local pressure is considered with the simplified Clausius–
Clapeyron relation. The coupled velocity pressure equation
is solved using the PIMPLE algorithm. The new solver is
validated and examined with (i) Stefan problem, (ii) two-
dimensional film boiling, (iii) the film condensation on a
horizontal plate, (iv) the laminar film condensation over a
vertical plate, and (v) bubble condensation in subcooled

Nomenclature: 𝐶𝛼 , Compression factor (−); 𝐶𝑝 , Specific heat (J/kgK); 𝑑𝐿 , Liquid thermal diffusivity: 𝑘∕𝜌𝐶 (m2 /s); 𝐷, Equiv-
𝑝
alent diameter (M); 𝐷0 , Bubble Initial diameter (M); 𝑔,⃗ gravity acceleration (m/s2 ); 𝐻𝐿𝐺 , Latent heat (J/kg); 𝑘, Thermal con-
ductivity (W/mK); 𝑀, Molar mass (kg/Kmol); ṁ , mass flow rate per unit area (kg/m2 s); ṁ ′′′ , mass flow rate per unit volume
′′

(kg/m3 s); 𝑃 , Pressure (Pa); 𝑃𝑟𝑔ℎ , Dynamic Pressure (Pa); 𝑟, Mass transfer intensity factor (1/s); 𝑅, Universal gas constant (J/mol
K); 𝑇 , Temperature (K); 𝑊 , width (M); 𝑈⃗ , velocity (m/s); 𝑈⃗𝑐 , Compressive velocity (m/s)
Greek symbols: 𝛼, Volume fraction factor; 𝜅, Interface curvature (m−1 ); 𝜇, Dynamic viscosity (Pa.s); 𝜆, Taylor unstable wave
(M); 𝜌, Density (kg/m3 ); 𝛾, Tanasawa Coefficient (−); 𝜎, Surface Tension (N/m)
Subscripts: 𝑐, Condensation; 𝑒, Evaporation; 𝐺, Gas (vapor) phase; 𝐿, Liquid phase; sat, Saturation; sub, Subcooled; sup,
Superheat; W, wall

Heat Transfer—Asian Res 2017; 00: 1–31 [Link]/journal/htj © 2017 Wiley Periodicals, Inc. 1
2 SAMKHANIANI AND ANSARI

boiling. The present study shows the capability of a dif-

fuse interface method in accurate simulation of the phase
change process and it is expected to be instructive for
further numerical studies in this area.

KEYWORDS
boiling, condensation, interFoam, OpenFOAM, volume of fluid

1 I N T RO D U C T I O N

1.1 Importance of predictive tools for phase change phenomena

The diabatic flows have significant industrial importance in heating and cooling systems. Two-phase
flow heating and cooling systems have gained significant popularity in recent years due to an urgent
need to overcome heat dissipation challenges in aerospace and commercial electronic devices. Heat
dissipation rates at such devices are no longer manageable with single phase cooling system. Therefore,
more aggressive two-phase solutions are being sought. Recent studies concerning the development of
two-phase cooling solutions have focused on heat acquisition or rejection by boiling and condensation,
respectively. While empirical formulations constitute the most popular means for predicting transport
behavior in designing cooling system, these formulations are limited by working fluids, geometries,
and operating conditions of databases upon which they are based. In order to achieve more universal
predictive tools, there is a desire to develop computational techniques.

1.2 Application of numerical models in phase change system

Several types of approaches such as Lagrangian, Eulerian, and Eulerian–Lagrangian have been devel-
oped for simulation of two-phase systems in recent decades.
The Lagrangian or Eulerian–Lagrangian approaches such as smoothed-Particle Hydrodynamic
(SPH)1 and front tracking (FT)2 are limited to simple cases due to the complexity of applying compli-
cated grid topologies. Hence, Eulerian methods are more popular due to their simplicity and feasibility
of tackling multiple interfaces.
The two main representative examples of Eulerian approach are Level-Set (LS)3 and volume of fluid
(VOF).4 There are also several Hybrid methods5 such as CLSVOF or VOSET which are couple of level
set and volume of fluid approaches.6–10 The numerical method in this approach is classified based on
the thickness of interface. An interface is the thin boundary layer that separates two distinct phases
of matter. The interface thickness in the numerical method as shown in Fig. 1 is zero (sharp) or finite
(diffusive). In sharp interface methods, the physical interface is a functional interface of zero thickness
and physical quantities such as density and viscosity are discontinuous at the interface. Volume of fluid
method with interface reconstruction such as PLIC11 and PROST12 are two examples of the sharp
interface method. In diffuse interface methods such as CICSAM13 and THINC,14 the interface has a
finite thickness and physical quantities vary continuously across the interface. However, the numerical
interface thickness in the diffuse method is much larger than the actual physical thickness which is for a
liquid–fluid interface typically a few nanometers.15 The main disadvantage of interface reconstruction
volume of fluid (IR-VOF) is the complexity of the interface reconstruction, especially in 3D simulations
SAMKHANIANI AND ANSARI 3

FIGURE 1 Comparison of sharp and diffuse interface methods

and physics involving extensive topology changes. Therefore, volume of fluid method without interface
reconstruction which called color function volume of fluid has gained more popularity.
The Eulerian approach has been developed to simulate flow with phase change in the last decades.
Welch and Wilson16 introduced a PLIC-VOF approach for simulation of boiling. Hardt and Wondra17
proposed an evaporation model that relies on a continuum field representation of the source terms in
the mass conservation equation and applied on PLIC-VOF in Fluent software. In order to simulate
evaporating droplets, Schlottke and Weigand18 developed a PLIC-VOF method which involves two
VOF variables, one for the Liquid phase and one for the vapor phase. Son and Dhir19 and Son et al.20
modified the Level set method to account for phase change. They assumed that the vapor phase is at
saturation temperature, so they solved the energy equation only for the liquid phase. Level set methods
for incompressible flow with phase change are also presented by Tanguy et al.21 and Gibou et al.22 who
applied the ghost fluid approach to impose the jump conditions at the interface. A CLSVOF method for
simulation of bubble growth in film boiling was developed by Tomar et al.23 Guo et al.64 simulated film
boiling using VOSET method. Kunkelmann and Stephan24 simulated nucleating boiling with special
treatment for heat flux in contact angle.
There are several Eulerian codes to simulate multiphase flows, a color function volume of fluid
(CF-VOF) solver named interFoam in open source CFD packge OpenFOAM25 has been receiving
an increased amount of attention and usage.26–28 This VOF solver was implemented by Ubbink29
at first, then it has undergone several modifications. The present code is a part of C++ libraries of
OpenFOAM. This software was designed for finite volume discretization of a generic set of partial
differential equations. However, several useful features such as ease of parallelization, availability
of preprocessing and post processing utilities, error checks and choice of spatial and temporal dis-
cretization schemes make the code attractive for research. InterFoam has been recently contributed to
consider phase change phenomena (boiling and condensation). Table 1 provides a comprehensive sum-
mery of literature concerning the implementation of different computational phase change model using
OpenFOAM.
4 SAMKHANIANI AND ANSARI

TABLE 1 Summary of computational phase change simulation in OpenFOAM frame work

Author Applications Model Remarks
Kunkelmann and Nucleate boiling Tanasawa axisymmetric grid, 3-dimensional grid with
Stephan30 adaptive grid refinement, contact line model for
evaporation, VOF
Kunkelmann and Nucleate boiling Sharp interface axisymmetric grid, couple of level set and VOF,
Stephan24 model
Zeng et al.31 Bubble condensation Tanasawa 3-dimensional grid, simple couple of level set and
VOF
Doro32 Falling film evaporation Tanasawa 2-dimensional grid, laminar, VOF
Rattner and Film condensation – 2-dimensional grid, laminar, VOF, new phase
Garimella33 change model
Bahreini et al.34 Bubble condensation Lee 2-dimensional grid, laminar, VOF

1.3 Phase change models

1.3.1 Tanasawa model
Scharge35 proposed a phase change model based on the Hertz–Knudsen equation,36 with the interfacial
jump in the temperature and pressure 𝑇𝑠𝑎𝑡 (𝑃𝐿 ) = 𝑇𝐿 ≠ 𝑇𝑠𝑎𝑡 (𝑃𝐺 ) = 𝑇𝐺 . Mass flux at the interface is
given by

√ [ ]
′′ 2 𝑀 𝑃𝐺 𝑃𝐿
𝑚̇ = 𝛾 √ − 𝛾𝑒 √ , (1)
2 − 𝛾𝑐 2𝜋𝑅 𝑐 𝑇𝐺 𝑇𝐿

where 𝑅 = 8.314 J/mol K is the universal gas constant, 𝛾 the fraction of molecules transferred from one
phase to the other during phase change. The subscripts 𝑐 and 𝑒 in the equation refer to condensation
and evaporation, respectively. Here𝛾𝑒 = 1 and 𝛾𝑐 represent complete evaporation and complete con-
densation, respectively. In phase change simulation, usually equal values of 𝛾𝑐 and 𝛾𝑒 are considered.
While 𝛾 = 0.1 − 1 for dynamically renewing water surfaces such as jets and moving films and 𝛾 < 0.1
for stagnant surfaces are recommended.37 Tanasawa38 simplified Eq. (1) using extra assumptions. He
assumed that the interface is at saturation temperature and the heat flux linearly dependent on the tem-
perature jump between the interface and the vapor. The mass flux rate in Tanasawa’s model is given by

√
′′2𝛾 𝑀 𝜌𝐺 𝐻𝐿𝐺 (𝑇 − 𝑇𝑠𝑎𝑡 )
𝑚̇ = , (2)
2−𝛾 2𝜋𝑅 𝑇
3∕2
𝑠𝑎𝑡

where 𝑇𝑠𝑎𝑡 is based on local pressure. The volumetric transferred mass rate in (kg/m3 s) is defined as

′′′ ′′
𝑚̇ = 𝑚̇ |∇𝛼𝐿 |. (3)
SAMKHANIANI AND ANSARI 5

1.3.2 Lee model

Lee39 assumed that mass is transferred at a constant pressure in phase change flow and derived a model
for a quasi-thermo-equilibrium state. The volumetric transferred mass rate is given by

′′′ 𝑇 − 𝑇𝑠𝑎𝑡
𝑚̇ = 𝑟𝑐 (1 − 𝛼𝐿 )𝜌𝐺 , for condensation 𝑇 < 𝑇𝑠𝑎𝑡 (4)
𝑇𝑠𝑎𝑡

and

′′′ 𝑇 − 𝑇𝑠𝑎𝑡
𝑚̇ = 𝑟𝑒 𝛼𝐿 𝜌𝐿 , for evaporation 𝑇 > 𝑇𝑠𝑎𝑡 , (5)
𝑇𝑠𝑎𝑡

where 𝑟𝑐 and 𝑟𝑒 in 𝑠−1 are empirical coefficients called the mass transfer intensity factor. The extensive
ranges from 0.1 to 107 have been conducted in literature for 𝑟.40–45

1.3.3 Sharp interface model

The sharp interface model is based on the Rankine–Hungoniot jump condition for energy conservation
where all the heat transferred at the interface is assumed to be consumed by latent heat due to phase
change and given by

′′′ 𝑘𝑒𝑓 𝑓 (∇𝛼𝐿 .∇𝑇 )

𝑚̇ = , (6)
𝐻𝐿𝐺

where 𝑘𝑒𝑓 𝑓 is the effective thermal conductivity, which can be determined from the volume fractions
and thermal conductivities of the vapor and liquid.46

1.4 Objectives of study

The sharp interface method is the preference in the implementation of the phase change process, the
capability of the diffuse interface method in correct simulation of phase change phenomena is not
well defined. In addition, despite several implementations of the phase change model (boiling or con-
densation) in the framework of OpenFOAM, there is no exclusive and detailed documented resource,
and proposed numerical methods are often validated for user specific interest. Therefore, in the present
study a solver called “phaseChangeHeatFoam” is developed to simulate boiling and condensation. This
article provides detailed information on governing equation, code structure, and numerical method for
implementing the phase change process in the framework of OpenFOAM. Moreover, this study pro-
vides several problems to investigate the performance of the present implementation.

2. COMPUTAT I O NA L ME T H O D

2.1 Governing equations

In this study, both vapor liquid phases are assumed incompressible and immiscible. The volume of
fluid method uses a discontinuous scalar color function to resolve the interface among two phases in
6 SAMKHANIANI AND ANSARI

fixed grids. This scalar function is the ratio of one fluid volume to the volume of cell, and it is defined
as

⎧1 𝑥⃖⃗ ∈ Liquid
𝑉𝐿𝑖𝑞𝑢𝑖𝑑 ⎪
𝛼𝐿 (⃖𝑥,
⃗ t) = = ⎨0 < 𝛼𝐿 < 1 𝑥⃖⃗ ∈ interface . (7)
𝑉 ⎪
⎩0 𝑥⃖⃗ ∈ vapor

The thermo-physical properties of two-phase flow such as viscosity (𝜇), density (𝜌), and thermal
conductivity (𝑘) are defined as

𝑦 = 𝛼𝐿 𝑦𝐿 + (1.0 − 𝛼𝐿 )𝑦𝐺 , 𝑦 ∈ [𝜇, 𝜌, 𝑘], (8)

where 𝑦 has liquid physical properties where 𝛼𝐿 = 1.0, and it has vapor physical properties where
𝛼𝐿 = 0.0. The interface is moving along the flow; therefore the volume fraction is advected to keep
the interface. In order to derive the transport equation for volume fraction, density (𝜌) from Eq. (8) is
placed in Eq. (9) and rearranged to obtain Eq. (10). Eq. (9) indicates the global continuity equation

𝜕 ⃖⃖⃗) = 0,
(𝜌) + ∇.(𝜌𝑈 (9)
𝜕𝑡

𝜕𝛼𝐿 𝜌 ∇.𝑈⃖⃖⃗
+𝑈 ⃖⃖⃗ = − 𝐺
⃖⃖⃗.∇𝛼𝐿 + 𝛼𝐿 ∇.𝑈 , (10)
𝜕𝑡 (𝜌𝐿 − 𝜌𝐺 )

where ∇.𝑈⃖⃖⃗ is zero in an incompressible adiabatic two-phase flow. In the phase change process, it is
determined from the local continuity for each phase:

𝜕 ⃖⃖⃗) = −𝑚̇ ′′′ ,

(𝜌 ) + ∇.(𝜌𝐿 𝑈 (11)
𝜕𝑡 𝐿

𝜕 ⃖⃖⃗) = +𝑚̇ ′′′ .

(𝜌𝐺 ) + ∇.(𝜌𝐺 𝑈 (12)
𝜕𝑡
As both phases are incompressible, the first term on the left hand of Eqs. (11) and (12) is zero.
Then, by using the integral approach and divergence theorem, Eqs. (11) and (12) are converted to the
following equations:

⃖⃖⃗.⃖⃖⃖ 𝑚̇
𝑈 𝑛𝐿⃗𝑑𝑆𝐿 = − , (13)
∫ 𝜌𝐿
𝐿

⃖⃖⃗.⃖⃖⃖⃖ 𝑚̇
𝑈 𝑛𝐺⃗𝑑𝑆𝐺 = + . (14)
∫ 𝜌𝐺
𝐺

The unit normal vector at the interface is displayed in Fig. 2. As 𝑛⃖⃗𝐺 = 𝑛⃖⃗𝐼 = −⃖𝑛⃗𝐿 , then the summation
of the former equations results in
( )
⃖⃖⃗.⃖𝑛⃗𝑑𝑆 = 𝑚̇ 1 1
𝑈 − . (15)
∫ 𝜌𝐺 𝜌𝐿
𝑆
SAMKHANIANI AND ANSARI 7

FIGURE 2 Two-phase flow system, S denotes surface and n denotes unit normal, L: liquid, G: vapor, and
I: interface

The differential form of Eq. (15) is

( )
⃖⃖⃗ = 𝑚̇ ′′′ 1 1
∇.𝑈 − . (16)
𝜌𝐺 𝜌𝐿

By replacing Eq. (16) into Eq. (10), the following transport equation for 𝛼𝐿 is derived:
[ ( )]
𝜕𝛼𝐿 1 1 1
⃖⃖
⃗
+ 𝑈 .∇𝛼𝐿 = −𝑚̇
′′′
− 𝛼𝐿 − . (17)
𝜕𝑡 𝜌𝐿 𝜌𝐿 𝜌𝐺

In OpenFOAM, the extra divergence term is added to Eq. (17). This term contributes only in the
region of the interface (0 < 𝛼𝐿 < 1.0). It limits the smearing of the interface because of the compen-
sation of the diffusive fluxes28
[ ( )]
𝜕𝛼𝐿 1 1 1
⃖⃖
⃗ ⃖⃖⃖⃗
+ 𝑈 .∇𝛼𝐿 + ∇.(𝛼𝐿 (1 − 𝛼𝐿 )𝑈𝑐 ) = −𝑚̇
′′′
− 𝛼𝐿 − . (18)
𝜕𝑡 𝜌𝐿 𝜌𝐿 𝜌𝐺

⃖⃖⃖⃗𝑐 is compressive velocity. It is calculated in the normal direction to the interface to avoid any
Here 𝑈
dispersion.47 It is defined as
∇𝛼𝐿
⃖⃖⃖⃗𝑐 = min{𝐶𝛼 |𝑈 |, max(|𝑈 |)}
𝑈 . (19)
|∇𝛼𝐿 |

The coefficient 𝐶𝛼 controls the weight of the compression flux and should usually be in the range
of unity (1.0 < 𝐶𝛼 < 4.0).28,48,49 In the present study, the compression factor 𝐶𝛼 = 1.0 is applied. The
large compression factor increases the non-physical spurious current in the vicinity of the interface.
The momentum equation is given by
⃖⃖⃗)
𝜕(𝜌𝑈 ⃖⃖⃗𝑈
+ ∇.(𝜌𝑈 ⃖⃖⃗) − ∇.(𝜇(∇𝑈 ⃖⃖⃗)) = −∇𝑃 + 𝜌⃖𝑔⃗ + 𝜎𝜅∇𝛼𝐿 .
⃖⃖⃗𝑇 + ∇𝑈 (20)
𝜕𝑡
The last term in the right hand of Eq. (20) is considering surface tension using the continuous surface
force (CSF) model.50 Where 𝜎 is the surface tension, and 𝜅 is the interface curvature and is defined as

( )
𝛼𝐿
∇̃
𝜅 = −∇. , (21)
|∇̃
𝛼𝐿 |
8 SAMKHANIANI AND ANSARI

where 𝛼̃𝐿 is a smoothed VOF function calculated from 𝛼𝐿 by smoothing it over a finite region around the
interface. In the VOF method, the fluid interface sharply changes over a thin region. This abrupt change
of the VOF function creates errors in calculating the normal vectors and the curvature of the interface,
which are used to evaluate the interfacial forces. These errors induce spurious currents in the interfacial
region.28,51 The spurious currents create extra heat convection around the interface which increases
local mass transfer. An easy way to suppress these artifacts is to compute the interface curvature from
a smoothed VOF function. In the present study, the smoother proposed by Lafaurie et al.52 is applied:
∑𝑛
𝑓 =1 𝛼𝐿𝑓 𝑆𝑓
𝛼
̃𝐿𝑝 = ∑𝑛 , (22)
𝑓 =1 𝑆𝑓

where 𝑆𝑓 is the magnitude of face area, the subscript 𝑝 denotes the cell index, and 𝑓 denotes the
face index. The interpolated value (𝛼𝐿𝑓 ) at the face center is calculated using linear interpolation.
The application of this filter can be repeated 𝑚 times to get a smoothed field. It should be noted that
smoothing tends to wipe out high curvature regions and should be applied only up to the level that is
strictly necessary to sufficiently suppress parasitic currents. Hoang et al.53 found that the magnitude
of parasitic current decreases up to one order from 𝑚 = 0 to 𝑚 = 2 and only a slight further decrease
was observed for 𝑚 > 2. Therefore, in present study 𝑚 = 2 is employed in all simulations.
Energy equation is given by

𝜕 ′′′
(𝜌𝐶𝑝 𝑇 ) + ∇.(𝜌𝐶𝑝 𝑈 𝑇 ) − ∇.(𝑘∇𝑇 ) = −𝑚̇ 𝐻𝐿𝐺 , (23)
𝜕𝑡
where 𝐻𝐿𝐺 is latent heat coefficient. The last term in Eq. (23) is added to consider latent heat during
phase change.
The interface temperature is often at saturation, it is almost constant during phase change. However,
the slight variation in the quantity of saturation temperature as a function of local pressure is determined
using a simplified version of the Clausius–Clapeyron relation which is based on ideal gas assumption
for the vapor phase
( )
𝑃𝑠𝑎𝑡,1 𝑀𝐻𝐿𝐺 1 1
ln =− − . (24)
𝑃𝑠𝑎𝑡,0 𝑅 𝑇𝑠𝑎𝑡,1 𝑇𝑠𝑎𝑡,0

It is well known that numerical interfacial flow computations become more challenging when the
imbalance of material properties between the phases increases.51 In order to reduce the material prop-
erties imbalance across the interface, and increase numerical robustness, Eq. (23) is rearranged as
𝜕 ′′′
(𝑇 ) + ∇.(𝑈 𝑇 ) − ∇.(𝐷𝑘 ∇𝑇 ) = −𝐷𝑐 𝑚̇ 𝐻𝐿𝐺 , (25)
𝜕𝑡
where 𝐷𝑘 and 𝐷𝑐 are defined as
1.0
𝐷𝑐 = , (26)
𝜌𝐿 𝐶𝐿 𝛼𝐿 + 𝜌𝐺 𝐶𝐺 (1.0 − 𝛼𝐿 )

𝑘𝐿 𝛼𝐿 + 𝑘𝐺 (1.0 − 𝛼𝐿 )
𝐷𝑘 = . (27)
𝜌𝐿 𝐶𝐿 𝛼𝐿 + 𝜌𝐺 𝐶𝐺 (1.0 − 𝛼𝐿 )
To close the above set of equations, an appropriate phase change model is necessary to calculate
′′′
transferred mass flux rate (𝑚̇ ). In the present study, Lee and Tanasawa mass transfer model are
employed.
SAMKHANIANI AND ANSARI 9

2.2 Numerical details

In the present study, OpenFOAM which is an open source software is employed. It is a flexible
and efficient C++ library for the customization and extension of applications (solvers and utili-
ties) and operating on scalar, vectorial, and tensorial fields.25 It is based on an unstructured mesh
formulation with a collocated cell-centered variable arrangement and provides numerical methods
for discretization of partial differential equations along with solvers for the corresponding numer-
ical solution of the resulting system of algebraic equations. There are two main tensor derivative
namespaces: the 𝑓 𝑣𝑚 (finiteVolumeMethod) namespace and the 𝑓 𝑣𝑐 (finiteVolumeCalculus) names-
𝜕
pace to consider differential operators such as ∇., ∇2 , 𝜕𝑡 . These differential operators correspond to
fvm::div(), fvm::laplacian()and fvm::ddt() or alternatively fvc::div(), fvc::laplacian()and fvc::ddt() in
OpenFOAM. The fvc functions perform an explicit evaluation of predetermined data and they map
from one field to another. In contrast, fvm functions construct appropriate matrices using the finite vol-
ume discretization. They enable to create entire matrix representations of differential equations using
fvScalarMatrix or fvVectorMatrix and their implicit numerical solution. These equation objects hold
the matrices that represent the equations and handle the numerical solution. Moreover, there are geo-
metric tensor field classes for example: volScalarField and volVectorField that contain field data on cell
centers to represent field’s variables such as velocity and pressure. The tensor field classes also have a
reference to the mesh and include boundary information, previous time steps data which are necessary
for temporal descritization, and SI dimension set information. Additionally, surfaceScalarField and
surfaceVectorField are available which hold data at the center of surfaces of control volumes.
The present solver phaseChangeHeatFoam is developed on interFoam solver with OpenFOAM220.
InterFoam is a solver for two incompressible, isothermal immiscible fluids using a volume of fluid
phase-fraction based interface capturing approach. The momentum and other fluid properties are of the
mixture and a single momentum equation is solved. In order to simulate boiling and condensation, the
thermal energy transport equation (TEqn.H) is added to interFoam. Additionally, the volume fraction
transport equation (alphaEqn.H) and momentum equations (pEqn.H) are modified with appropriate
source terms. It provides support for generic phase-change models. Thus, phase-change models can be
selected by the user during runtime, and new models can be rapidly incorporated. A summary of the
phaseChangeHeatFoam algorithm is presented in Fig. 3.
In the beginning of the simulation, the solver loads the mesh and reads tensor fields (𝛼𝐿 , 𝑇 , 𝑈⃖⃖⃗ and
𝑃𝑟𝑔ℎ ) and boundary conditions. 𝑃𝑟𝑔ℎ = 𝑃 − 𝜌𝑔ℎ is dynamic pressure and ℎ is the fluid height. 𝑃𝑟𝑔ℎ is
applied to avoid any sudden changes in the pressure at the boundaries for hydrostatic terms.
The main solver loop is then initiated. First, the time step is dynamically modified to ensure numer-
ical stability. The adaptive time step is adjusted to satisfy the CFL condition and user-defined time
step (maxDeltaT).27 User-defined time step is employed to apply manually other temporal scales in the
problem for instance the thermal diffusion time scale.33
Next the discretized volume fraction Eq. (18) is solved for a user defined number of sub-time steps.
Four sub-cycles are considered here. The volume fraction equation (alphaEqn.H) is modified as shown
in Fig. 4. The volume fraction advection equation is solved using the multidimensional universal limiter
with explicit solution (MULES) method which is based on the method of flux corrected transport
(FCT) where an additional limiter is used to cutoff the face-fluxes at the critical values.54 This solver
is included in the OpenFOAM library, and performs the conservative transport equation of hyperbolic
transport equations with defined bounds (0 and 1 for 𝛼𝐿 ).
Once the update phase field is obtained, the program updates thermo physical properties according
to Eqs. (8), (26) and (27). Then the program enters the pressure-velocity loop. The process of correcting
the pressure and velocity fields in sequence is known as pressure implicit with splitting of operators
10 SAMKHANIANI AND ANSARI

FIGURE 3 PhaseChangeHeatFoam algorithm summary

FIGURE 4 Modified alphaEqn.H in phaseChangeHeatFoam

SAMKHANIANI AND ANSARI 11

FIGURE 5 Modified pEqn.H in phaseChangeHeatFoam

FIGURE 6 TEqn.H in phaseChangeHeatFoam

(PISO).55 In the OpenFOAM frame work, PISO can be repeated for multiple iterations at each time
step. This process is referred to as the PIMPLE loop. Here, only one outer correction is applied and the
PIMPLE loop works in PISO mode. First, the matrix equation for the momentum equation is formed.
Then the inner pressure–velocity correction process is initiated. In PISO, an intermediate velocity
field is first obtained, and the cell-face volume fluxes (𝜙) are evaluated and corrected for gravitational
forces, the continuum surface-tension force, and boundary conditions. The pressure-Poisson equation
is then formed and solved. Following the approach of Ref. 25 the coefficients of the pressure equation
are obtained from the diagonal entries of the momentum matrix equation (𝐴𝐷 ). For condensing and
evaporating flow, the pressure equation would be

( ) ( )
1 ′′′ 1 1
∇. ∇𝑝 = ∇.𝜙 − 𝑚̇ − . (28)
𝐴𝐷 𝜌𝐺 𝜌𝐿

The last term on the left hand of Eq. (28) is considering the phase change process. The modified pres-
sure equation (pEqn.H) is shown in Fig. 5. In order to avoid a checker-boarding effect in the momentum
equation, the Rhie–Chow momentum interpolation56 is applied.
Finally, the thermal energy transport equation (Eq. (25)) is solved. The TEqn.H is shown in Fig. 6.
For numerical stability, source terms in TEqn.H and pEqn.H are linearized.
The terms vDotAlpha, vDotP and vDotT are mass flux rates which are calculated based on selected
mass transfer models. Lee mass transform model implemented in Lee.C is shown in Fig. 7.
The discretization scheme is selectable by the user; in this series of simulations the schemes dis-
played in Table 2 are chosen. For convenience, the corresponding terminology of OpenFOAM is given.
12 SAMKHANIANI AND ANSARI

FIGURE 7 Lee.C in phaseChangeHeatFoam

3. VA L I DAT I O N A N D V E R I F I CAT I O N

In this section, for validation of the present solver, five problems are considered. Among them, the
Stefan problem, horizontal film condensation on a horizontal plate, and laminar film condensation on
a vertical plate have analytical solutions. For the others including 2D film boiling and bubble conden-
sations, the CFD result is compared with experimental data and correlations or previous numerical
simulations using the sharp interface method. Since there is no interface reconstruction in the present
VOF method, an iso-contour 𝛼𝐿 = 0.5 is assumed as the interface in all following simulations.

3.1 Stefan problem

The one-dimensional Stefan problem is a well-known benchmark for boiling flow.16,17,60,61 A vapor
film separates saturated liquid from a super-heated wall. Liquid and vapor are initially in quiescent
SAMKHANIANI AND ANSARI 13

TABLE 2 Discretization and interpolation schemes of the numerical model

Term Discretization scheme Method
𝜕 ⃖⃖⃗),
(𝜌𝑈 𝜕 ⃖⃖⃗𝑇 )
(𝜌𝑈 Euler The first-order bounded implicit scheme
𝜕𝑡 𝜕𝑡
∇.(𝜌𝑈⃖⃖⃗𝑈
⃖⃖⃗) vanLeerV Similar to VanLeer scheme57 modified for vector field
⃖⃖⃗𝛼𝐿 ), ∇.(𝜌𝑈
∇.(𝑈 ⃖⃖⃗𝑇 ) vanLeer See Ref. 57
⃖⃖⃖⃗𝑐 𝛼𝐿 (1 − 𝛼𝐿 ))
∇.(𝑈 InterfaceCompression See Ref. 58
∇𝜒 ∗ Linear Central difference schemes
∗∗
∇ 𝑓1 𝜒 corrected Surface normal gradient with correction on non-orthogonal
meshes59
∇.(𝜒1 ∇𝜒2 ) Linear corrected Face values (𝜒1 ) approximated by central difference scheme,
and the resulting surface normal gradient is calculated using
central difference schemes with non-orthogonal correction
Term Interpolation scheme Method
𝜒𝑓 Linear Default interpolation schemes for getting face values from cell
values
𝜒 : volScalarField or volVectorField
𝜒𝑓 : surfaceScalarField or surfaceVectorField

equilibrium. Evaporation pushes liquid away from the heated wall. The analytical solution for interface
position 𝛿(𝑡) and temperature 𝑇 (𝑥, 𝑡) on the vapor side are given by16
√
𝛿𝑎𝑛 (𝑡) = 2𝜂 𝑑𝐺 𝑡, (29)

( )
𝑇 − 𝑇𝑤 𝑥
𝑇 (𝑥, 𝑡) = 𝑇𝑤 + 𝑠𝑎𝑡 𝑒𝑟𝑓 √ , (30)
𝑒𝑟𝑓 (𝜂) 2 𝑑𝐺 𝑡

where 𝑑𝐺 is vapor thermal diffusivity, and 𝜂 is determined from

𝑐𝑝𝐺 (𝑇𝑤 − 𝑇𝑠𝑎𝑡 )

𝜂 exp(𝜂 2 )𝑒𝑟𝑓 (𝜂) = √ . (31)
𝜋𝐻𝐿𝐺

A quasi 1D computational domain with only one grid cell in the direction of translational invariance
is considered. A very thin vapor film is inserted in the computational domain in the initial time near the
hot wall. The liquid-vapor thermo physical properties for water at saturation pressure 1.0 Mpa are cho-
𝜌
sen which correspond to the moderate density ratio of ( 𝜌𝐿 ≈ 170). In order to ensure that the coefficient
𝐺
of mass flux rate in the energy equation is constant in the CFD model during the phase change process,
liquid and vapor phases specific heat are assumed equal (𝐶𝑝𝐺 = 𝐶𝑝𝐿 (𝑃𝑠𝑎𝑡 )). A schematic of Stefan
problem is illustrated in Fig. 8. A no slip boundary condition is employed for the velocity boundary
condition at the hot wall. Temperature of the super-heated wall is 10K higher than saturation temper-
ature.
The integrated simulation error is estimated as the film thickness error (|𝛿𝑠𝑖𝑚 − 𝛿𝑎𝑛 |) summed over
time steps 𝑖, weighted by time step (Δ𝑡):
∑
𝐸= |𝛿𝑠𝑖𝑚 − 𝛿𝑎𝑛 |Δ𝑡. (32)
𝑖
14 SAMKHANIANI AND ANSARI

FIGURE 8 Schematic of Stefan problem, boundary conditions, 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎, 𝑇𝑠𝑎𝑡 = 453.03𝐾, Δ𝑇sup = 10𝐾

TABLE 3 Stefan problem convergence

Numerical Simulation
Cell Number 128 256 512 Exact
Lee Final Thickness (mm) 1.88514 1.88228 1.87743 1.86722
Error (mm.s) 0.83902 0.65766 0.48976
Tanasawa Final Thickness (mm) 1.88093 1.88001 1.88084
Error (mm.s) 0.51143 0.53373 0.62118

FIGURE 9 Interface position in Stefan problem, compare exact solution with present numerical simulations (Lee
and Tanasawa), 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎 and 𝑟𝑒 = 105 , 𝛾 = 1.0

Simulation is conducted on three different grids with 128, 256, and 512 cells. The final interface
position is at 𝑡 = 50(𝑠) and the integrated errors are display in Table 3. It shows the numerical error in
the Tanasawa model is smaller than the Lee model on coarse grids. Furthermore, the numerical error
becomes smaller on finer grids in the Lee model, but remains constant in the Tanasawa model.
𝑇 −𝑇
The interface position (𝛿) and dimensionless temperature (𝜃 = 𝑇 −𝑇𝑤 ) are displayed in Figs. 9 and
𝑠𝑎𝑡 𝑤
10, respectively. It can be seen there is excellent agreement between CFD results and the exact solution.
SAMKHANIANI AND ANSARI 15

FIGURE 10 Dimensionless temperature in Stefan problem, compare exact solution with present numerical sim-
ulations (Lee and Tanasawa), 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎 and 𝑟𝑒 = 105 , 𝛾 = 1.0

FIGURE 11 Hydrodynamic transition in bubble release pattern observed by Reiman and Grigull63

3.2 2D-film boiling

Boiling of quiescent liquid near a hot solid surface is known as pool boiling. Some major regimes
in pool boiling are convective heat transfer, nucleation boiling, and film boiling. In film boiling, a
vapor layer separates saturated or sub-cooled liquid from the hot solid plate. Due to the liquid’s higher
density compared to vapor, Rayleigh–Taylor instability occurs which amplifies small perturbations at
the interface and leads to bubble growth. Figure 11 indicates that the bubble is released individually
and sequentially for the low heat flux situation, and columns are found for the high heat flux one.
This problem has no analytical solution. Therefore, numerical simulations are compared with previous
numerical work using the sharp interface method and the widely accepted correlations of Berenson62
for the Nusselt number
( )
𝜌 (𝜌 − 𝜌𝐺 )𝑔𝐻𝐿𝐺 1∕4 3∕4
𝑁𝑢𝐵𝑒 = 0.425 𝐺 𝐿 𝜆 , (33)
𝑘𝐺 𝜇𝐺 [𝑇𝑤 − 𝑇𝑠𝑎𝑡 ]

where 𝜆 is the characteristic length and given by

√
𝜎
𝜆= . (34)
(𝜌𝐿 − 𝜌𝐺 )𝑔
16 SAMKHANIANI AND ANSARI

FIGURE 12 Two-dimensional film boiling, dash line is computational domain, 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎

For numerical simulation of two-dimensional film boiling, two scenarios have been reported in the
literature:

1. Quasi-periodic

Bubble pinch off occurs individually in the node16,64 or happens sequentially between node and
anti-node23,64 and creates a quasi-periodic pattern.

2. Bubble column

The vapor bubble is connected to the vapor film with a thin neck and creates a mushroom
shape.16,17,22,65 Some researchers have noted that 2D numerical simulations cannot be expected to
give an accurate representation of the 3D physics of film boiling. Furthermore, the two-dimensional
simulation only considers the components of the interface curvature lying in the computational plane,
whereas the out-of-plane components are neglected. In other word, surface tension driving hydrody-
namic instabilities usually cannot be captured in a realistic way17 and Rayleigh–Plateau instability
which plays a role in bubble pinch off may never show up in the 2D model.17,22
The sharp interface method was conducted in all previous studies on the numerical simulation of
film boiling,16,17,22,23,64,66 hence, the present numerical study using the diffuse interface method is
compared with them. Here, the film boiling problems are set similar to Tomar23 and Guo64 works
which are for saturated water at near critical pressure𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎. They simulated film boiling
using the couple of level set and VOF method (CLSVOF23 and VOSET64 ) with a sharp mass transfer
model. Additionally, Guo64 applied a PLIC method for interface reconstruction to gain more accuracy.
A schematic of two-dimensional film boiling is illustrated in Fig. 12. The super-heated wall temperature
is set 10k and 30K above the saturation temperature. The wall boundary condition is

𝜕𝑃𝑟𝑔ℎ 𝜕𝛼𝐿
⃖⃖⃗ = (0, 0),
at 𝑦 = 0 ∶ 𝑈 𝑇 = 𝑇𝑠𝑎𝑡 + Δ𝑇sup , = = 0. (35)
𝜕𝑦 𝜕𝑦
SAMKHANIANI AND ANSARI 17

The outflow boundary condition is

𝜕𝑇 𝜕𝛼 ⃖⃖⃗
𝜕𝑈
at 𝑦 = 𝐻 ∶ 𝑃𝑟𝑔ℎ = 𝑃𝑠𝑎𝑡 , = 𝐿 = = 0. (36)
𝜕𝑦 𝜕𝑦 𝜕𝑦
At the beginning of simulation, a vapor film is inserted in the computational domain. Vapor film
thickness is given by66
( )
𝜋𝑁𝑥
𝑦 = 𝑦𝑐 + 𝜀𝐶𝑜𝑠 , (37)
𝑊
where 𝑦𝐶 , 𝜀 and 𝑁 are unperturbed film thickness, perturbation amplitude and perturbation wave,
respectively. The following values are chosen for the present simulation:
4.0 1 /
𝑦𝐶 = 𝜆2𝑑 , 𝜀 = 𝜆2𝑑 , 𝑁 = 1, 𝑊 = 𝜆2𝑑 2, (38)
128 128
where 𝜆2𝑑 is the most unstable Taylor in viscid wavelength for two-dimensional horizontal flat surfaces
and is defined as
√
3𝜎
𝜆2𝑑 = 2𝜋 . (39)
𝑔(𝜌𝐿 − 𝜌𝐺 )

Nusselt number is computed as the dimensionless heat flux through the wall and is given by
𝑊( )
𝜆 𝜕𝑇 |
∫ 𝑇𝑤 −𝑇𝑠𝑎𝑡
× |
𝜕𝑦 |𝑦=0
𝑑𝑥
0
𝑁𝑢2𝑑 = , (40)
𝑊
where 𝑊 is the computational domain width which is 𝜆2𝑑 ∕2. The computational domain height is 𝜆2𝑑
in the present simulations.
The grid independence test is conducted with different grid sizes using the Tanasawa mass transfer
model. Figure 13 depicts the interface as the first bubble is going to leave the vapor film. The bubble
detachment time from the thin vapor neck is dependent on grid size, because the interface breaks
whenever the thickness of neck falls below the specific resolved grids. Therefore, in the coarse grid
(120 × 240) the bubble leaves the vapor film sooner. There is no significant change in bubble shape
between the two last finer grids, so a grid with 240 × 480 cells is chosen for the film boiling simulations.
Figure 14 shows the quasi-periodic release of bubble from node and anti-node. The temperature of
the solid wall is set at 10K above saturation temperature. With the given initial interface, a bubble
is formed at the node (center). Then, the bubble leaves the film under the effect of buoyancy. After
detachment, the vapor neck returns to the wall due to surface tension and moves to the anti-node lead
bubble growth on the outer side. The present numerical result is in good agreement with the numerical
results of Refs. 64 and 23 in predicting the bubble release pattern, bubble shape, and first detachment
time of bubble. The space-averaged Nusselt number over time is shown in Fig. 15. The Nusselt number
depends strongly on the film thickness. Heat flux is greater where the vapor film is thin and less when
the film is the thick. Therefore, average heat flux and Nusselt number increase when the vapor rushes
to fill the bubble and the remaining film becomes thin. To the contrary, after detachment, the vapor
returns to the superheated wall and increased film thickness leads to a smaller Nusselt number. The
Nusselt number predicted by Berenson’s correlation is 3.69. The computed time and space averaged
Nusselt number is 2.57 and 2.64 for the Tanasawa and Lee models, respectively. The difference is
about 28% and 30%. Guo et al.64 and Tomar et al.23 obtained average Nusselt numbers of 4.67 and
18 SAMKHANIANI AND ANSARI

F I G U R E 1 3 Bubble shape at five grids before the first detachment at t = 0.32 s, the grids 210 × 420 and 240 ×
480 have almost the same shape, 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎, Δ𝑇sup = 10𝐾

FIGURE 14 Bubble shapes at different instants. 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎, Δ𝑇sup = 10𝐾

5.37 for a similar problem, corresponding to 27% and 45% difference, respectively. This shows the
present numerical errors are in the range of previous numerical simulations.
The second series of simulation were performed with a wall superheat of 30K above saturation
temperature. Figure 16 shows the interface evolution at different time instants. For this case, Lee’s mass
transfer model still predicts the quasi periodic release of the bubble, but Tanasawa’s mass transfer model
predicts a tall and stable vapor column which is similar with numerical studies of Refs. 23 and 64. The
difference may exist due to how these models calculate transferred mass flux rate during the boiling
process. In Lee’s model, phase change does not need any initial interface. The evaporation takes place
wherever in the liquid phase with the temperature above saturation temperature, thus several small
bubbles can also be seen in the CFD result in Fig. 16 using Lee’s model. In contrast, Tanaswa’s model
is similar to the sharp interface model, mass transfer only occurs across the existing interface. The space
SAMKHANIANI AND ANSARI 19

FIGURE 15 Simulation results (Lee and Tanasawa) for space-averaged Nusselt number compare with Berenson
correlation, 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎, Δ𝑇sup = 10𝐾

FIGURE 16 Bubble shapes at different instants. 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎, Δ𝑇sup = 30𝐾

averaged Nusselt number over time for this case is depicted in Fig. 17. It can be found that when vapor
film is stable, the space averaged Nusselt number is almost uniform. Lee’s model predicts a higher
Nusselt number, because after bubble detachment, the vapor rushes back to the wall and increases the
heat convection mode near the heated wall. The computed time and space averaged Nusselt number
is 2.04 and 2.71 for Tanasawa and Lee models, respectively. The Nusselt number from Berenson’s
correlation is 2.81, therefore, the difference is 44% and 26%, respectively.

3.3 Horizontal film condensation

There are two forms of condensation process: dropwise and filmwise condensation.67 Dropwise con-
densation usually occurs when the condensate does not wet the surface and therefore it takes the shape
of droplets instead of spreading into a thin film. The filmwise condensation generally prevails on
hydrophilic surfaces where a stable thin condensing film will be formed.
20 SAMKHANIANI AND ANSARI

FIGURE 17 Simulation results (Lee and Tanasawa) for space-averaged Nusselt number compare with Berenson
correlation, 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎, Δ𝑇sup = 30𝐾

FIGURE 18 Horizontal film condensation, 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎, Δ𝑇sup = 30𝐾

In this case, vapor condenses to form a liquid film on the top surfaces of an isothermal plate at 𝑇𝑤 .
By assuming a linear temperature profile from 𝑇𝑤 at the lower wall to 𝑇𝑠𝑎𝑡 at the interface, a control
volume analysis leads to an analytical solution for the film thickness 𝛿(𝑡):33
[ ( )−1 ]1∕2
𝑘 1 𝐻𝐿𝐺
𝛿𝑎𝑛 (𝑡) = 2𝑡 𝐿 + . (41)
𝜌𝐿 𝐶𝑝𝐿 2 𝐶𝑝𝐿 Δ𝑇

Relevant thermo physical properties are for saturated water at 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎. Simulation of fluids
with greater liquid to vapor property ratios and phase change enthalpies require finer meshes and
smaller time steps to adequately resolve the increased range of scales and heat and mass flux mag-
nitudes. Thus, the domain represents a section of an infinitely wide condensing film, a quasi 1D com-
putational domain with only one grid cell in the direction of the translational invariance is considered.
The top free stream boundary is set to saturation temperature 𝑇𝑠𝑎𝑡 = 453.03𝐾 and the bottom wall
is set 𝑇𝑤 = 423.03𝐾 which corresponds to a 30K subcooled temperature. All fluid entering through
SAMKHANIANI AND ANSARI 21

FIGURE 19 Film thickness in horizontal film condensation, compare exact solution with present numerical sim-
ulations (Lee and Tanasawa), 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎 and 𝑟𝑐 = 105 , 𝛾 = 1.0

FIGURE 20 Relative error in horizontal film condensation simulations (Lee and Tanasawa), 𝑃𝑠𝑎𝑡 = 1𝑀𝑝𝑎 and
𝑟𝑐 = 105 , 𝛾 = 1.0

the free stream boundary is vapor. In the beginning of simulation, a very thin film is inserted in the
computational domain.
The development of the liquid film thicknesses and relative errors are presented in Figs. 19 and 20,
respectively. The simulation relative error is defined as

𝛿𝑠𝑖𝑚 − 𝛿𝑎𝑛
𝑟𝐸 = . (42)
𝛿𝑎𝑛

Relative errors in predicted film thickness are higher during the beginning of simulations, because
the liquid films are initially under resolved on fixed meshes. Relative error is smoother and higher in
Lee’s model than Tanasawa’s model. Integrated simulation error was defined in Eq. (32) and presented
in Table 4. These results indicate an excellent agreement between both mass transfer models with the
analytical solution.
22 SAMKHANIANI AND ANSARI

TABLE 4 Convergence study on horizontal film condensation

Numerical Simulation
Cell Number 128 256 512 Exact
Lee film Thickness (mm) 0.46845 0.467253 0.467361 0.467829
Error (𝜇m.s) 11.3 8.82 9.09
Tanasawa film Thickness (mm) 0.469314 0.471797 0.47051
Error (𝜇m.s) 28.1 24.4 21.6

FIGURE 21 Laminar film condensation on a vertical plate, boundary conditions and grid

3.4 Laminar film condensation on a vertical plate

A sketch of the physical problem of laminar film condensation on a vertical plate is shown in Fig. 21.
Where z-axis is along the plate and the y-axis is perpendicular to the plate. The hydrophilic flat plate
maintained at a uniform subcooled temperature 𝑇𝑤 is suspended in a large volume of vapor at a sat-
uration temperature 𝑇𝑠𝑎𝑡 . Nusselt68 studied the steady laminar film condensation by (i) assuming the
vapor and liquid film are divided by a sharp interface, (ii) assuming a linear temperature profile across
the condensate, (iii) neglecting effects of inertia forces and interface shearing stress. The obtained
analytical expression for the film thickness (𝛿𝐹 ) is:
[ ]1∕4
4𝜇𝐿 𝑘𝐿 (𝑇𝑠𝑎𝑡 − 𝑇𝑤 )𝑧
𝛿𝐹 = . (43)
𝑔𝐻𝐿𝐺 𝜌𝐿 (𝜌𝐿 − 𝜌𝐺 )
SAMKHANIANI AND ANSARI 23

FIGURE 22 Comparison of film thickness given by present numerical simulations (Tanasawa and Lee) and Nus-
selt’s solution

And the velocity distribution in the condensate is given by:

𝑔(𝜌𝐿 − 𝜌𝐺 ) ( 1
)
𝑈 (𝑦) = 𝛿𝑦 − 𝑦2 . (44)
𝜇𝐿 2

It is assumed in Nusselt’s solution that the vapor velocity is zero.

In this section, laminar film condensation on the vertical plate is simulated and numerical results
are compared with Nusselt’s solution. The material property is for saturated water at 𝑃𝑠𝑎𝑡 = 21.9𝑀𝑝𝑎.
A 2D rectangular domain is employed with horizontal and vertical dimensions of 0.5 mm and 10 mm.
The mesh is defined with 100 × 300 cells which is grading in the horizontal and vertical direction with
24 SAMKHANIANI AND ANSARI

FIGURE 23 Comparison velocity distribution ax z = 5 mm given by present numerical simulations (Lee and
Tanasawa) and Nusselt’s solution

cells in a row near the wall or inlet is 0.1 the height of the cells row near free stream or outlet, respec-
tively. The wall is 20 K subcooled. The boundary conditions and mesh is illustrated at Fig. 21. For the
vapor side, total pressure is assigned for pressure and velocity is calculated from the pressure boundary
condition.
Simulation is conducted for 2 seconds to reach steady state condition. Assuming negligible inter-
shear, the wall shear stress can be evaluated as 𝜏𝑤 = 𝜌𝐿 𝑔𝛿 and dimensionless wall distance as
facial √
𝑔𝛿𝑦
𝑦+ = 𝜈𝐿
where 𝑦 is the cell center distance to the nearest wall and 𝜈𝐿 is the liquid kinematic vis-
cosity. In the present simulation 𝑦+ is smaller than 3. This low value of y+ indicates that near wall
fluid and heat transfer physics are well resolved. Figure 22 shows the liquid film thickness at steady
state determined from numerical result using Tanasawa’s and Lee’s models. It is seen that the liquid
film thickness obtained from the present numerical simulation agrees well with the classical Nusselt
analytical prediction.
The velocity profile from the numerical simulation is compared with Nusselt’s solution in Fig. 23.
The velocity increases from zero at the wall to a maximum value away from the wall and then decreases
monotonically to zero in the vapor at infinity. Thus, the velocity profile is continuous from the sub-
cooled liquid at the wall to dry vapor at infinity. The velocity profile given by the Nusselt model is
discontinuous. It can be seen that the velocity obtained by both mass transfer models are the same and
they are smaller than that of the classical Nusselt model given by Eq. (44). Because the classical Nus-
selt model did not take into consideration the effects of inertia forces, convection terms and interfaces
shear stress. In contrast, these three parameters are taken into consideration in the numerical simula-
tions. Therefore, the film liquid induces the vapor flow and at the same time, liquid film gives up some
of its momentum, and therefore its velocity decreases as compared to that given by Eq. (44). A similar
result for the velocity profile is pointed out in numerical results of Liu and Cheng69 who simulated
laminar film condensation with the Lattice Boltzman method.
𝑇 −𝑇
The dimensionless temperature defined as 𝜃 = 𝑇 −𝑇𝑤 is presented in Fig. 24. Nusselt’s solution
𝑠𝑎𝑡 𝑤
assumed a linear variation of temperature within the liquid film and constant saturation temperature
outside of the condensate liquid film. Numerical simulation shows for the present thermophysical prop-
erties and film thickness the linear variation of temperature in the liquid film is not valid.
SAMKHANIANI AND ANSARI 25

FIGURE 24 Temperature profiles at z = 5 mm in film condensation on a vertical flat plate, compare present
numerical simulations (Tanasawa and Lee) and Nusselt’s solution

FIGURE 25 Schematic of bubble condensation in subcooled boiling, initial and boundary conditions

3.5 Bubble condensation

Bubble condensation is one of the fundamental issues in a subcooled flow boiling to describe the heat
and mass transfer. In this section, the rising of a single vapor bubble in quiescent subcooled water
is simulated. The geometry of the considered problem and applied boundary condition are illustrated
26 SAMKHANIANI AND ANSARI

FIGURE 26 Bubble shape at 𝑡 = 1𝑚𝑠 for different grids, 𝑃𝑠𝑎𝑡 = 0.130[𝑀𝑝𝑎] and Δ𝑇𝑠𝑢𝑏 = 25[𝐾]

FIGURE 27 Bubble diameter history, compare between present numerical simulations and experimental data,70
𝐷0 = 1.008𝑚𝑚, 𝑃𝑠𝑎𝑡 = 0.130[𝑀𝑝𝑎] and Δ𝑇𝑠𝑢𝑏 = 25[𝐾]

in Fig. 25 which is similar to Samkhaniani and Ansari.71 The 2D space domain is set as 2𝐷0 × 4𝐷0
where 𝐷0 is the initial diameter of vapor bubble. The bubble is located in the position of (𝐷0 , 𝐷0 ) at
the beginning of the simulation.
A mesh resolution analysis is performed on four different grids using Tanasawa’s model to figure
out the appropriate grid size for simulation of bubble condensation. These grids are 50 × 100, 75 ×
150, 100 × 200, and 150 × 300 which, respectively, represent 25, 37, 50, and 75 cells across the initial
diameter. It can be seen obviously in Fig. 26 that bubble shape becomes stable with meshes of 100 ×
200 and 150 × 300. The variation between the two lateral grids is negligible. Therefore, the bubble is
resolved using 100 × 200 cells in this study.
In order to validate the present numerical simulation, the shape and the bubble diameter history
is compared with published experimental results70 under saturation pressure 0.130 Mpa for water
SAMKHANIANI AND ANSARI 27

FIGURE 28 Bubble shape, comparison between present numerical results (Lee and Tanasawa) and experimental
data,70 𝐷0 = 1.008𝑚𝑚, 𝑃𝑠𝑎𝑡 = 0.130[𝑀𝑝𝑎] and Δ𝑇𝑠𝑢𝑏 = 25[𝐾]

𝜌
corresponding to the severe density ratio of ( 𝜌𝐿 ≈ 1263). The comparison of bubble diameter history
𝐺
is shown in Fig. 27 and the bubble shape in Fig. 28. Bubble diameter in 2D simulation is the diameter
of an equivalent circle. It can be seen that numerical result is in reasonable good agreement with the
experimental results.

4. CONC LU SI ON

In the present study, a new solver called phaseChangeHeatFoam is implemented in the OpenFOAM
cfd package to simulate boiling and condensation processes. For now, two mass transfer models (Lee
and Tanasawa) are implemented in the solver. However, adding a new mass transfer model is straight
forward. The governing equations and programming information are discussed in detail. This solver is
based on the CF-VOF method, but spurious current is reduced by using a filter method for calculation
of the interface curvature. The present numerical simulations are compared with analytical solutions,
previous sharp interface method results and experimental data for different thermo-physical properties.
Our investigation shows the present diffuse interface method is capable of accurate simulation of the
𝜌
phase change process even for a severe ratio of thermophysical properties ( 𝜌𝐿 ≈ 1000). Additionally,
𝑣
the performance of Lee and Tanasawa mass transfer models are investigated here through simulations,
it was found both mass transfer models have similar performance in simple simulations (Stefan and
horizontal film condensation). However, in the simulation of complex phenomena, for instance film
boiling, the difference in the numerical results obtained with Lee and Tanasawa mass transfer models
can be considerable. Our study shows numerical results with Tanasawa mass transfer model are very
similar to the sharp interface model.
28 SAMKHANIANI AND ANSARI

APPENDIX A : THE R MO PH YS ICA L PRO P E RT I E S O F SAT U R AT E D WAT E R

Liquid
𝒌𝒈
𝑻 [𝑲] 𝑷 [𝑴𝒑𝒂] 𝝆[ 𝒎𝟑 ] 𝑯[ 𝒌𝑱
𝒌𝒈
] 𝒌𝑱
𝑪𝒑 [ 𝒌𝒈.𝑲 ] 𝝁[𝑷 𝒂.𝒔] 𝒘
𝒌[ 𝒎.𝑲 ] 𝝈[ 𝑵
𝒎
]
380.26 0.13 953.13 449.19 4.2244 2.620E-04 0.68106 0.05753
453.03 1 887.13 762.52 4.4045 1.5024E-4 0.67337 0.042217
646 21.9 402.4 1963.5 218 4.67E-05 0.545 0.00007
Vapor
𝒌𝒈
𝑻 [𝑲] 𝑷 [𝑴𝒑𝒂] 𝝆[ 𝒎𝟑 ] 𝑯[ 𝒌𝑱
𝒌𝒈
] 𝒌𝑱
𝑪𝒑 [ 𝒌𝒈.𝑲 ] 𝝁[𝑷 𝒂.𝒔] 𝒘
𝒌[ 𝒎.𝑲 ]
380.26 0.13 0.7545 2686.6 2.1107 1.251E-05 0.02590
453.03 1 5.1450 2777.1 2.7114 1.502E-05 0.03643
646 21.9 242.7 2240 352 3.23E-06 0.538

REFE RENCES
1. Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon
Not R Astron Soc. 1977;181:375–389.
2. Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys.
1992;100:25–37.
3. Sussman M, Smereka P, Osher S. A level set approach for computing solutions to incompressible two-phase flow.
J Comput Phys. 1995;114:146–159.
4. Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys.
1981;39:201–225.
5. Worner M. Numerical modeling of multiphase flow in microfluidics and micro process engineering: a review of
methods and applications. Microfluid NanoFluid. 2012;12:841–886.
6. Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric
incompressible two-phase flows. J Comput Phys. 2000;162:301–337.
7. Sussman M. A second order coupled level set and volume-of-fluid method for computing growth and collapse of
vapor bubbles. J Comput Phys. 2003;187:110–136.
8. Van der Pijl S, Segal A, Vuik C, Wesseling P. A mass conserving level set method for modelling of multi-phase
flows. Int J Numer Methods Fluids. 2005;47:339–361.
9. Lv X, Zou Q, Zhao Y, Reeve D. A novel coupled level set and volume of fluid method for sharp interface capturing
on 3D tetrahedral grids. J Comput Phys. 2010;229:2573–2604.
10. Sun D, Tao W. A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase
flows. Int J Heat Mass Transfer. 2010;53:645–655.
11. Youngs DL. Time-dependent multi-material flow with large fluid distortion. Num Methods Fluid Dyn. 1982;24:273–
285.
12. Renardy Y, Renardy M. PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method.
J Comput Phys. 2002;183:400–421.
13. Ubbink O, Issa R. A method for capturing sharp fluid interfaces on arbitrary meshes. J Comput Phys. 1999;153:26–
50.
14. Xiao F, Honma Y, Kono T. A simple algebraic interface capturing scheme using hyperbolic tangent function. Int
J Numer Methods Fluids. 2005;48:1023–1040.
15. Lyklema J. Fundamentals of Interfaces and Colloid Science, Vol. I: Fundamentals. Academic Press; 1991.
16. Welch SWJ, Wilson J. A volume of fluid based method for fluid flows with phase change. J Comput Phys.
2000;160:662–682.
SAMKHANIANI AND ANSARI 29

17. Hardt S, Wondra F. Evaporation model for interfacial flows based on a continuum-field representation of the source
terms. J Comput Phys. 2008;227:5871–5895.
18. Schlottke J, Weigand B. Direct numerical simulation of evaporating droplets. J Comput Phys. 2008;227 5215–5237.
19. Son G, Dhir V. Numerical simulation of film boiling near critical pressures with a level set method. J Heat Transfer.
1998;120:183–192.
20. Son G, Dhir V, Ramanujapu N. Dynamics and heat transfer associated with a single bubble during nucleate boiling
on a horizontal surface. J Heat Transfer. 1999;121:623–631.
21. Tanguy S, Ménard T, Berlemont A. A level set method for vaporizing two-phase flows. J Comput Phys.
2007;221:837–853.
22. Gibou FDR, Chen L, Nguyen D, Banerjee S. A level set based sharp interface method for the multiphase incom-
pressible Navier–Stokes equations with phase change. Comput Phys. 2007;222:536–555.
23. Tomar G, Biswas G, Sharma A, Agrawal A. Numerical simulation of bubble growth in film boiling using a coupled
level-set and volume-of-fluid method. Phys Fluids. 2005;17.
24. Kunkelmann C, Stephan P. Numerical simulation of the transient heat transfer during nucleate boiling of refrigerant
HFE-7100. Int J Refrig. 2010;33:1221–1228.
25. Weller HG, TaboraI G, Jasak H, Fureby C. A tensorial approach to computational continuum mechanics using
object-oriented techniques. Comput Phys. 1998;12:620–632.
26. Deshpande SS, Anumolu L, Trujillo MF. Evaluating the performance of the two-phase flow solver interFoam.
Comput Sci Discovery. 2012;5:014016.
27. Berberović E, Hinsberg NPv, Jakirlić S, Roisman IV, Tropea C. Drop impact onto a liquid layer of finite thickness:
dynamics of the cavity evolution. Am Phys Soc. 2009;79:036306 (15).
28. Klostermann J, Schaake K, Schwarze R. Numerical simulation of a single rising bubble by VOF with surface com-
pression. Int J Numer Methods Fluids. 2013;71:960–982.
29. Ubbink O. Numerical Prediction of Two Fluid Systems with Sharp Interfaces. London: University of London; 1997.
30. Kunkelmann C, Stephan P. CFD simulation of boiling flows using the volume-of-fluid method within OpenFOAM.
Num Heat Transfer 2009;56:631–646.
31. Zeng Q, Cai J, Yin H, Yang X, Watanabe T. Numerical simulation of single bubble condensation in subcooled flow
using OpenFOAM. Progress Nucl Energy. 2015;83:336–346.
32. Doro EO. Computational modeling of falling liquid film free surface evaporation. 2012.
33. Rattner AS, Garimella S. Simple mechanistically consistent formulation for volume-of-fluid based computations of
condensing flows. J Heat Transfer. 2014;136:071501.
34. Bahreini M, Ramiar A, Ranjbar AA. Numerical simulation of bubble behavior in subcooled flow boiling under
velocity and temperature gradient. Nucl Eng Des. 2015;293:238–248.
35. Schrage RW. A Theoretical Study of Interphase Mass Transfer. Columbia University Press; 1953.
36. Knudsen M. The Kinetic Theory of Gases: Some Modern Aspects. Metheun & Company; 1950.
37. Marek R, Straub J. Analysis of the evaporation coefficient and the condensation coefficient of water. Int J Heat Mass
Transfer. 2001;44:39–53.
38. Tanasawa I. Advances in condensation heat transfer. Adv Heat Transfer. 1991;21:55–139.
39. Lee WH. A pressure iteration scheme for two-phase flow modeling. Multiph Transp Fundam React Saf Appl.
1980;1:407–431.
40. Sun D-L, Xu J-L, Wang L. Development of a vapor–liquid phase change model for volume-of-fluid method in
FLUENT. Int. Commun Heat Mass Transfer. 2012;39:1101–1106.
41. Wu H, Peng X, Ye P, Gong YE. Simulation of refrigerant flow boiling in serpentine tubes. Int J Heat Mass Transfer.
2007;50:1186–1195.
42. De Schepper SC, Heynderickx GJ, Marin GB. Modeling the evaporation of a hydrocarbon feedstock in the convec-
tion section of a steam cracker. Comput Chem Eng. 2009;33:122–132.
30 SAMKHANIANI AND ANSARI

43. Alizadehdakhel A, Rahimi M, Alsairafi AA. CFD modeling of flow and heat transfer in a thermosyphon. Int Commun
Heat Mass Transfer. 2010;37:312–318.
44. Yang Z, Peng X, Ye P. Numerical and experimental investigation of two phase flow during boiling in a coiled tube.
Int J Heat Mass Transfer. 2008;51:1003–1016.
45. Goodson K, Rogacs A, David M, Fang C. Volume of fluid simulation of boiling two-phase flow in a vapor-venting
microchannel. Frontiers Heat Mass Transfer. 2010;1.
46. Lee H, Kharangate CR, Mascarenhas N, Park I, Mudawar I. Experimental and computational investigation of vertical
downflow condensation. Int J Heat Mass Transfer. 2015;85:865–879.
47. Albadawi A, Donoghue D, Robinson A, Murray D, Delauré Y. Influence of surface tension implementation in
volume of fluid and coupled volume of fluid with level set methods for bubble growth and detachment. Int J Multiph
Flow. 2013;53:11–28.
48. Samkhaniani N, Gharehbaghi A, Ahmadi Z. Numerical simulation of reaction injection molding with polyurethane
foam. J Cell Plast. 2013;49:405–421.
49. Marschall H, Hinterberger K, Schuler C, Habla F, Hinrichsen O. Numerical simulation of species transfer across
fluid interfaces in free-surface flows using OpenFOAM. Chem Eng Sci. 2012;78:111–127.
50. Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. J Comput Phys.
1992;100:335–354.
51. Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial flow. Ann Rev Fluid Mech.
1999;31:567–603.
52. Lafaurie B, Nardone C, Scardovelli R, Zaleski S, Zanetti G. Modelling merging and fragmentation in multiphase
flows with SURFER. J Comput Phys. 1994;113:134–147.
53. Hoang DA, van Steijn V, Portela LM, Kreutzer MT, Kleijn CR. Benchmark numerical simulations of segmented
two-phase flows in microchannels using the volume of fluid method. Comput Fluids. 2013;86:28–36.
54. Zalesak ST. Fully multidimensional flux-corrected transport algorithms for fluids. J Comput Phys. 1979;31:335–
362.
55. Issa RI. Solution of the implicitly discretised fluid flow equations by operator-splitting. J Comput Phys. 1986;62:40–
65.
56. Rhie C, Chow W. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J.
1983;21:1525–1532.
57. Van Leer B. Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in
a second-order scheme. J Comput Phys. 1974;14:361–370.
58. Weller H. A new approach to VOF-based interface capturing methods for incompressible and compressible flow.
OpenCFD Ltd, Report TR/HGW/04. 2008.
59. Jasak H. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows: University
of London; 1996.
60. Sun D-L, Xu J-L, Wang L. Development of a vapor–liquid phase change model for volume-of-fluid method in
FLUENT. International Communications in Heat and Mass Transfer 2012;39:1101–1106.
61. Tsujimoto K, Kambayashi Y, Shakouchi T, Ando T. Numerical simulation of gas-liquid two-phase flow with phase
change using Cahn-Hilliard equation. Turbulence Heat Mass Transfer 2009;6.
62. Berenson P. Film-boiling heat transfer from a horizontal surface. J Heat Transfer. 1961;83:351–356.
63. Reimann M, Grigull U. Warmeubergang bei freier konvektion und flimsieden im kritischen gebiet von wasser und
kohlendioxid. Warmeund Stof-fubertragung. 1975;8:229–239.
64. Guo DZ, Sun DL, Li ZY, Tao WQ. Phase change heat transfer simulation for boiling bubbles arising from a vapor
film by the VOSET method. Num Heat Transfer. 2011;59:857–881.
65. Juric D, Tryggvason Gt. Computations of boiling flows. Int J Multiphase Flow. 1998;24:387–410.
66. Esmaeeli A, Tryggvason G. Computations of film boiling. Part I: numerical method. Int J Heat Mass Transfer.
2004;47
SAMKHANIANI AND ANSARI 31

67. Ghiaasiaan SM. Two-Phase Flow, Boiling and Condensation in Conventional and Miniature System. Cambridge
University Press; 2008.
68. Nusselt W. The surface condensation of water vapour. VDI Z. 1916;60:541–546.
69. Liu X, Cheng P. Lattice Boltzmann simulation of steady laminar film condensation on a vertical hydrophilic sub-
cooled flat plate. Int J Heat Mass Transfer. 2013;62:507–514.
70. Kamei S, Hirata M. Condensing phenomena of a single vapor bubble into subcooled water. Exp Heat Transfer.
1990;3:173–182.
71. Samkhaniani N, Ansari M. Numerical simulation of bubble condensation using CF-VOF. Progress in Nuclear
Energy. 2016;89:120–131.

How to cite this article: Samkhaniani N, Ansari MR. The evaluation of the diffuse inter-
face method for phase change simulations using OpenFOAM. Heat Transfer–Asian Research.
2017;00:1–31. doi:10.1002/htj.21268