Document Type : Full Lenght Research Article
Authors
Faculty of Mechanical Engineering, Babol Noshirvani University of Technology, Babol, Iran
Abstract
Keywords
Main Subjects
History: Submitted 6 September 2015 Revised 16 November 2016 Accepted 28 December 2016
Keywords:
Phasechange model Volumeoffluid Boiling Condensation OpenFOAM 
A B S T R A C T In this study, the volumeoffluid method in the OpenFOAM opensource CFD package will be extended to consider the phasechange phenomena with a modified model due to the condensation and boiling processes. This model is suitable for the case, in which the unsaturated and saturated phases both are present and for which an initial interface is not needed to begin boiling and condensation processes. Both phases (liquid and vapor) are incompressible and immiscible. The interface between two phases is tracked with the colorfunction volumeoffluid (CFVOF) method. Surface tension is taken into consideration by the continuous surface force (CSF) model. Pressurevelocity coupling will be solved with the PISO algorithm in the collocated grid. The accuracy of this phasechange model is verified by two evaporation problems (a onedimensional Stefan problem and a twodimensional filmboiling problem) and two condensation problems (a onedimensional Stefan problem and filmwise condensation). The simulation results of this model show good agreement with the classical analytical or numerical results, proving its accuracy and feasibility. © 2016 Published by Semnan University Press. All rights reserved.

Development of a PhaseChange Model for the VolumeofFluid Method in OpenFOAM
Mohammad Bahreini*, Abbas Ramiar, Ali Akbar Ranjbar Faculty of Mechanical Engineering, Babol Noshirvani University of Technology, Babol, Iran 

Journal of Heat and Mass Transfer Research
Journal homepage: http://jhmtr.journals.semnan.ac.ir 

Corresponding Author: M. Bahreini, Faculty of Mechanical Engineering, Babol Noshirvani University of Technology, Babol, Iran Email: m.bahreini1990@gmail.com

Many methods have been proposed to simulate the vapor–liquid phasechange problems, such as the volumeoffluid (VOF) method [13] and the level set (LS) method [46]. The VOF method has an inherent massconservation property, more easily capturing the interface with the phase change’s heat transfer [7, 8]. The massconservation feature is particularly important when solving phasechange problems [8]. Therefore, the VOF method is a good choice. At present, the VOF method has been employed in the OpenFOAM 220 code to solve twophase flows. However, the default VOF method employed in OpenFOAM cannot simulate heat and mass transfer through the phase interface. To overcome this shortcoming, the phasechange model needs to be added to the source terms in the governing equations and implemented in the code, and a new solver needs to be proposed. There are many different kinds of phasechange models in the references. In the existing models, the vapor–liquid phasechange model proposed by Hardt [9] has been most widely used. Its expressions are shown below.
The evaporating mass flux at the liquid–vapor interface is calculated as
. 
(1) 
Herein, is the heat resistance of the liquid–vapor interface. Based on the considerations of Schrage [10], this variable can be expressed by
. 
(2)

is the evaporation coefficient. It represents the fraction of the fraction of the molecules that depart the interface during the evaporation, and its value depends on the fluid. Its range is from 0.04 to 1.0 [11]. Details on the procedure are given by Hardt and Wondra [9].
The interfacial temperature usually is set as the saturation temperature in the phasechange process. Based on Fourier’s law, the interfacial heatflux jump can be calculated by the following expression:
, 
(3) 
where is the interfacial unit normal vector and it points toward the vapor phase. Nichita and Thome [12] and Mao[13] Worked on the following vapor–liquid phasechange model was derived according to Eq. (3):
. 
(4) 
Due to unreasonable assumptions in the process of derivation, there is a large deviation between the computational results and the actual physical phenomena. For example, the bubble growth rate is not relevant to the vapor thermal conductivity in the growing process of a saturated bubble in a superheated liquid. However, Eq. (4) contains the information of , which does not match the actual physical process.
Based on the VOF method, many authors, such as Welch and Wilson [7] and Guo et al. [14], have developed program codes for solving phasechange problems. The key point of these methods is how to calculate the interfacial heat flux on both sides accurately, for example, and in Eq. (3). These methods can simulate the boiling and condensation problems that have an initial interface, but with the Lee phasechange model, boiling and condensation can be simulated without a primary interface, and when the temperature is higher than the saturation temperature, the phase change occurs.
In this study, a multiphase solver capable of predicting the phase change is developed by extending the VOF method currently implemented in the OpenFOAM220[15] software package (interFoam) and by adding an energy equation while considering proper source terms. The implemented masstransfer model has been modified to improve the simulation accuracy. Also, using the color function in conjunction with the VOF method eliminates the need for an initial interface to start the boiling and condensation processes. The phasechange process is governed by the saturation temperature. Furthermore, developing new models using the OpenFOAM platform has the advantage of employing unstructured grids and parallel processing. Thus, the proposed method is not limited to any geometry, physics, or flow regime. The new solver is validated against three different cases, including the onedimensional phase change, twodimensional film boiling, and condensation film flow.
2. Computational method mathematical models
2.1. Governing equations
The VOF method uses a discontinuous scalar function to resolve the interface in fixed grids. This scalar function is the ratio of one fluid volume to the volume of the cell, and it is defined as The physical properties of the twophase flow, such as viscosity ( ), density ( ), thermal conductivity ( ), and specific heat capacity ( ),
(5)

are defined as
, 
(6) 
where .
The following equations must be solved in the boiling simulation. These equations are the conservations for mass, momentum, energy, and volume fraction:
, 
(7) 
, 
(8) 
, 
(9) 
, 
(10) 
where defines the surface tension between two phases, without the density average, determined by the continuum surface force (CSF) model, proposed by [16]. is the surface tension, and is the interface curvature given by
(11) 
. 
The extra divergence term in Eq. (6) contributes only at the region of the interface ( ). It limits the smearing of the interface because of the compensation of the diffusive fluxes [17]. is the compressive velocity, which is calculated in the normal direction of the interface to avoid any dispersion [18]. Moreover, a compressive factor is used to increase compression as
(12) 
. 
In the present study, the compression factor is considered. The coefficient controls the weight of the compression flux and usually should be in the range of unity ( ) [17, 19] .
The Lee masstransfer model is implemented to account for evaporation and condensation. The mass transfer rate in volume (kg/m^{3}s) is given by [20, 21]
, 
(13) 
, 
(14) 
where and are the masstransfer intensity factors for evaporation and condensation, respectively. Different values can be assigned to and .[21, 22] An exclusive discussion about the masstransfer intensity factor can be found in the Appendix.
2.2. Numerical details
All computations are performed using OpenFOAM [23]. In each time step, the following tasks are performed by the solver:
1 advection of the VOF field and update of the thermophysical variables;
2 an initial solution of the velocity field;
3 calculation of the energy equation and an update of the masstransfer rate; and
4 a solution of the pressurevelocity coupling (pressure implicit with splitting of operators (PISO) algorithm).
The governing equations are discretized based on a finite volume (FVM) formulation. The discretization is performed on a structured grid. All the variables are stored on the cell centers, where a collocated grid arrangement is used. In order to avoid a checkerboarding effect in the momentum equation, the Rhie–Chow momentum interpolation [24] is applied.
In Table 1, the corresponding terminology of OpenFOAM also is given.
Table. 1 Discretization and interpolation schemes of the numerical model.
Term 
Discretization scheme 
Method 
Euler 
the firstorder bounded implicit scheme 

vanLeerV 
See text 

vanLeer 
See text 

Interface compression 
See text 

vanLeer 
See text 

Linear 
central difference schemes 

Corrected 
surfacenormal gradient with correction on nonorthogonal meshes [25] 

Linear corrected 
face values ( ) approximated by central difference scheme, and the resulting surfacenormal gradient is calculated using central difference schemes with nonorthogonal correction 

Linear 
default interpolation schemes for getting face values from cell values 
The spatial derivatives are discretized using secondorder schemes. The compressive term in Eq. (10) is discretized using the interface compression scheme [26] with the limiter function defined as
(15) 
where and are the flux at the evaluated and neighbor cells, respectively.
The discretization of other convective terms are performed using the vanLeer scheme [27] with the following flux limiter function:
(16) 
, 
where the smoothness parameter is defined by the ratio of consecutive gradients [28].
3. Result and discussion
3.1. Stephan Problem
3.1.1. Evaporation test
The onedimensional Stephan problem has become a benchmark for the phasechange model [7, 9, 29, 30]. In this problem, a vapor film separates the liquid from the superheated wall; the liquid and vapor are initially in quiescent equilibrium. Evaporation pushes the liquid away from the heated wall. The interface is at the saturation temperature. Fig. 1 shows the boundary conditions for this problem, and the physical properties of the Stephan problem are shown in Table 2. A vapor film is inserted in the computational domain in the initial time, which consists of four first cells near the hot wall.
The analytical solution of this problem is as follows [7]:
, 
(17) 
, 
(18) 
where x (t) is the exact interface position with time, T (x,t) is the temperature of the vapor side, is the vapor thermal conductivity, and is calculated from the following equation [7]:
, 
(19) 
or can be estimated with the following formula:
(20) 
In this simulation, a onedimensional computational domain is set up with only one grid cell in the direction of the translational invariance. The integrated simulation error is estimated as the vapor filmthickness error summed over time steps , weighted by the time step :
. 
(21) 
Fig. 1 Schematic shape of the Stefan problem and boundary conditions for evaporation.
Table. 2 Physical properties used in the Stefan problem.
Density (gas, liquid) 
1 
Viscosity (gas, liquid) 

Specific Heat (gas, liquid) 
1000 
Liquid Heat Conductivity 
1 
Gas Heat Conductivity 
0.1& 0.01 
Latent Heat 

Surface Tension 
0.1 
r (masstransfer intensity) 
1000 
(a) 
(b) 
Fig. 2 The interface positions at different times in the Stephan problem for the evaporation condition (a: b: ).
Fig. 2 shows the volume fraction distributions for two different thermal conductivities of the vapor phase. There is excellent agreement between the present numerical result and the exact solution based on the illustrated result in Fig. 3 and Table 3. Fig. 3 shows the interface position as a function of time for two different thermal conductivities of the vapor phase. The numerical results are displayed as lines and the analytical results as symbols. The curves indicate how the evaporation process is sped up with the increasing thermal conductivity of the vapor phase, whereupon excellent agreement between the analytical and the numerical results is obtained. Table 3 shows different mesh sizes—M1 (64), M2 (128), M3 (512), and M4 (1000)—and compares convergence with an exact solution. This study used 1000 grids for the Stephan problem.
Fig. 3 Comparison of the interface position between the analytical solution and computations.
Fig. 4 Temperature distribution along the domain for the evaporation condition at time = 100[s].
Fig. 4 shows the temperature distribution at t = 100[s] along the computational domain. On the vapor side, it predicts an analytical solution (Eq. (18)), and on the liquid side, it is correctly in the saturation temperature [7].
3.1.2. Condensation Test
We can use the Stephan problem for the condensation condition. In this situation, a liquid film separates the vapor from the cold wall; liquid and vapor are initially in quiescent equilibrium. Condensation pushes the vapor away from the cold wall. The interface is at the saturation temperature. Fig. 5 shows the boundary conditions for this problem. At the initial time, a liquid film is inserted in the computational domain, which consists of four first cells near the cold wall.
Fig. 6 shows the volume fraction distribution in the condensation condition for .
Table 3. Stefan problem convergence

Present Simulation (Grid Number) 



M1(64) 
M2(128) 
M3(512) 
M4(1000) 
Exact 
0.01443 
0.01426 
0.01414 
0.01411 
1.410 

0.0292 
0.0141 
0.0051 
0.0012 
 
Fig. 5 Stephan problem for condensation conditions.
Also, Fig.7 shows the temperature distribution at t = 100[s] along the computational domain. On the liquid side, it predicts the analytical solution (Eq. (18)), and on the vapor side, it is correctly in the superheat temperature [7].
Fig. 6 Interface positions at different times in the Stephan problem for condensation condition
3.2. Twodimensional film boiling problem
Boiling of a quiescent liquid near a hot solid surface is called pool boiling. Some major regimes in pool boiling are convective heat transfer, nucleation boiling,
Fig. 7 Temperature distribution along the domain for the condensation condition at time = 100[s].
and film boiling. The target of this section is the last one. In film boiling, the vapor layer separates a saturated or subcooled liquid from the hot solid plate. Since liquid has a higher density than vapor, Rayleigh–Taylor instability [31] occurs, which amplifies the small perturbation at the interface and leads to bubble growth:
. 
(22) 
Rayleigh–Taylor’s correlation assumes a uniform film thickness and a bubble diameter and height that are proportional to the bubble spacing. Conservation of mass and momentum in the vapor film, as well as the assumption that heat is transferred across the film due to conduction, enabled him to arrive at an expression for the Nusselt number [7]:
, 
(23) 
where the characteristic length is given by [7].
. 
(24) 
This correlation and the physical situation modeled by it are inherently threedimensional. We will present twodimensional simulations with the understanding that the numerical results cannot completely capture the physics of threedimensional film boiling, but these simulations do present useful test problems with which to develop methods directed toward simulating flows with mass transfer. We note also that there are correlations in the literature thought to be more accurate than the Berenson correlation [7]. We do not expect that twodimensional simulations can predict Nusselt numbers with the fidelity necessary to require a comparison to the more accurate correlations.
We consider a twophase fluid with the properties given in Table 4 that also are used in [7], and at the beginning, the interface has an initial position that is given by the following equation [11]:
. 
(25) 
Table. 4 Simulation details used in the filmboiling problem.
Density Liquid 
200 
Density Gas 
5 
Viscosity Liquid 

Viscosity Gas 

Specific Heat Liquid 
400 
Specific Heat Gas 
200 
Liquid Heat Conductivity 
40 
Gas Heat Conductivity 
1 
Latent Heat 

Surface Tension 
0.1 

0.0786 m 
Gravity 
9.81 
r (masstransfer intensity) 
10000 
Computational Grid (coarse) 

Computational Grid (refined) 
Fig. 8 shows the schematic diagram and boundary condition of the twodimensional boiling problem.
Fig. 8 Computational domain and boundary condition for filmboiling simulations.
Fig. 9 Quasiperiodic pattern in a filmboiling problem in coarse mesh.
Fig. 9 shows a configuration similar to that seen by Welch and Wilson [7], who observed a quasiperiodic pattern, but we observed a mushroom shape (Fig. 10). It should be mentioned that a group of authors observed a quasiperiodic pattern for bubble growth [7], but Hardt and Wondra [9] correctly noted it as a numerical artifact. This pattern can be observed only in coarse mesh, and in a mesh that is refined enough, a mushroom structure will be observed.
In our simulations, the local Nusselt number is calculated as the dimensionless heat flux through the wall:
(26) 


Fig. 10 Mushroom structure in a filmboiling problem in refined mesh.
Fig. 11 Nusselt number for a twodimensional film boiling problem in coarse mesh.
The Nusselt number in Fig. 11 is an averaged value over the width of the solid wall. The maximum and minimum Nusselt numbers correspond to minimum and maximum average film thickness, respectively. The corresponding timeaveraged Nusselt number is 2.7, which is only a 2% deviation from the result of the Berenson correlation (Eq. 24) [7, 9, 27].
3.3. Filmwise condensation on vertical channel
In this section, a computational model has been developed to predict the condensation heat transfer on a vertical channel. Filmwise condensation occurs when a cold wall surface is in contact with a vapor near saturation conditions [32, 33]. The condensation process begins with vapor condensing directly on the wall surface. However, in contrast with dropwise condensation, after the wall is initially wetted, it remains covered by a thin film of condensate. After that point, condensation occurs only at the liquid–vapor interface. Therefore, the condensation rate is directly a function of the rate at which heat is transported across the liquid film from the liquid–vapor interface to the wall. Fig. 12 shows three distinct regimes of filmwise condensation on a vertical wall. These regimes proceed in order from the top of the wall (x = 0): laminar, wavy, and turbulent. The Reynolds number is defined as , where Γ is the mass flow rate of condensate per unit width. At the top of the wall, where the film is thinnest, the laminar regime exists, and the liquid film is the laminar as the film’s Reynolds number is below 30 [34, 35]. As the condensation process proceeds, more and more condensation appears on the surface and the liquid condensate is pulled downward by gravity. As the condensate moves downward, the film becomes thicker. The first sign of transition to a non laminar regime appears as a series of regular ripples or waves of condensate. This regime is called the wavy regime ( ) and is considered neither laminar nor turbulent. It is characterized by a consistent, regular series of waves in time. Finally, if the wall is long enough, the film thickness becomes so great that irregular ripples in both time and space will appear, which are identified as a turbulent flow regime [35].
The saturated water vapor enters into the channel and the condensation takes place because the wall temperature is below the saturation temperature of the working fluid. Channel geometries are employed as follows: channel length L = 40 cm and width 2W = 4 cm. The flow both in vapor and liquid phases is laminar and incompressible. The geometric configuration and boundary condition of the condensation in a vertical channel are shown in Fig. 13, and the details of this problem are shown in Table 5.
Table. 5 Simulation details used in the filmwise condensation problem.
Density Liquid 
971.8 
Density Gas 
0.5981 
Viscosity Liquid 

Viscosity Gas 

Specific Heat Liquid 
4.1963 
Specific Heat Gas 
200 
Liquid Heat Conductivity 
0.67 
Gas Heat Conductivity 
0.0251 
Latent Heat 

Surface Tension 
0.0582 
Gravity 
9.81 
r (masstransfer intensity) 

Computational Grid 

Fig. 12 Flow regimes of film condensate on a vertical wall [35].
Fig. 13 Geometric configuration and boundary condition of the vertical plate.
The calculations were carried out primarily with three different grid densities: 0.2 million, 0.3 million, and 0.4 million elements to perform a grid independency test. The total condensation rate for the two fine grids is almost the same and 3.5% higher than that of the coarse grid. For the finest grids, the solution is too time consuming to be used practically. In this work, a mesh around 0.3 million quadrilateral grid elements is used for twodimensional simulations. In the liquid film flow regions and vapor–liquid interface regions, quadrilateralshaped grids are either parallel or normal to the flow and masstransfer direction. This means that the effect on the numerical diffusion is minimal and the simulation discrepancy is minimized. Other shaped grids, such as triangular grids, will lead to a difficulty in convergence. The accuracy of the simulations is determined by the grid shape and size. In the liquid film flow region, the momentum gradients are larger than those in other regions. Thus, a fine mesh is needed to reduce the discretization errors. A finer mesh can lead to a more accurate simulation but additional computing cost [36]. To balance simulation accuracy and computing cost, a mesh that is 0.5 mm along the streamwise direction and 0.05 mm along the mass transfer direction is adopted in this work, while the mesh is much coarser in the vapor flow region. In general, the convergence criterion of the constant is applied for all equations, where R_{n} refers to the maximum residual value summed over all the computation cells after the nth iteration. For all reported computational results, the convergence criterion is less than for velocities and less than for enthalpy.
Fig. 14 shows the volume fraction factor field of the liquid in the considered domain, and Fig. 15 shows the distribution along the Xaxis at the middle of the channel for the liquid void fraction and temperature. The sharp variations are found in the interfacial cells near the wall. The void fraction is equal to unity near the wall due to the liquid film flowing along the wall. A sharp decrease in the void fraction along the X axis indicates the transition from a liquid to a vapor through the thin interface region, as shown in Fig. 15 a. The temperature is shared among the phases, and there is a sharp increase in the temperature along the axis because the vapor’s temperature is fixed at the saturation temperature and the variation only exists in the thin liquid film region (see Fig. 15b).
Fig. 16 shows a comparison between the results of the present calculations and the results using the source terms adopted in Zhenyu et al. [34] and Aghanajafi et al. [37]. The source terms based on the kinetic gas theory predict a similar distribution of the heattransfer coefficient using source terms based on the energy balance at the interface [37] for the whole wall length. At this position, the film thickness is the main influential factor for the thermal resistance between the saturated vapor and the cold wall. A high heattransfer coefficient is achieved because the film thickness is the thinnest near the inlet of the channel. As the condensate moves downward, the film becomes thicker; therefore, the heattransfer coefficient decreases, as expected. As ripples or waves of the condensate appear, the heat transfer will be enhanced due to the increase in the interface area and the fluctuation of the condensate film.
Fig. 14 Volume fraction factor of liquid and temperature field in the channel.
(a) 

(b) 
Fig. 15 a: Liquidphase void fraction factor; b: temperature
Fig. 16 Comparison of heattransfer coefficient between present numerical simulations
and experimental data of [34, 37].
4. Conclusion
The modified evaporation–condensation model was implemented into the VOF scheme with a piecewise linear interface construction as available within the commercial CFD code OpenFOAM. The model is suitable for the case in which one phase is unsaturated and the other is saturated. Calculations were performed for three analytically or semianalytically solvable test cases, a onedimensional Stefan problem, and filmwise condensation. In addition, a twodimensional film boiling problem was considered that has been analyzed by other
authors based on the VOF and the levelset technique.
Good agreement between the numerical and the analytical results was found for the Stefan problems and filmwise condensation on the vertical channel. In the filmboiling problem, a predominant evaporation mode with mushroomshaped bubbles connected to the vapor film through a thin neck was identified. The same evaporation mode has been found by other authors [9] but stands in conflict with quasiperiodic bubble release patterns determined in earlier work [7]. The application of the model to the filmboiling problem demonstrates its usefulness in more complex situations than in the other validation cases.
Appendix
One of the wellknown equations for calculating net mass flux in boiling and condensation through the interface is the HertzKnudsenScharge equation [10], which is defined as
(27) 

Impinging molecules from the vapor phase can be reflected back into the vapor or can be absorbed by the liquid. Therefore, is considered to be a function of the vapor molecules numbers that collide interface and condense or a function of liquid molecules numbers that collide interface and evaporate [38].
Pressure and temperature are related in a saturation condition. While the temperature is close to a saturation temperature, the Clausius–Clapeyron equation is used to define it as follows:
(28) 
Substitution of Eq. (28) in Eq. (27) leads to the following equation for net mass flux:
(29) 
The mass transfer should be expressed in kg/m^{3}s in conservative equations; therefore,
(30) 
where is the interfacial surface area. Hence, for evaporation would be
(31) 

(32) 
It can be concluded that and depend on material properties, and they are functions of the interfacial area. Therefore, the fixed mass transfer intensity factor might be considered as a source of error in the present numerical model [21].
References
[1] D.L. Youngs, Timedependent multimaterial flow with large fluid distortion, Numerical methods for fluid dynamics, 24 (1982) 273285.
[2] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics, 39 (1981) 201225.
[3] M. Huang, L. Wu, B. Chen, A piecewise linear interfacecapturing volumeoffluid method based on unstructured grids, Numerical Heat Transfer, Part B: Fundamentals, 61 (2012) 412437.
[4] S. Osher, R.P. Fedkiw, Level set methods: an overview and some recent results, Journal of Computational Physics, 169 (2001) 463502.
[5] P. Wang, H. Sun, P.Y. Wong, H. Fukuda, T. Ando, Modeling of dropletbased processing for the production of highperformance particulate materials using the level set method, Numerical Heat Transfer, Part A: Applications, 61 (2012) 401416.
[6] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible twophase flow, Journal of Computational physics, 114 (1994) 146159.
[7] S.W. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, Journal of computational physics, 160 (2000) 662682.
[8] M.W. Akhtar, S.J. Kleis, Boiling flow simulations on adaptive octree grids, International Journal of Multiphase Flow, 53 (2013) 8899.
[9] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuumfield representation of the source terms, Journal of Computational Physics, 227 (2008) 58715895.
[10] R.W. Schrage, A theoretical study of interphase mass transfer, Columbia University Press, 1953.
[11] D. Juric, G. Tryggvason, Computations of boiling flows, International Journal of Multiphase Flow, 24 (1998) 387410.
[12] B.A. Nichita, J.R. Thome, A level set method and a heat transfer model implemented into FLUENT for modeling of microscale two phase flows, in: AVT178 Specialists" Meeting on System Level Thermal Management for Enhanced Platform Efficiency, 2010.
[13] W. Mao, Numerical Simulation of Vapor–liquid Phase Change Heat Transfer and Micromixing in Microfluidic Systems, Master's thesis, GuangZhou Institute of Energy Conversion Chinese Academy of Sciences, China, (2009).
[14] D. Guo, D. Sun, Z. Li, W. Tao, Phase change heat transfer simulation for boiling bubbles arising from a vapor film by the VOSET method, Numerical Heat Transfer, Part A: Applications, 59 (2011) 857881.
[15] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using objectoriented techniques, Computers in physics, 12 (1998) 620631.
[16] J.U. Brackbill, D.B. Kothe, C. Zemach, A Continuum Method for Modeling Surface Tension, J. Comput. Phys., 100 (1992) 335354.
[17] J. Klostermann, K. Schaake, R. Schwarze, Numerical simulation of a single rising bubble by VOF with surface compression, International Journal for Numerical Methods in Fluids, 71 (2013) 960982.
[18] A. Albadawi, D. Donoghue, A. Robinson, D. Murray, Y. Delauré, Influence of surface tension implementation in volume of fluid and coupled volume of fluid with level set methods for bubble growth and detachment, International Journal of Multiphase Flow, 53 (2013) 1128.
[19] N. Samkhaniani, A. Gharehbaghi, Z. Ahmadi, Numerical simulation of reaction injection molding with polyurethane foam, Journal of Cellular Plastics, 49 (2013) 405421.
[20] W.H. Lee, Pressure iteration scheme for twophase flow modeling, IN" MULTIPHASE TRANSPORT: FUNDAMENTALS, REACTOR SAFETY, APPLICATIONS". (1980) 407432.
[21] M. Bahreini, A. Ramiar, A.A. Ranjbar, Numerical simulation of bubble behavior in subcooled flow boiling under velocity and temperature gradient, Nuclear Engineering and Design, 293 (2015) 238248.
[22] M. Bahreini, A. Ramiar, A. Ranjbar, Numerical simulation of subcooled flow boiling under conjugate heat transfer and microgravity condition in a vertical mini channel, Applied Thermal Engineering, 113 (2017) 170185.
[23] H.G. Weller, G. TaboraI, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using objectoriented techniques, Comput. Phys., 12 (1998) 620632.
[24] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA journal, 21 (1983) 15251532.
[25] H. Jasak, Error analysis and estimation for the finite volume method with applications to fluid flows, in: Imperial College, Vol. PhD Thesis, University of London, 1996.
[26] H. Weller, A new approach to VOFbased interface capturing methods for incompressible and compressible flow, OpenCFD Ltd., Report TR/HGW/04, (2008).
[27] B. Van Leer, Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a secondorder scheme, Journal of Computational Physics, 14 (1974) 361370.
[28] M. Darwish, F. Moukalled, TVD schemes for unstructured grids, International Journal of Heat and Mass Transfer, 46 (2003) 599611.
[29] B. Shu, F. Dammel, P. Stephan, Phase change model for twophase fluid flow based on the volume of fluid method, ICHMT DIGITAL LIBRARY ONLINE, 13 (2008).
[30] G. Son, V. Dhir, Numerical simulation of film boiling near critical pressures with a level set method, Journal of Heat Transfer, 120 (1998) 183192.
[31] P. Berenson, Filmboiling heat transfer from a horizontal surface, Journal of Heat Transfer, 83 (1961) 351356.
[32] A. Gavrish, V.G. Rifert, Investigation of transition mode phenomenon from dropwise to filmwise steam condensation, Heat Transfer Research, 30 (1999).
[33] V. Semenov, N. Nikitin, Condensation heat transfer on noncircular pipes in stationary vapor, Heat Transfer Research, 39 (2008).
[34] Z. Liu, B. Sunden, J. Yuan, VOF modeling and analysis of filmwise condensation between vertical parallel plates, Heat Transfer Research, 43 (2012).
[35] S.M. Ghiaasiaan, Convective heat and mass transfer, Cambridge University Press, 2011.
[36] Z. Liu, B. Sunden, J. Yuan, Numerical Modelling of Condensation of Multiple Bubbles in Subcooled Flow Boiling With VOF Method, in: ASME 2013 Heat Transfer Summer Conference collocated with the ASME 2013 7th International Conference on Energy Sustainability and the ASME 2013 11th International Conference on Fuel Cell Science, Engineering and Technology, American Society of Mechanical Engineers, 2013, pp. V001T003A014V001T003A014.
[37] C. Aghanajafi, K. Hesampour, Heat transfer analysis of a condensate flow by VOF method, Journal of Fusion Energy, 25 (2006) 219223.
[38] R. Marek, J. Straub, Analysis of the evaporation coefficient and the condensation coefficient of water, International Journal of Heat and Mass Transfer, 44 (2001) 3953.
