Development of a phase change model for volume-of-fluid method in OpenFOAM

Many methods have been proposed to simulate the vapor–liquid phase-change problems, such as the volume-of-fluid (VOF) method [1-3] and the level set (LS) method [4-6]. The VOF method has an inherent mass-conservation property, more easily capturing the interface with the phase change’s heat transfer [7, 8]. The mass-conservation feature is particularly important when solving phase-change 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 two-phase flows. However, the default VOF method employed in OpenFOAM cannot simulate heat and mass transfer through the phase interface. To overcome this shortcoming, the phase-change 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 phase-change models in the references. In the existing models, the vapor–liquid phase-change 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

Herein,  is the heat resistance of the liquid–vapor interface. Based on the considerations of Schrage [10], this variable can be expressed by

 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 phase-change process. Based on Fourier’s law, the interfacial heat-flux jump can be calculated by the following expression:

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 phase-change model was derived according to Eq. (3):

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 phase-change 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 phase-change 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 mass-transfer 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 phase-change 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 one-dimensional phase change, two-dimensional 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 two-phase flow, such as viscosity ( ), density ( ), thermal conductivity ( ), and specific heat capacity ( ),

are defined as

where .

The following equations must be solved in the boiling simulation. These equations are the conservations for mass, momentum, energy, and volume fraction:

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

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

 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 mass-transfer model is implemented to account for evaporation and condensation. The mass transfer rate in volume (kg/m³s) is given by [20, 21]

where and are the mass-transfer intensity factors for evaporation and condensation, respectively. Different values can be assigned to and .[21, 22] An exclusive discussion about the mass-transfer 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 mass-transfer rate; and

4- a solution of the pressure-velocity 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 checker-boarding 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.

The spatial derivatives are discretized using second-order schemes. The compressive term in Eq. (10) is discretized using the interface compression scheme [26] with the limiter function defined as

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:

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 one-dimensional Stephan problem has become a benchmark for the phase-change 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]:

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]:

 or can be estimated with the following formula:

In this simulation, a one-dimensional 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 film-thickness error  summed over time steps , weighted by the time step :

Fig. 1   Schematic shape of the Stefan problem and boundary conditions for evaporation.

Table. 2  Physical properties used in the Stefan problem.

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

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. Two-dimensional 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 sub-cooled 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:

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]:

where the characteristic length is given by [7].

This correlation and the physical situation modeled by it are inherently three-dimensional. We will present two-dimensional simulations with the understanding that the numerical results cannot completely capture the physics of three-dimensional 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 two-dimensional simulations can predict Nusselt numbers with the fidelity necessary to require a comparison to the more accurate correlations.

We consider a two-phase 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]:

Table. 4  Simulation details used in the film-boiling problem.

Fig. 8 shows the schematic diagram and boundary condition of the two-dimensional boiling problem.

Fig. 8  Computational domain and boundary condition for film-boiling simulations.

Fig. 9  Quasi-periodic pattern in a film-boiling problem in coarse mesh.

Fig. 9 shows a configuration similar to that seen by Welch and Wilson [7], who observed a quasi-periodic pattern, but we observed a mushroom shape (Fig. 10). It should be mentioned that a group of authors observed a quasi-periodic 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:

Fig. 10  Mushroom structure in a film-boiling problem in refined mesh.

Fig. 11  Nusselt number for a two-dimensional 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 time-averaged 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.

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 two-dimensional simulations. In the liquid film flow regions and vapor–liquid interface regions, quadrilateral-shaped grids are either parallel or normal to the flow and mass-transfer 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 heat-transfer 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 heat-transfer 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 heat-transfer 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.

Fig. 15  a: Liquid-phase void fraction factor; b: temperature

Fig. 16  Comparison of heat-transfer 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 semi-analytically solvable test cases, a one-dimensional Stefan problem, and filmwise condensation. In addition, a two-dimensional film boiling problem was considered that has been analyzed by other
authors based on the VOF and the level-set 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 film-boiling problem, a predominant evaporation mode with mushroom-shaped 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 quasi-periodic bubble release patterns determined in earlier work [7]. The application of the model to the film-boiling problem demonstrates its usefulness in more complex situations than in the other validation cases.

Appendix

One of the well-known equations for calculating net mass flux in boiling and condensation through the interface is the Hertz-Knudsen-Scharge equation [10], which is defined as

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:

Substitution of Eq. (28) in Eq. (27) leads to the following equation for net mass flux:

The mass transfer should be expressed in kg/m³s in conservative equations; therefore,

where is the interfacial surface area. Hence,  for evaporation would be

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, Time-dependent multi-material flow with large fluid distortion, Numerical methods for fluid dynamics, 24 (1982) 273-285.

[2] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics, 39 (1981) 201-225.

[3] M. Huang, L. Wu, B. Chen, A piecewise linear interface-capturing volume-of-fluid method based on unstructured grids, Numerical Heat Transfer, Part B: Fundamentals, 61 (2012) 412-437.

[4] S. Osher, R.P. Fedkiw, Level set methods: an overview and some recent results, Journal of Computational Physics, 169 (2001) 463-502.

[5] P. Wang, H. Sun, P.Y. Wong, H. Fukuda, T. Ando, Modeling of droplet-based processing for the production of high-performance particulate materials using the level set method, Numerical Heat Transfer, Part A: Applications, 61 (2012) 401-416.

[6] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, Journal of Computational physics, 114 (1994) 146-159.

[7] S.W. Welch, J. Wilson, A volume of fluid based method for fluid flows with phase change, Journal of computational physics, 160 (2000) 662-682.

[8] M.W. Akhtar, S.J. Kleis, Boiling flow simulations on adaptive octree grids, International Journal of Multiphase Flow, 53 (2013) 88-99.

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

[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) 387-410.

[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:  AVT-178 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) 857-881.

[15] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Computers in physics, 12 (1998) 620-631.

[16] J.U. Brackbill, D.B. Kothe, C. Zemach, A Continuum Method for Modeling Surface Tension, J. Comput. Phys., 100 (1992) 335-354.

[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) 960-982.

[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) 11-28.

[19] N. Samkhaniani, A. Gharehbaghi, Z. Ahmadi, Numerical simulation of reaction injection molding with polyurethane foam, Journal of Cellular Plastics, 49 (2013) 405-421.

[20] W.H. Lee, Pressure iteration scheme for two-phase flow modeling, IN" MULTIPHASE TRANSPORT: FUNDAMENTALS, REACTOR SAFETY, APPLICATIONS". (1980) 407-432.

[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) 238-248.

[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) 170-185.

[23] H.G. Weller, G. TaboraI, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys., 12 (1998) 620-632.

[24] C. Rhie, W. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA journal, 21 (1983) 1525-1532.

[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 VOF-based 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 second-order scheme, Journal of Computational Physics, 14 (1974) 361-370.

[28] M. Darwish, F. Moukalled, TVD schemes for unstructured grids, International Journal of Heat and Mass Transfer, 46 (2003) 599-611.

[29] B. Shu, F. Dammel, P. Stephan, Phase change model for two-phase 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) 183-192.

[31] P. Berenson, Film-boiling heat transfer from a horizontal surface, Journal of Heat Transfer, 83 (1961) 351-356.

[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. V001T003A014-V001T003A014.

[37] C. Aghanajafi, K. Hesampour, Heat transfer analysis of a condensate flow by VOF method, Journal of Fusion Energy, 25 (2006) 219-223.

[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) 39-53.