A volume of fluid based method for vapor-liquid phase change simulation with numerical oscillation suppression

International Journal of Heat and Mass Transfer 110 (2017) 348–359

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

A volume of fluid based method for vapor-liquid phase change simulation with numerical oscillation suppression Shui-Ting Ding, Bin Luo, Guo Li ⇑ Aircraft/Engine Integrated System Safety Beijing Key Laboratory, School of Energy and Power Engineering, Beihang University, Beijing 100191, China

a r t i c l e

i n f o

Article history: Received 2 September 2016 Received in revised form 25 February 2017 Accepted 5 March 2017

Keywords: Phase change Volume of fluid (VOF) Piecewise linear interface calculation (PLIC) Numerical oscillation Energy source donor-acceptor scheme

a b s t r a c t A method of modeling phase change and suppressing numerical oscillation in vapor-liquid phase change simulation based on volume of fluid (VOF) is proposed in this study. The vapor-liquid interface is tracked by using VOF and reconstructed by applying a piecewise linear interface calculation (PLIC) scheme. Mass and heat transfer across the interface are modeled by employing a mass source term and an energy source term in the governing equations. In phase change model, the interface area needs to be calculated explicitly by assuming that the interface is linear in a cell. Expressions of two-dimensional interface length in the phase change model are also provided. The cause of the numerical oscillation, which may appear in the simulation, is analyzed. The analysis shows that small interfacial thermal resistance, large time step and ratio of interface area to cell volume are the causes of the numerical oscillation. An energy source donor-acceptor scheme is presented to suppress the numerical oscillation. One-dimensional Stefan problem, two-dimensional film boiling phenomenon and two-dimensional film condensation process were simulated using the model and scheme. Results validate the feasibility of using this method to simulate vapor-liquid flow with phase change. The numerical oscillations are also suppressed effectively in the simulations. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Boiling and condensation can effectively transfer large amounts of heat because of the large amounts of energy in the form of latent heat. The two phase change processes are very important in extracting energy from various sources, such as solar, fossil, and nuclear fuels. Hydrodynamic and heat transfer processes in boiling and condensation are linked very closely. This coupling is much closer than that in single-phase flows [1]. In addition, it is extremely difficult to conduct accurate and effective experimental measurements because the spatial scales are small and the time constants are rapid in the phase change process. The importance and complexity have led to decades of experimental and theoretical studies. Substantial experimental and analytical studies in the past have resulted in numerous empirical relations or correlations associated with phase change. Workers have successfully matched experimental data to the relations or correlations. However, empirical finding or correlation is only applicable to the specific conditions under which it was developed. Most of these findings or correlations cannot be applied to new situations because of their

⇑ Corresponding author. E-mail address: [email protected] (G. Li). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.03.015 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

limited range of applicability. Numerical methods for predicting phase change problems have become popular and reliable as a result of the rapid development of computer hardware and numerical algorithms [2]. Accurate numerical simulations can reveal details of heat and mass transfer processes that are not measurable in experiments. Such procedures will facilitate understanding of underlying physical behavior. Two main difficulties exist in accurate phase change simulations. One difficulty is the accurate tracking of vapor-liquid interfaces in phase change problems. Different interface tracking methods have been used to simulate the phase change problems. Son and Dhir [3] used a coordinate transformation technique supplemented by a numerical grid generation method to construct a grid system aligned with the interface. This method was used to simulate saturated film boiling on a horizontal surface. The same method was also applied by Banerjee and Dhir [4] in studying subcooled film boiling. Welch [5] used an interface tracking method in conjunction with a finite volume method on a moving unstructured mesh to simulate the axisymmetric vapor bubble growth. When the interface changes, these methods require remeshing to be adapted to the interface deformation. Thus, these methods show a limited applicability because of their inability to cope with large interfacial distortion or change in topology. This limitation was later overcome by Juric and Tryggvason [6]. They used front track-

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

349

Nomenclature Ai cp E ! F !CSF g Ga hfg k lcr lm _ m ! n Nu p Pr q QE r SE Sm t Dt T ! u Vc x; y

interface area in a cell ðm2 Þ specific heat ½J=ðkg KÞ energy per unit mass ðJ=kgÞ surface tension induced volume force ðN=m3 Þ gravity acceleration vector ðm=s2 Þ Galileo number latent heat ðJ=kgÞ phase change mass transfer coefficient ½kg=ðm2 s KÞ critical wavelength ðmÞ most unstable wavelength ðmÞ mass flux across the interface ½kg=ðm2 sÞ unit normal vector Nusselt number pressure ðPaÞ Prandtl number heat flux ðW=m2 Þ produced heat per second ðWÞ mass transfer intensity factor ðs1 Þ energy source term ðW=m3 Þ mass source term ½kg=ðm3 sÞ time ðsÞ time step ðsÞ temperature ðKÞ velocity vector ðm=sÞ volume of a cell ðm3 Þ Cartesian coordinates

ing method, which is based on finite difference method, to simulate boiling flow. A fixed grid was used in the method. A separate front marked the interface, and was only modified near the front to make a grid line follow the interface. This method was also used by Esmaeeli and Tryggvason [7] to simulate three-dimensional film boiling. A level set method can also overcome this limitation. The interface in the level-set method is captured and tracked by the level-set function, which is defined as a signed distance from the interface [8]. This method was used by Son and Dhir to simulated axisymmetric film boiling flows [8] and the bubble merger process on a single nucleation site during pool nucleate boiling [9]. Mukherjee and Dhir [10] used this method to simulate the lateral merger of vapor bubbles during nucleate pool boiling. Luo et al. [11] combined the variable density projection method with levelset method to simulate immiscible interfacial flows with phase change. Wu et al. [12] simulated subcooled nucleate boiling coupling level-set method with a moving-mesh method. Despite its extensive application, level-set method is found to have a deficiency in preserving volume conservation [13]. Volume of Fluid (VOF) method is another method that can model the vapor-liquid interface and it does not have the deficiency in preserving volume conservation. Volume fraction is employed in VOF to track the interface. Various technologies for simulating phase change problems have been used in VOF to obtain an accurate curvature and to smoothen the discontinuous physical quantities near the interfaces, for example, the PLIC scheme [13–22], the modified high resolution interface capturing (HRIC) scheme [23,24], and the conservative interpolation scheme for interface tracking (CISIT) [25]. The PLIC scheme is one of the most popular methods used in the applications. The other important difficulty is the need to model the mass and heat transfer of phase change at the interface. The mass and heat flux across the interface needs first to be modeled. Two approaches were commonly used to obtain the interfacial mass

Greek symbols volume fraction dx cell width ðmÞ d thickness ðmÞ k thermal conductivity ½W=ðm KÞ r surface tension ðN=mÞ j surface curvature q density ðkg=m3 Þ l dynamic viscosity ðPa sÞ mv kinematic viscosity ðm2 =sÞ

a

Subscripts and superscripts a acceptor c condensing d donor e evaporating i interface l liquid max maximum min minimum sat saturation v vapor w wall

flux in the existing references. One approach involves direct or indirect application of energy jump condition at the interface [5, 7–13,15,17,19,20,22,25,26]. The condition can be expressed as

_ l hfg ¼ q ¼ _ v hfg ¼ m m


 
  ! @T
@T
 n; kl

 kv

@n l @n v

ð1Þ

!

where n is a unit vector normal to the interface and pointing toward the vapor. In Eq. (1), q > 0 means evaporation and q < 0 means condensation. However, this method does not consider the interfacial thermal resistance for phase change. The other approach is based on the interfacial thermal resistance model proposed by Schrage [21,27–31]. In the model, the mass flux at the interface can be determined by

_l¼ _ v ¼ m m

2b 2b

rffiffiffiffiffiffiffiffiffi M qv hfg ðT  T sat Þ ; 2pR T 3=2 sat

ð2Þ

where coefficient b represents the fraction of molecules transferred from one phase to the other during phase change. Then, volumetric source terms should be derived from the flux model to apply it to the phase change simulations. Two commonly methods are available by which the source terms can be constructed in the references. The first method is using the feature of the volume faction [31]. The feature can be expressed as

Z

X

jrajdX ¼

Z

dS;

ð3Þ

S

where a can be either the liquid or vapor volume fraction. Thus, the integration represents the interface area in theory. After using this feature, the mass source terms in the first method can be written as [8,19,20,28,30,31]:



_ v jraj Sm;v ¼ m : _ v jraj Sm;l ¼ m

ð4Þ

350

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

That means jrajV c could represent the interface area in a cell, which contains an interface. The accuracy of the sources terms highly depends on the calculation of the gradient of the volume fraction. In the second method, which is only applicable to interfacial thermal resistance based flux model, the source terms are obtained based on the assumption that the flow regime is dispersed with a constant diameter [29,32], and they are known as Lee’s model [33]. They can be expressed as follows [16,18,21,23,24,29,32–34]:

Sm;v ¼ Sm;l

ðT  T sat Þ ¼ re ql al ; T sat

T > T sat ðevaporating processÞ; ð5Þ

Sm;l ¼ Sm;v ¼ rc qv av

ðT sat  TÞ ; T sat

T < T sat ðcondensing processÞ; ð6Þ

where r stands for the mass transfer intensity factor with unit s1 . The factor should theoretically be different for evaporation and condensation expression. Different values were empirically used to numerically maintain small difference between the interfacial temperature and the saturation temperature in simulating phase change problems. For example, Riva et al. [24] used factors ranging from r ¼ 7:5  105 s1 to r ¼ 5  106 s1 to maintain the interface temperature within approximately 1 K of the saturation temperature. Meanwhile, Yang et al. [34] only used r ¼ 100 s1 to achieve the same purpose. Using this model to simulate the phase change problems, which do not satisfy the assumption, may not be reasonable. To the best of authors’ knowledge, the method has not been verified by the basic problems of phase change, such as the Stephan problem and the film boiling problem that were widely used for verification. When interfacial thermal resistance based phase change models are applied to the phase change simulations, it may encounter some numerical problems in some cases. Hardt and Wondra [31] found that source terms located in a very narrow region may lead to numerical instabilities. Researchers who used Lee’s model found that too large value of r, which represents a small interfacial thermal resistance, would induce a numerical oscillation and it can cause the solution to diverge, whereas too small value lead to a significant deviation between the interfacial temperature and the saturation temperature [16,21,23,24,29,34]. However, the cause of the numerical oscillation is yet to be analyzed in these references. Another method of modeling the phase change is proposed in this study. In this method, the interface area needs to be calculated directly by assuming that the interface is linear in the cell. The cause of the numerical oscillation is analyzed. An energy source donor-acceptor scheme is presented to suppress the numerical oscillation. The method is tested in FLUENT 14.5 by using userdefined functions (UDFs). One-dimensional Stefan problem, twodimensional film boiling phenomenon, and two-dimensional film condensation process were simulated to validate the effectiveness of the method. 2. Mathematical formulation 2.1. Governing equations The volume fraction a of phase is introduced in VOF method. The value of a lies between 0 and 1. If the liquid’s volume fraction of the cell is denoted as al . al ¼ 1 indicates that the cell is full of liquid. al ¼ 0 indicates that the cell is full of vapor, and 0 < al < 1 indicates that the cell contains an interface between the vapor

and liquid. In each cell, the sum of the volume fractions of vapor and liquid is 1, as follows:

al þ av ¼ 1:

ð7Þ

The tracking of the interface between the liquid and vapor is accomplished by solving a continuity equation for the volume fraction of liquid or vapor. For the liquid and vapor phase, the equations are provided as ! @ ðal ql Þ þ r  ðal ql u Þ ¼ Sm;l ; @t

ð8Þ

! @ ðav qv Þ þ r  ðav qv u Þ ¼ Sm;v ; @t

ð9Þ

where S is the mass source, which will be discussed in detail in the next subsection. As a relation [Eq. (7)] exists between the volume fraction of liquid and vapor, one of the volume fraction equations will not be solved. The liquid’s volume fraction equation is chosen to be solved in this study. It is assumed that the flows in both phases are incompressible. Thus the mass conservation equation in VOF model can be written as follows: ! @q þ r  ðq u Þ ¼ 0; @t

ð10Þ

where

q ¼ al ql þ av qv :

ð11Þ

Surface tension is taken into account by using the continuum surface force model proposed by Brackbill et al. [35], which is given as: !

F CSF ¼ 2r

al ql jv rav þ av qv jl ral ; ql þ qv

ð12Þ

where r is the surface tension, and j is the surface curvature. Thus, the momentum conservation equation for VOF model is described by the following: ! ! ! ! ! ! ! @ ðq u Þ þ r  ðq u u Þ ¼ rp þ r  ½lðr u þr u T Þ þ q g þF CSF ; @t

ð13Þ

where

l ¼ al ll þ ð1  al Þlv :

ð14Þ

The energy equation for VOF model is given as follows: ! @ ðqEÞ þ r  ½u ðqE þ pÞ ¼ r  ðkrTÞ þ SE ; @t

ð15Þ

where

E¼

al ql El þ av qv Ev ; al ql þ av qv

ð16Þ

El ¼ cp;l ðT  298:15 KÞ;

ð17Þ

Ev ¼ cp;v ðT  298:15 KÞ;

ð18Þ

k ¼ al kl þ ð1  al Þkv ;

ð19Þ

SE ¼ hfg Sm;l :

ð20Þ

2.2. Phase change model A phase change model is presented in this subsection. The mass and heat transfer across the interface is modeled by employing a mass source term and an energy source term in the governing

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

equations. Interfacial thermal resistance is taken into account. However, disagreement exists on the formulation of the resistance. Different expressions, such as kinetic theory based model [27], statistical rate model [36] and non-equilibrium model [37], that describe the resistance have been proposed. Despite the disagreement, the heat flux across the interface can be simplified as follows:

_ v hfg ¼ m _ l hfg ¼ khfg ðT  T sat Þ; q¼m

ð21Þ

where k is a mass transfer coefficient for phase change with unit kg=ðm2 s KÞ, and 1=khfg represents the interfacial thermal resistance. The formulation of k depends on the selective resistance model and the value of k is always greater than 0. Eq. (21) shows that when the heat flux is constant, temperature difference jT  T sat j decreases as k increases, but the limit of the difference is 0. Thus, a coefficient kl should exist. The difference can be ignored in numerical simulations when k > kl . Therefore, the simulation results should slightly differ in the use of different values of k when k > kl . The results should be consistent with those obtained when interfacial thermal resistance is ignored. The mass source terms in Eqs. (8) and (9) can be expressed as follows:

Sm;v ¼ Sm;l ¼

kðT  T sat ÞAi ; Vc

ð22Þ

where Ai and V c are the interface area and cell volume respectively in three-dimensional case. They are the interface length and the cell area respectively in two-dimensional case. A brief introduction on how to calculate the interface area Ai is provided herein. For the purpose of illustration, two-dimensional case is considered. The interface is assumed to be linear in the interface cell in this study. Any interface cell can be rotated to one of the interface configurations shown in Fig. 1. The angle the interface makes with the x-axial could be expressed as

u ¼ arctan

  @ al =@x ðp < u 6 pÞ: @ al =@y

The angle



ð23Þ

v is defined as

dx tan u v ¼ arctan dy



ð0 6 v 6 p=2Þ:

ð24Þ

The liquid fractions of the top, right, bottom and left sides of the cell are defined as st , sr , sb and sl respectively. Then, the interface length in Fig. 1(a-d) can be respectively calculated by the following equations:

Ai;a

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ðsb dxÞ2 þ ðsr dyÞ2 ¼ ð2al cot vÞdx2 þ ð2al tan vÞdy2 ;

ð25Þ

Ai;b

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ dx2 þ ½ðsr  sl Þdy2 ¼ dx2 þ ðdy tan vÞ2 ;

ð26Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½ðsb  st Þdx2 þ dy2 ¼ ðdx cot vÞ2 þ dy2 ;

ð27Þ

Ai;c ¼

Ai;d

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ½ð1  st Þdx2 þ ½ð1  sl Þdy2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ½2ð1  al Þ cot vdx2 þ ½2ð1  al Þ tan vdy2 :

351

ð28Þ

2.3. Energy source donor-acceptor scheme This subsection discusses the cause of the numerical oscillation which may occur in the simulation. An energy source donoracceptor scheme is proposed to suppress the oscillation. Neglecting the advection and diffusion terms, Eq. (15) can be expanded as follows:

hfg kAi @ ½ðT  298:15 KÞðal ql cp;l þ av qv cp;v Þ ¼ ðT sat  TÞ: @t Vc

ð29Þ

The equation can be changed into the following form when the first-order explicit scheme is applied:

 n hfg kAi T nþ1  T n ¼ Dt; ðal ql cp;l þ av qv cp;v ÞV c T nsat  T n

ð30Þ

where n and n þ 1 represent the index for the previous and new time step respectively. If the right side of Eq. (30) is greater than 1, the following expressions are true:

(

T nþ1 > T nsat ; if T n < T nsat T nþ1 < T nsat ; if T n > T nsat

:

ð31Þ

Inequality (31) indicates that if the interface is condensing at the previous time step, it will change to evaporating at the new time step, and vice versa. It is a numerical oscillation phenomenon. Thus, numerical oscillation is attributed to the large value of the right side of Eq. (30). An analysis of Eq. (30) shows that numerical oscillation may easily occur in two cases. The first case is when too large values of k  Dt are used. The second case is illustrated in Fig. 2: The value of al is relatively small, and the value of Ai =V c is relatively large. As ql  cp;l is always higher than qv  cp;v for most materials, the second case may also cause numerical oscillation. To avoid the numerical oscillation, the right side of Eq. (30) should satisfy the following inequality:



hfg kAi ðal ql cp;l þ av qv cp;v ÞV c

n

Dt < 1:

ð32Þ

Then, Dt should satisfy the following inequality:

q cp;v dx Dt < v : hfg k

ð33Þ

Inequality (33) shows that a very small time step Dt must be used if no treatment is used. This procedure may lead to an expensive simulation. Thus, an energy source donor-acceptor scheme is proposed to reduce the requirement of the time step. In this scheme, the cell, wherein the value of al  av and SE is not zero, is defined as a source donor as shown in Fig. 3. The liquid in the source donor can only be in the shape of a triangle, quadrilateral

Fig. 1. Four possible interface configurations.

352

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

The maximum energy source in the donor should be expressed as follows:

SE;d;max ¼

Q E;max al;d ql cp;l ðT sat  TÞ ¼ : Vd Dt

ð39Þ

The energy source in the donor-acceptor scheme without an acceptor can be expressed as follows:



Fig. 2. Cell with small al and large Ai =V c .

SE;d ¼

SE ; jSE j < jSE;d;max j : SE;d;max ; jSE j P jSE;d;max j

ð40Þ

The use of Eqs. (35), (36), and (40) can lower the requirement of the time step, but the required time step may remain small enough to cause expansive simulation costs because the time step is constrained by inequality (37). To use a large time step that will not satisfy the inequality, the energy source in the donor is modified to guarantee that the numerical oscillation will not occur. The modified energy sources in the energy source donor-acceptor scheme are formulated as follows:



SE;d;lim ited ¼

SE;a;limited ¼

Fig. 3. Definition of source donor and acceptor. Liquid in the source donor can be in the shape of a (a) triangle, (b) quadrilateral or (c) pentagon.

SE;d ;

jSE;d j < jSE;d;max j

SE;d;max ; jSE;d j P jSE;d;max j

SE;d;limited

al;d

Ta  Td > 0; T d  T sat

ð34Þ

where T a and T d are the temperature of the acceptor and donor, respectively. Energy sources at the donor and acceptor in the energy source donor-acceptor scheme are evaluated as follows:

SE;d ¼

SE V d Pn

al;d V d þ

SE;a;i ¼

i¼1 V a;i



al;d V d Vd

SE V d Pn

al;d V d þ

¼

al;d V d P S ; al;d V d þ ni¼1 V a;i E

V d;i Vd SE;d Pn  ¼ SE ¼ ; V al;d V a V þ V d;i l;d d i¼1 a;i i¼1 a;i

ð35Þ

ð36Þ

where n is the total number of the acceptor, and i is the index of the acceptor. After using this scheme, the required time step Dt can be estimated as follows:

Dt <

ql cp;l dx hfg k

:

ð37Þ

When inequality (33) is compared with inequality (37), the required time step Dt in the scheme can be ql cp;l =qv cp;v times larger than that without any treatment. When an energy source donor has no source acceptor, the phase change energy can be provided only by the donor itself. The maximum energy that the donor can provide for the phase change in a time step can be evaluated by the following expression:

Q E;max ¼

al;d ql cp;l V d ðT sat  TÞ Dt

:

ð38Þ

ð41Þ

:

ð42Þ

The time step constraint can then be removed after using the modified sources. When time step is sufficiently small, Eq. (41) is the same as Eq. (35). Therefore, Eqs. (41) and (42) are recommended for the energy source donor-acceptor scheme. Then, the mass sources in the scheme can be expressed as follows:

Sm;l;limited ¼ Sm;v ;limited ¼ or pentagon. The source acceptor is defined according to the shape of the liquid in the source donor, as shown in Fig. 3. The source acceptor must be full of liquid and the following inequality must be satisfied

;

SE;d;limited V d þ

Pn

i¼1 SE;a;limited V a;i

Vd

:

ð43Þ

3. Verification of the method In this section, the method is applied to solve a onedimensional Stefan problem, a two-dimensional film boiling problem, and a two-dimensional film condensation problem. A userdefined code is developed to implement the energy source donor-acceptor scheme based on the FLUENT platform. The second-order upwind scheme is adopted for the convective terms in the momentum equation and energy equation. The PLIC georeconstruct method is applied to solve the volume fraction equation. The double-precision pressure based solver is selected to solve the governing equations. In the verification, different ranges of k are used in different cases. In a specific case, a mass transfer coefficient k1% , which is defined as

k1% ¼

q ; hfg jDT max j  1%

ð44Þ

should fall in the range. That means, when k ¼ k1% , the ratio of the difference between the interface temperature and saturation temperature jT i  T sat j to the maximum temperature difference jDT max j is within 1%. 3.1. One-dimensional Stefan problem The schematic of the one-dimensional Stefan problem can be seen in Fig. 4. The length of computation domain is 0.2 m. The vapor and liquid are assumed to be incompressible with constant properties. The vapor properties used are qv ¼ 0:001 kg=m3 , kv ¼ 0:005 W=ðm KÞ, and cp;v ¼ 200 J=ðkg KÞ. The liquid properties used are ql ¼ 1 kg=m3 , kl ¼ 0:6 W=ðm KÞ, and cp;l ¼ 400 J=ðkg KÞ. The latent heat is hfg ¼ 10 kJ=kg. The vapor and liquid are initially in quiescent equilibrium. The saturation temperature is

353

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

Fig. 4. Schematic of one-dimensional Stefan problem.

T sat ¼ 300 K, and the wall temperature is T w ¼ 310 K. The vapor experiences an increase in temperature on the solid boundary and becomes superheated. Heat is transferred from the wall to the interface through the vapor. The liquid at the interface is heated to achieve a temperature higher than the saturation temperature before evaporation. Thus, the generated vapor pushes the liquid away from the wall. The theoretical interface position or vapor thickness can be described as follows [38]:

sffiffiffiffiffiffiffiffiffiffiffiffiffi kv t dðtÞ ¼ 2e ; qv cp;v

e is a solution to the following transcendental equation:

e expðe2 Þerf ðeÞ ¼

cp;v ðT w  T sat Þ pffiffiffiffi : hfg p

ð46Þ

Different grids and time steps were used to ensure independent results, as shown in Fig. 5. In the independence test, the value of k is 100. As a result, using mesh number n ¼ 800 and time step Dt ¼ 0:0001 s is adequate in the simulation. Thus, the mesh length in the calculation is dx ¼ 2:5  104 m. At t ¼ 0, the vapor thickness was initialized to be dx to make sure that the initial heat flux at the interface can be calculated correctly. As can be seen from Eq. (45), the vapor length increases with time. Thus the phase change flux at the interface decreases with time. Therefore, the heat flux at t ¼ 0 at the interface is the maximum. It was used to estimated k1% . The heat flux was estimated by the following:

q ¼ kv

Fig. 5. Grid and time step independence test for Stephan problem simulation with k ¼ 100 kg=ðm2 s KÞ.

T sat  T w : dx

0.12 Theory k=0.0005 k=0.00075 k=0.001 k=0.002

0.10 0.08

δ/m

where

ð45Þ

k=0.005 k=0.01 k=0.1 k=1 k=10 k=100

0.06 0.04 0.02 0.00 0.0

0.2

0.4

0.6

0.8

1.0

t/s ð47Þ

The maximum temperature difference DT max is jT w  T sat j. Thus, k1% can be calculated from Eq. (44), and the value of k1% is 0.02 in this case. When k > k1% , the temperature difference jT i  T sat j can be less than 0.1 K. So using the range 0.0005–100 for k to simulate different interfacial thermal resistances in this problem is reasonable. Fig. 6 shows the simulation results for different values of k. When k 6 0:002 kg=ðm2 s KÞ, the interfacial thermal resistance decreases as k increases, thereby increasing the moving velocity of the interface. However, when k P 0:005 kg=ðm2 s KÞ, the result does not change as k increases, and it is highly consistent with the theoretical solution. It is because a large value of k indicates a small resistance. The resistance exhibits little influence to the result, and can be ignored when k P 0:005 kg=ðm2 s KÞ. In addition, the resistance is not considered in the theoretical solution. This law is highly consistent with the analysis of Eq. (21) in Section 2.2.

Fig. 6. Results of Stephan problem with different k ½kg=ðm2 s KÞ.

These results show that the method proposed in this study can simulate accurately the one-dimensional phase change with or without the interfacial thermal resistance. Table 1 shows the comparison of vapor thickness and minimal temperature after three time steps under different conditions. All the values of k in the table are greater than 0.005. Thus, the results in different cases should be the same because the interfacial thermal resistances are small enough. If k P 1 kg=ðm2 s KÞ and the energy source donor-acceptor scheme is not used, the minimum temperature would be smaller than the saturation temperature. As k increases, the temperature deviation increases. It is a temperature oscillation phenomenon, and it is unrealistic. The phenomenon can be explained by Eq. (30) and inequality (31). In addition, the vapor thickness increases as k increases. When

354

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

Table 1 Comparison of vapor thickness and minimal temperature after three time steps in Stefan problem simulations. Case

k ½kg=ðm2 s KÞ

Using the energy source donor-acceptor scheme?

Vapor thickness d (m)

T min (K)

1 2 3 4 5 6

0.1 1 10 0.1 1 10

Yes Yes Yes No No No

0.000875 0.000875 0.000875 0.000875 0.006503 0.062776

300 300 300 300 299.94 297.33

k ¼ 10 kg=ðm2 s KÞ, the vapor thickness is 0.06277 m, which should be reached after t ¼ 0:4 s in theory. However, when the energy source donor-acceptor scheme is used, the unreasonable phenomenon disappears and the vapor thicknesses are the same, although the value of k is different. The comparison shows the following: (1) numerical oscillation can easily occur in the Stefan problem simulation if the mass transfer coefficient is too large; (2) the numerical oscillation causes the results to be inaccurate; and (3) the energy source donor-acceptor scheme can suppress the numerical oscillation and guarantee the results to be accurate when the mass transfer coefficient is relatively large. 3.2. Two-dimensional film boiling Fig. 7 shows the schematic of the two-dimensional film boiling problem. A two-phase fluid with latent heat hfg ¼ 10 kJ=kg, surface tension r ¼ 0:1 N m, saturation temperature T sat ¼ 300 K, and saturation pressure P sat ¼ 101325Pa was considered. The vapor and liquid were assumed to be incompressible with constant properties as well. The vapor properties used were qv ¼ 5 kg=m3 , lv ¼ 0:005 Pa s, kv ¼ 1 W=ðm KÞ, and cp;v ¼ 200 J=ðkg KÞ. The liquid ql ¼ 200 kg=m3 , ll ¼ 0:1 Pa s, properties used were cp;l ¼ 400 J=ðkg KÞ, and kl ¼ 40 W=ðm KÞ. The gravity acceleration was g ¼ 9:81 m=s2 . Initially, the liquid was saturated, and the wall with a temperature of T w ¼ 305 K was superheated. The distance between the interface and the wall is provided by the following equation:

dðx; 0Þ ¼

   lm 2px ; 4:0 þ cos lm 128

ð48Þ

where lm is the most unstable Taylor wavelength that can be expressed as follows:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3r : lm ¼ 2p ðql  qv Þg

ð49Þ

The vapor temperature increased linearly from the interface to the wall. The velocity was initialized at 0. Because of the symmetry, a computational domain with a width of lm =2 and a height with lm was used. As the bubble arising from the wall can be considered as axisymmetric, an axisymmetric two-dimensional model was used to simulate the problem. A non-axisymmetric model was also simulated to study the difference between the method of this study and the method in which the energy jump condition was applied. The film boiling flow was considered to be laminar in the simulation. The average heat flux at the interface at t ¼ 0 was used to determine the value of k1% in either axisymmetric or non-axisymmetric model. In axisymmetric model, the average heat flux was estimated by

T w  T sat qaxi ¼ kv R lm =2 : 0

2pxdðx;0Þdx

ð50Þ

2

pðlm =2Þ

In non-axisymmetric model, the average heat flux was estimated by

T w  T sat qnonaxi ¼ kv R lm =2 : 0

dðx;0Þdx

ð51Þ

lm =2

Applying Eq. (44), the values of k1% in axisymmetric and nonaxisymmetric models are 10 and 8 respectively. When the value of k is greater than 10 in either model, the temperature difference jT i  T sat j could be less than 0.05 K most of the time. Thus, using the range 0.1–1000 for k to simulate different interfacial thermal resistances in this problem is reasonable. A convergence study about grid resolution and time step was conducted in the axisymmetric model. Three grids of resolution, namely, 32ðxÞ  64ðyÞ, 64ðxÞ  128ðyÞ, and 128ðxÞ  256ðyÞ, were considered in the study, as shown in Fig. 8(a). Three time steps, i.e., Dt ¼ 1  105 s, Dt ¼ 2:5  105 s, and Dt ¼ 5  105 s, were also considered in the study as shown in Fig. 8(b). The ratio of the vapor void to the initial vapor void av =av ;0 was used to determine convergence. The results in the figures show that using grid

Fig. 7. Schematic of the two-dimensional film boiling problem.

of resolution 64ðxÞ  128ðyÞ with time step Dt ¼ 2:5  105 s is adequate in simulating the case under consideration. Different interfacial thermal resistances were considered by assigning different values to k. Fig. 9 shows the trends of av =av ;0 over time with different values of k. In axisymmetric model, when k 6 10 kg=ðm2 s KÞ, interfacial thermal resistance decreased the velocity of the vapor production, which increases with the increase in k value. The result did not change with change in k when k P 100 kg=ðm2 s KÞ. A similar law can also be observed in nonaxisymmetric model. However, the threshold values of k were different. The laws also agree well with the analysis of Eq. (21) in Section 2.2. Fig. 9(b) illustrates that the trend in non-axisymmetric model, when k P 100 kg=ðm2 s KÞ and t 6 0:35 s, is consistent with

355

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

7

5

αv/αv,0

k=10 k=100 k=1000

k=0.1 k=0.4 k=0.7 k=1

6

4 3 2 1 0

(a)

Fig. 8. Grid and time step independence test for two-dimensional film boiling problem with k ¼ 100 kg=ðm2 s KÞ.

the trend in Ref. [26] wherein the energy jump condition was used. Only a slight difference existed when t P 0:35 s. Fig. 10 shows the shape of vapor-liquid interface at different times with k ¼ 100 kg=ðm2 s KÞ when the interfacial thermal resistance was sufficiently small and negligible. A mushroom-shaped vapor bubble growth from the vapor layer can be observed in both axisymmetric and non-axisymmetric models. This phenomenon has been experimentally observed in situations with a high heat flux [39,40]. Klimenko [41] obtained several correlations to predict heat transfer in film boiling on a horizontal plate. These correlations were the best for comparing the results of the experiments on the film boiling of nine different liquids with an accuracy of 25%. The correlations are expressed as follows: ( 1=3 0:19½Gaðql =qv  1Þ Pr 1=3 Gaðql =qv  1Þ < 108 ðLaminarÞ v f1 Nu ¼ ; 1=2 1=3 0:0086½Gaðql =qv  1Þ Pr v f 2 Gaðql =qv  1Þ > 108 ðTurbulentÞ ð52Þ

where

f1 ¼

f2 ¼

8 > <1

h i1=3 hfg > : 0:89 cp;v ðT w T sat Þ 8 > <1

h i1=2 hfg > : 0:71 cp;v ðT w T sat Þ

hfg cp;v ðT w T sat Þ

6 1:4

hfg cp;v ðT w T sat Þ

> 1:4

hfg cp;v ðT w T sat Þ

6 2:0

hfg cp;v ðT w T sat Þ

> 2:0

;

ð53Þ

:

ð54Þ

In Eq. (52), Ga is the Galileo number and lcr is the critical wavelength. They are expressed as follows:

0.0

0.1

0.2

0.3

0.4

0.5

t/s

Fig. 9. Trends of av =av ;0 over time with different k ½kg=ðm2 s KÞ in (a) axisymmetric and (b) non-axisymmetric models. 3

Ga ¼

glcr 2m2v

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r lcr ¼ 2p ðql  qv Þg

ð55Þ

ð56Þ

The Nusselt number of the axisymmetric model is used for comparison with Klimenko’s correlations. The Nusselt number in the simulation is computed as follows:

Nu ¼

qav g;w lcr  ; ðT w  T sat Þ kv

ð57Þ

where qav g;w is the area-weight average heat flux at the wall. In this case, the value of Gaðql =qv  1Þ is 17,934. As can be known from Eq. (52), the film boiling flow in this study is laminar. The Nusselt number calculated by Klimenko’s correlations is 9.54. Fig. 11 shows the comparison between the results obtained from Klimenko’s correlations and the results from the simulation. The simulation result oscillates with time, but is close to Klimenko’s results. Simulation method in current study is considered reasonable. Fig. 12 shows the time-weight average Nusselt numbers with different phase change mass transfer coefficients k. The figure shows that only results with k P 10 kg=ðm2 s KÞ fall within the error band of Klimenko’s result. Therefore, the actual interfacial thermal resistance in the experiments conducted by Klimenko should be miniscule. Interfacial thermal resistance increases as k decreases when k 6 10 kg=ðm2 s KÞ, but the Nusselt number of the film boiling increases. This phenomenon occurs because large interfacial thermal resistance reduces the velocity of vapor production under the same temperature difference between the wall and the liquid.

356

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

Fig. 10. Vapor-liquid interface shapes at different times in (a) axisymmetric and (b) non-axisymmetric models.

Fig. 11. Comparison of Nusselt number with k ¼ 100 kg=ðm2 s KÞ.

20

Simulation Klimenko

18 16

Nuavg

14

+25% line

12 10 8

-25% line

6 4 2 0

0.1

1

10

100

1000

k/(kg•m •s •K ) -2

-1

-1

Fig. 12. Comparison of Nusselt number with different k ½kg=ðm2 s KÞ.

Fig. 13 shows the comparison of minimum temperature with different k when 0:2 s 6 t 6 0:4 s. When k P 10 kg=ðm2 s KÞ and the energy source donor-acceptor scheme is not used, the minimum temperature oscillates with time and is always smaller than T sat . This phenomenon is unrealistic. The increasing of k intensifies the temperature oscillation. The phenomenon can also be explained by Eq. (30) and inequality (31). When the scheme is used or k 6 1 kg=ðm2 s KÞ, the oscillation disappears and the minimum temperature would not smaller than T sat . The results indicate the

Fig. 13. Comparison of the minimum temperature between using and not using the energy source donor-acceptor scheme in film boiling simulations.

following: (1) numerical oscillation can easily occur in the film boiling problem simulation if the mass transfer coefficient is too large; and (2) the energy source donor-acceptor scheme can suppress the numerical oscillation and guarantee the results to be realistic.

357

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

3.3. Two-dimensional film condensation

2.0

105 Pa s, kv ¼ 0:023 W=ðm KÞ, and cp;v ¼ 1892 J=ðkg KÞ. The liquid

properties used are ql ¼ 972 kg=m3 , ll ¼ 3:6  104 Pa s, cp;l ¼ 4182 J=ðkg KÞ, and kl ¼ 0:668 W=ðm KÞ. The gravity acceleration was g ¼ 9:81 m=s2 . The vapor was assumed to be at a saturation temperature of 363 K and the temperature of the cold wall was set to 340 K. The inlet vapor velocity was assumed to be 0:1 m=s. Film thickness can be estimated using Nusselt theory [42]. Nusselt theory states that the thickness can be estimated using the following equation:

 1=4 4ll kl ðT sat  T wall Þx d¼ : g q2l hfg

ð58Þ

The wall heat flux can be estimated as follows:

q¼

kl ðT sat  T wall Þ3=4 :  1=4

ð59Þ

4ll kl x g q2l hfg

The average wall heat flux was used to determine the value of k1% in this model. The average wall heat flux was estimated by

Rl qav g ¼

0

qdx kl ðT sat  T wall Þ3=4 4 1=4 : l ¼  1=4 3 l 4l k

ð60Þ

l l g q2l hfg

According to the Eq. (44), the value of k1% is 0.5. When the value of k is greater than 0.5 in this case, the temperature difference jT i  T sat j can be less than 0.23 K most of the time. Thus, using different values in the range of 0.0005–1 for k to simulate different interfacial thermal resistances in this problem is reasonable. According to Eq. (58), the maximum film thickness can reach 7:13  105 m. Thus, in the area of y 6 1:2  104 m, mesh numbers ny2 of 12, 24, 36 and 48 in the y-direction were used to conduct a grid independence test. In the area of y > 1:2  104 m, 20 rows of mesh were used in the y-direction. One hundred columns of mesh were used in the x-direction in the entire area. The variable time step scheme, which guarantees a global Courant number less than 1, was used to set the size of time step. Fig. 15 shows the result of grid independence test. The liquid void histories of different meshing strategies were used to evaluate the grid independence. Fig. 15 shows that ny2 ¼ 36 is suitable for the film condensation simulation and that the fluid can be considered as being in quasi-equilibrium when t > 0:4 s. Different phase change mass transfer coefficients k were used to simulate the film condensation problem. In the simulation, the liquid film thickness and wall heat flux profiles changed with time. Thus, the time-weight average film thickness and wall heat flux

1.6

αl ×106

The film condensation on a vertical plate was simulated as shown in Fig. 14. A two-dimensional geometry with a length of 0:025 m and a width of 0:0025 m was used. A two-phase fluid with latent heat hfg ¼ 2309 kJ=kg was considered. The vapor and liquid were assumed to be incompressible with constant properties as well. The vapor properties used were qv ¼ 29 kg=m3 , lv ¼ 1:19

1.2

ny2=12 ny2=24 ny2=36 ny2=48

0.8 0.4 0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

t/s Fig. 15. Grid independence test for two-dimensional film condensation simulation with k ¼ 0:5 kg=ðm2 s KÞ.

profiles were used to compare with Nusselt theory, as shown in Fig. 16. When k 6 0:1 kg=ðm2 s KÞ, the film thickness and heat flux increase as the coefficient increases. The lager the coefficient is, the smaller the thermal resistance between the saturated vapor and wall is, the larger the heat flux will be, the larger the amount of condensation will be produced, and then the larger the liquid film thickness will be. However, film thickness and heat flux nearly remain unchanged as the coefficient increases when k P 0:5 kg=ðm2 s KÞ. This result indicates that the interfacial thermal resistance can be ignored when k P 0:5 kg=ðm2 s KÞ in the simulation. Under this condition, the result is in good agreement with Nusselt theory. However a slight difference exists near the inlet and outlet. This difference is attributed to the maximum vapor velocity at the inlet during the simulation. By contrast, the vapor is assumed to be stable in Nusselt theory. The film near the inlet is thinned by the vapor flow thereby reducing the total thermal resistance between the vapor and wall. Heat flux near the inlet is larger than that obtained from the theory. The total heat transferred to the wall is increased. Therefore, total condensation is increased as well. Thus, the film near the outlet is thicker than that obtained from the theory. Fig. 17 shows the comparison of the maximum temperature between using and not using the energy source donor-acceptor scheme. In the comparison, the initial conditions and time steps used for different numerical calculations were the same. The calculation results at t ¼ 0:8 s and k ¼ 0:5 kg=ðm2 s KÞ, which can be obtained from the grid independence test, were used as the initial condition. The time step Dt was 2  105 s, and it could ensure the global Courant number to be less than 1. The results in 0.008 s were used for discussion. Fig. 17(b) shows that, when k P 0:05 kg=ðm2 s KÞ and the energy source donor-acceptor scheme is not used, the maximum temperature is greater than the saturation temperature all the time and it oscillates with time. The oscillation intensifies as k increases. This phenomenon can also be explained by Eq. (30) and inequality (31). When the scheme is used, the oscillation can be suppressed as can be seen from Fig. 17(a). Although the oscillation still exists when k ¼ 0:5 kg=ðm2 s KÞ, it is much weaker than the former. This com-

Fig. 14. Schematic diagram of two-dimensional film condensation problem.

358

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359

δ×10 5/m

10 9 8 7 6 5 4 3 2 1 0 0.000

(a)

parison also shows that the energy source donor-acceptor scheme can suppress the numerical oscillation in the phase change simulation.

k=0.0005 k=0.001 k=0.005 k=0.1 k=0.5 k=1 Theory

0.005

4. Conclusion

0.010

0.015

0.025

x/m

700

k=0.0005 k=0.001 k=0.005 k=0.1 k=0.5 k=1 Theory

600 500

q/(kw•m-2)

0.020

This study proposed a VOF based method for vapor-liquid phase change simulation with numerical oscillation suppression. In this method, the interface area need to be calculated directly by assuming that the interface is linear in the cell. The cause of the numerical oscillation is analyzed and the energy source donor-acceptor scheme is proposed to suppress the oscillation. One-dimensional Stephan problem, two-dimensional film boiling and film condensation were simulated by using different phase change mass transfer coefficients. The key conclusions of the study are as follows:

400 300 200 100 0

(b)

0.000

0.005

0.010

0.015

0.020

0.025

x/m

Fig. 16. Time-weight average (a) film thickness and (b) wall heat flux profiles for different k ½kg=ðm2 s KÞ in film condensation simulations.

363.8

Tmax (K)

363.6 363.4

Acknowledgments

k=0.005, the scheme is used. k=0.05, the scheme is used. k=0.5, the scheme is used.

The work was supported by Innovation Team of Complex System Safety and Airworthiness of Aeroengine from Co-Innovation Center for Advanced Aeroengine of China. The work was funded by the National Natural Science Foundation of China (Grant No. 51676009).

363.2 363.0 362.8 0.000

References 0.002

(a)

0.004

0.006

0.008

t (s)

Tmax (K)

420 405

k=0.005, the scheme is not used. k=0.05, the scheme is not used. k=0.5, the scheme is not used.

390 375 360 0.000

(b)

(1) The results with a relatively large mass transfer coefficient for phase change were highly consistent with theoretical solutions. (2) Large mass transfer coefficient, time step and ratio of interface area to cell volume are the main causes of the numerical oscillation. (3) The numerical oscillation intensifies as the mass transfer coefficient increases if the energy source donor-acceptor scheme is not used. It can also reduce the accuracy of the numerical results. The energy source donor-acceptor scheme can suppress the oscillation effectively. (4) The method proposed in this study is highly suitable for modeling phase change problems with or without interfacial thermal resistance. Moreover, it would not easily encounter a numerical oscillation problem.

0.002

0.004

0.006

0.008

t (s)

Fig. 17. Comparison of the maximum temperature between using and not using the energy source donor-acceptor scheme in film condensation simulations.

[1] J.G. Collier, J.R. Thome, Convective Boiling and Condensation, third ed., Oxford University Press, 1994. [2] V.K. Dhir, Mechanistic prediction of nucleate boiling heat transfer-achievable or a hopeless task, J. Heat Transf. 128 (1) (2006) 1–12, http://dx.doi.org/ 10.1115/1.2136366. [3] V.K. Dhir, G. Son, Numerical simulation of saturated film boiling on a horizontal surface, J. Heat Transf. 119 (3) (1997) 525–533, http://dx.doi.org/ 10.1115/1.2824132. [4] D. Banerjee, V.K. Dhir, Study of subcooled film boiling on a horizontal disc: Part I-analysis, J. Heat Transf. 123 (2) (2001) 271–284, http://dx.doi.org/10.1115/ 1.1345889. [5] S.W.J. Welch, Direct simulation of vapor bubble growth, Int. J. Heat Mass Transf. 41 (12) (1998) 1655–1666, http://dx.doi.org/10.1016/S0017-9310(97) 00285-8. [6] D. Juric, G. Tryggvason, Computations of boiling flows, Int. J. Multiph. Flow 24 (3) (2010) 387–410, http://dx.doi.org/10.1016/S0301-9322(97)00050-5. [7] A. Esmaeeli, G. Tryggvason, Computations of Film Boiling. Part I: numerical method, Int. J. Heat Mass Transf. 47 (25) (2004) 5451–5461, http://dx.doi.org/ 10.1016/j.ijheatmasstransfer.2004.07.027. [8] G. Son, V. Dir, Numerical simulation of film boiling near critical pressures with a level set method, J. Heat Transf. 120 (1) (1998) 183–192, http://dx.doi.org/ 10.1115/1.2830042. [9] G. Son, N. Ramanujapu, V.K. Dhir, Numerical simulation of bubble merger process on a single nucleation site during pool nucleate boiling, J. Heat Transf. 124 (1) (2002) 51–62, http://dx.doi.org/10.1115/1.1420713.

S.-T. Ding et al. / International Journal of Heat and Mass Transfer 110 (2017) 348–359 [10] V.K. Dhir, A. Mukherjee, Study of lateral merger of vapor bubbles during nucleate pool boiling, J. Heat Transf. 126 (6) (2004) 1023–1039, http://dx.doi. org/10.1115/1.1834614. [11] X. Luo, M. Ni, A. Ying, M.A. Abdou, Numerical modeling for multiphase incompressible flow with phase change, Numer. Heat Transf. Fund. 48 (5) (2005) 425–444, http://dx.doi.org/10.1080/10407790500274364. [12] J. Wu, V.K. Dhir, J. Qian, Numerical simulation of subcooled nucleate boiling by coupling level-set method with moving-mesh method, Numer. Heat Transf. Fund. 51 (6) (2007) 535–563, http://dx.doi.org/10.1080/10407790601177763. [13] S.W.J. Welch, J.A. Wilson, Volume of fluid based method for fluid flows with phase change, J. Comput. Phys. 160 (2) (2000) 662–682, http://dx.doi.org/ 10.1006/jcph.2000.6481. [14] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, Numer. Meth. Fluid Dynam. (1982) 273–285. [15] M.H. Yuan, Y.H. Yang, T.S. Li, Z.H. Hu, Numerical simulation of film boiling on a sphere with a volume of fluid interface tracking method, Int. J. Heat Mass Transf. 51 (7–8) (2008) 1646–1657, http://dx.doi.org/10.1016/j.ijheatmass transfer.2007.07.037. [16] C. Fang, M. David, A. Rogacs, K. Goodson, Volume of fluid simulation of boiling two-phase flow in a vapor-venting microchannel, Front. Heat Mass Transf. 1 (1) (2010) 1–11, http://dx.doi.org/10.5098/hmt.v1.1.3002. [17] R.S. Maurya, S.V. Diwakar, T. Sundararajan, S.K. Das, Numerical investigation of evaporation in the developing region of laminar falling film flow under constant wall heat flux conditions, Numer. Heat Transf., Part A 58 (1) (2010) 41–64, http://dx.doi.org/10.1080/10407782.2010.490174. [18] D. Xu, T. Chen, Y. Xuan, Thermo-hydrodynamics analysis of vapor-liquid twophase flow in the flat-plate pulsating heat pipe, Int. Commun. Heat Mass Transf. 39 (4) (2012) 504–508, http://dx.doi.org/10.1016/j.icheatmasstransfer. 2012.02.002. [19] D. Sun, J. Xu, L. Wang, Development of a vapor-liquid phase change model for volume-of-fluid method in FLUENT, Int. Commun. Heat Mass Transf. 39 (8) (2012) 1101–1106, http://dx.doi.org/10.1016/j.icheatmasstransfer. 2012.07.020. [20] H. Ganapathy, A. Shooshtari, K. Choo, S. Dessiatoun, M. Alshehhi, M. Ohadi, Volume of fluid-based numerical modeling of condensation heat transfer and fluid flow characteristics in microchannels, Int. J. Heat Mass Transf. 65 (5) (2013) 62–72, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.05.044. [21] Z. Liu, B. Sunden, J. Yuan, VOF modeling and analysis of filmwise condensation between vertical parallel plates, Heat Transf. Res. 43 (1) (2012) 47–68, http:// dx.doi.org/10.1615/HeatTransRes. 2012004376. [22] A.A. Ganguli, A.B. Pandit, J.B. Joshi, Bubble dynamics of a single condensing vapor bubble from vertically heated wall in subcooled pool boiling system: experimental measurements and CFD simulations, Int. J. Chem. Eng. 2012 (1) (2012) 1–11, http://dx.doi.org/10.1155/2012/712986. [23] E.D. Riva, D.D. Col, Numerical simulation of laminar liquid film condensation in a horizontal circular minichannel, J. Heat Transf. 134 (5) (2012) 807–824, http://dx.doi.org/10.1115/1.4005710. [24] E.D. Riva, D.D. Col, S.V. Garimella, A. Gavallini, The importance of turbulence during condensation in a horizontal circular minichannel, Int. J. Heat Mass Transf. 55 (13–14) (2012) 3470–3481, http://dx.doi.org/10.1016/j.ijheat masstransfer.2012.02.026.

359

[25] Y. Tsui, S. Lin, Y. Lai, F. Wu, Phase change calculations for film boiling flows, Int. J. Heat Mass Transf. 70 (3) (2014) 745–757, http://dx.doi.org/10.1016/j. ijheatmasstransfer.2013.11.061. [26] 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, Numer. Heat Transf. Appl. 59 (11) (2011) 857–881, http://dx.doi.org/10.1080/10407782. 2011.561079. [27] R.W. Schrage, A Theoretical Study of Interface Mass Transfer, Columbia University Press, 1953. [28] N. Samkhaniani, M.R. Ansari, Numerical simulation of bubble condensation using CF-VOF, Prog. Nucl. Energy 89 (2016) 120–131, http://dx.doi.org/ 10.1016/j.pnucene.2016.02.004. [29] M. Bahreini, A. Ramiar, A.A. Ranjbar, Numerical simulation of bubble behavior in subcooled flow boiling under velocity and temperature gradient, Nucl. Eng. Des. 293 (2015) 238–248, http://dx.doi.org/10.1016/j.nucengdes.2015.08.004. [30] Q. Zeng, J. Cai, H. Yin, et al., Numerical simulation of single bubble condensation in subcooled flow using OpenFOAM, Prog. Nucl. Energy 83 (2015) 336–346, http://dx.doi.org/10.1016/j.pnucene.2015.04.011. [31] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source terms, J. Comput. Phys. 227 (11) (2008) 5871–5895, http://dx.doi.org/10.1016/j.jcp.2008.02.020. [32] ANSYS FLUENT 14.5, Theory Guide, 2012. [33] Lee WH. A pressure iteration scheme for two-phase modeling. In: Multiphase transport fundamentals, reactor safety, applications, vol. 1; 1980. [34] Z. Yang, X. Peng, P. Ye, Numerical and experimental investigation of two phase flow during boiling in a coiled tube, Int. J. Heat Mass Transf. 51 (5–6) (2008) 1003–1016, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2007.05.025. [35] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354, http://dx.doi.org/ 10.1016/0021-9991(92)90240-Y. [36] C.A. Ward, G. Fang, Expression for predicting liquid evaporation flux: statistical rate theory approach, Phys. Rev. E 59 (1) (1999) 429–440, http://dx.doi.org/ 10.1103/PhysRevE.59.429. [37] D. Bedeaux, S. Kjelstrup, Transfer coefficients for evaporation, Physica A 270 (3–4) (1999) 413–426, http://dx.doi.org/10.1016/S0378-4371(99)00162-4. [38] V. Alexiades, A.D. Solomon, Mathematical modeling of melting and freezing processes, J. Sol. Energy Eng. 115 (2) (1993). [39] J.S. Ervin, H. Merte, R.B. Keller, K. Kirk, Transient pool boiling in microgravity, Int. J. Heat Mass Transf. 35 (3) (1992) 659–674, http://dx.doi.org/10.1016/ 0017-9310(92)90125-C. [40] O. Kunito, K. Yoshiyuki, I. Akira, A. Shigebumi, Transient boiling heat transfer characteristics of R113 at large stepwise power generation, Int. J. Heat Mass Transf. 31 (10) (1988) 2161–2174, http://dx.doi.org/10.1016/0017-9310(88) 90125-1. [41] V.V. Klimenko, Film boiling on a horizontal plate-new correlation, Int. J. Heat Mass Transf. 24 (1) (1981) 69–79, http://dx.doi.org/10.1016/0017-9310(81) 90094-6. [42] W. Nusselt, Die OberflächenKondensation des Wasserdamfes, Z. Vereins deutscher Ininuere 60 (1916) 541–575.