Guo Et Al 2014 A Review of Computational Modelling of Flow Boiling in Microchannels

79
A Review of Computational Modelling of

Flow Boiling in Microchannels
Z. Guo, D.F. Fletcher and B.S. Haynes
School of Chemical and Biomolecular Engineering
The University of Sydney,
NSW 2006, Australia
Recieved: 3 December 2013; Accepted: 1 March 2014
Abstract
Flow boiling in microchannels has received enormous interest in academic research
over the past few decades because of its importance in the thermal management of
micro-structured devices. A vast number of experimental and theoretical studies have
been reported and an increasing number of computational studies of microchannel flow
boiling have been performed in the past few years. This article provides a review of the
previous studies of flow boiling in microchannels, with a focus on the computational
work. We present the governing equations, boundary conditions, numerical methods
and the mathematical treatment of phase change used. Some important numerical
studies are reviewed in detail and their strengths and weaknesses are presented. Finally,
a summary of the current status of modelling in this area is provided.
1. INTRODUCTION
The application of micro-structured devices is widespread in cutting edge science and technologies,
for examples in electrical engineering, medical engineering, pharmaceutical engineering and in the
micro-structured units used in modern chemical plants [1, 2]. Development of micro-structured
devices for engineering applications has presented major challenges for the understanding of heat
transfer in microchannels [3]. Flow boiling has received increasing interest because of its potential
for addressing the thermal management challenges arising from the increase in heat flux density in
micro-structured devices, especially in microelectronic devices [4, 5]. In pool and flow boiling in
macrochannels, bubbles are observed to have diameters that are larger than the typical diameters of
a microchannel, which indicates that the flow regime in microchannel flow boiling is likely to be
very different from the boiling flow regime at the macro-scale. Therefore, the correlations
developed for flow boiling at the macro-scale are not applicable for microchannel flow boiling. In
the past few decades, many experimental, theoretical and numerical studies have been performed
to examine flow boiling in microchannels. This article reviews the previous research work on the
Computational Fluid Dynamics (CFD) simulation of such boiling flows.
In this review, Section 2 describes the experimental studies and theoretical modelling of flow
boiling in microchannels, referencing some excellent reviews covering the experimental and
theoretical work. Section 3 presents the CFD framework used for flow boiling in microchannels,
including the governing equations, boundary conditions, numerical approaches and the phase
change models. In Section 4, numerical work on microchannel flow boiling performed by various
research groups is examined in detail. Finally, Section 5 presents some conclusions for
microchannel flow boiling research and gives perspectives for future work.
2. EXPERIMENTAL AND ANALYTICAL WORK ON MICROCHANNEL FLOW BOILING

Over the past few decades, a vast number of experimental studies have been performed to examine
flow boiling in microchannels. At the same time, analytical studies have led to theoretical models
for prediction of boiling flow behaviour in microchannels. Comprehensive reviews of
microchannel flow boiling are available in Thome [6], Garimella and Sobhan [7], Bertsch et al. [5],
Kandlikar [8] and Baldassari and Marengo [9], where many datasets from experimental studies are
compared with the available theoretical models. However, not only do the experimental results
*Corresponding author: E-mail: david.fletcher@sydney.edu.au

80 A Review of Computational Modelling of Flow Boiling in Microchannels
show contradictory outcomes, but also the analytical correlations do not provide satisfactory
predictions. Thus, the discussion of the dominant heat transfer mechanism of flow boiling in
microchannels continues.
2.1 Discussion of the dominant mechanism of flow boiling heat transfer

Many researchers [10-12] have noted that nucleate boiling, which is mainly influenced by the heat
flux at the macro-scale, was apparently the dominant mechanism for flow boiling heat transfer in
microchannels. The main reason for this belief is that the measurements from their experiments
showed that the heat transfer coefficient was strongly dependent on heat flux and weakly dependent
on mass flow rate and vapour quality [11-13]. Some pool boiling correlations, such as Armstrong
[14] and Cooper [15] have been used to predict flow boiling heat transfer in microchannels with
some success [8]. Bertsch et al. [5] compared twenty five predicted correlations with ten sets of
experimental data and found Cooper’s correlation to provide the overall best prediction. However,
although Cooper’s correlation was shown to be the best among the twenty five tested, it predicted
only 48% of the datasets within a deviation of 30%. This result indicates that nucleate boiling, as
assumed in Cooper’s correlation, may not be the right model for heat transfer in flow boiling at the
micro-scale or another possibility is that there are mechanisms other than nucleate boiling having
an important impact on flow boiling heat transfer in microchannels.
In fact, nucleate boiling only occurs across a small quality range. Therefore, if the heat flux is
high there is only a small zone of nucleate boiling before the flow regime changes to Taylor and
then annular flow. Experiments showed that the two-phase flow regime in microchannels
transitioned from bubbly to slug, slug-annular and annular flow at very low vapour qualities
(x < 0.1) [16]. Data from Lee and Mudawar [17] showed that the heat transfer was still strong in
the annular regime (x > 0.1), where very few bubbles were observed and the heat transfer
coefficient had not decreased much compared with that in the bubbly region. The same observation
has been made in many other researchers’ works [10, 18, 19], see, for example, Figure 1 taken from
Bao et al. [10]. Particularly, data from Bertsch et al. [19], Lin et al. [20] and Yun et al. [21] yielded
the peak heat transfer coefficients at vapour qualities in the range between 0.1 and 0.6, in which the
two-phase flow regime is most-likely annular. This fact suggests that nucleation might not be
the dominant mechanism for flow boiling in microchannels.
Other conflicting trends reported in many works also caused concerns as to the dominant role of
nucleation for flow boiling in microchannels. Lee and Lee [22] performed flow boiling experiments
with refrigerant R113 using rectangular low-aspect-ratio channels. The heat transfer coefficient
increased with mass flux and quality; however, the effect of heat flux appeared to be minor, which
was opposite to that expected if nucleate boiling were the dominant mechanism. These authors [22]
created a film flow heat transfer correlation and found it applicable to their data in the low flow rate
range. Their high flow rate heat transfer data appeared to be in accord with pool boiling
correlations. Another model given by Qu and Mudawar [23] described flow boiling at the
6
hexp (kWm−2 K−1)
4 HCFC123
446
335
2 279 kg.m−2.s−1
223
167
0
−0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0
xth
Figure 1. Heat transfer coefficient versus vapour quality for a heat flux of 39 kW m-2.
The heat transfer coefficient does not decrease significantly with increasing quality as
the flow becomes annular at around x = 0.1. Taken from Bao et al. [1].
Journal of Computational Multiphase Flows

Z. Guo, D.F. Fletcher and B.S. Haynes 81
micro-scale as a steady annular flow with droplets entrained in the vapour core. This model
predicted the heat transfer coefficient assuming all the heat was used to evaporate liquid at the gas-
liquid interface, which completely ignored the effect of nucleation in the liquid film. All the data
provided by the author from their own experiments [24] were located within a ± 40% error band
with a mean absolute error (MAE) of 13.3%. However, the physical picture of a gas core containing
droplets is not one that has been reported in any other study of boiling in micro-channels.
2.2 Effects of upstream compressibility

Recently, an experimental study conducted by Yang et al. [25] has demonstrated that inlet
compressibility has a significant effect on microchannel flow boiling. Their data showed that with
upstream compressibility (termed a “soft inlet”) the flow boiling in microchannels appeared
“nucleate-like” across the entire vapour quality range. However, an incompressible inlet condition
(“hard inlet”) led to convective boiling heat transfer behaviour except at the very low quality
condition, which still showed the expected nucleate boiling behaviour. This observation is very
important in that it pointed out that the massive datasets from previous experiments have been
produced with inlets with varying degree of compressibility, arising from the amount of trapped gas
or vapour, within the experimental apparatus which had not been controlled nor accounted for in
the analysis. Moreover, the correlations or models generated from the data have not taken inlet
compressibility into account. This could possibly be one of the reasons for the large divergence of
the predictions made by the correlations and the poor performance when compared with large sets
of data from different experiments.
2.3 Thome’s three-zone model

Jacobi and Thome [26] presented an elongated bubble model in 2002, which was later refined in
2004 into a three-zone model by Thome et al. [27]. The conceptual picture behind the model is the
formation of Taylor-like flow. Due to heat transfer the film evaporates. If the thickness reaches a
minimum value (to be specified later) before the next bubble arrives dryout occurs locally and a
liquid slug, elongated bubble and vapour slug triplet is formed, as shown in Figure 2.
This three-zone model considers the heat transfer in each zone separately. It assumes that liquid
and vapour travel at the same velocity (homogenous flow) and a constant heat flux is applied at the
inner wall of the channel. Evaporation of the thin film between the gas bubble and the tube wall is
considered to be the dominant heat transfer mechanism. All energy that enters the fluid is used to
vaporize liquid. Based on this assumption, no fluid is ever superheated. The heat transfer
coefficients in the liquid slug and dry zone are calculated assuming single phase, steady, laminar
flow. The heat transfer coefficient across the thin film between the gas bubble and the tube wall is
calculated using a one-dimensional thermal conduction equation in the liquid film. The film
thickness at a fixed position is calculated from an initial film thickness and the evaporation rate of
the liquid, which is based on the heat conducted across the film.
Three important variables used in the model are not predicted theoretically. They are the
frequency of passage of the bubbles (f), the minimum liquid film thickness (dmin) at dryout and a
Lp
Lv Ll
R
Dry Elongated Liquid
zone bubble δ0 slug d
δmin
Ldry Lfilm
Figure 2. Thome’s three-zone heat transfer model for the elongated bubble flow regime
in microchannels. Taken from Thome et al. [27].
Volume 6 · Number 2 · 2014

correction factor (Cd0) used to predict the initial film thickness. The last two parameters are used
in the calculation of the local film thickness. As shown in Figure 2, the local film thickness starts
from the initial film thickness and jumps to zero when it reaches the minimum film thickness.
These three variables are defined by comparing the model outputs with the experimental database
which will be discussed below.
Dupont et al. [28] compared this three-zone model with nine sets of experimental data which
covered tube diameters ranging from 0.77–3.1 mm, mass velocities from 50 to 564 kg m-2 s-1, heat
fluxes from 5 to 178 kW m-2 and vapour qualities from 0.01 to 0.99. The authors showed that the
model predicted 70% of the data within a deviation of 30% by implementing an optimization
method for the parameters which specified an independent set of model parameter values for each
set of data. The optimization method for the free parameters used least-square error minimization
in which an objective function was created for each single dataset by curve fitting to the
experimental data. This method created different values for the same parameters for the various
datasets, and this brought two uncertainties.
First, there were five parameters being determined by the optimization process. Three of them
(aq, nq, nf) were created to express the relationship between the bubble frequency (f) and the wall
heat flux. The bubble frequency was defined as a power of the non-dimensional wall heat flux (the
wall heat flux divided by a reference heat flux), with nf being the exponent. nq and aq are parameters
used to define the reference heat flux as a function of saturation pressure, following the method
used in Cooper’s correlation [15]. The experimental analysis of Brutin et al. [29] showed
dependency of the bubble passage frequency on the wall heat flux and liquid mass flow rate. They
surmised the frequency to be a function of many parameters, such as wall heat flux, mass velocity,
fluid properties and the geometry of the test section. Hence, the prediction process for the bubble
frequency using the above three parameters which only included heat flux effects, is not physically
correct. The other two parameters, the correction factor (Cd0) on the initial film thickness and the
minimum film thickness (dmin) are associated with real physical values. These two thicknesses are
affected by many physical variables in the experiments, such as wall roughness, the surface tension
of the liquid and the evaporation rate. However, the optimization process gave values for the two
parameters by numerical fitting to the data without considering the above physics.
Second, the values of the parameters for each experimental dataset are very different for some
parameters, as shown by the range of values of the parameters shown in Table 1. This indicates that
the effect of each single parameter on the model predictions is not clear. The optimization method
applied to the model just performed mathematical fitting to each set of experimental data. Each
model with a specified set of values could only fit the dataset used to set the constants. Therefore,
the general model from Thome et al. [27] is not able to provide accurate prediction of flow boiling
in microchannels.
However, Thome’s model [27] is able to predict many datasets by changing the parameters for
each dataset [28], but these predictions cannot be correct without considering the inlet
compressibility which has been shown to have a large effect on microchannel boiling [25] and
definitely existed in many previous experiments.
Based on the above discussion it is clear that fitting empirical models is both difficult and leads
to the introduction of many uncertain parameters. This has led a number of groups to look at more
fundamental methods to study flow boiling in microchannels.
3. CFD FRAMEWORK OF FLOW BOILING

As experimental analysis at the micro-scale is limited by the difficulty of making measurements
and visualising rapid phase change phenomena, and theoretical predictions can only work for
simple mathematical models, numerical methods have been developed for engineering research.
Gupta et al. [35] provided a review of experimental and computational work on Taylor flow in
microchannels, where the CFD methods for two phase flow without phase change were
documented. Therefore, this section focuses on numerical approaches that can be used for flow
boiling in microchannels. In Sections 3.1 and 3.2, the governing equations and boundary conditions
for flow boiling in microchannels are presented. Section 3.3 describes numerical methods used for
microchannel flow boiling simulation. Finally, the phase change models used in the numerical
studies of microchannel flow boiling are discussed in Section 3.4.

Table 1. Values of parameters in optimization function for Thome’s three-zone model

with aq in W m-2 and dmin in microns. nq, nf and Cd0 are dimensionless. Adapted from
Dupont et al. [28].
Study Fluid Diameter Specific parameters for each study

(mm) aq nq nf Cd0 dmin
Lazarek and
Black [30] R113 3.1 3319 -0.54 1.8 1.176 0.01
Wambsganss
et al. [31] R113 2.92 3298 -0.47 1.72 0.77 0.47
Tran et al. [11] R12 2.46 3321 -0.57 1.71 0.87 0.36
Yan and
Lin [32] R-134a 2 3551 -0.44 1.27 1.2 1.92
Bao et al. [10] R11 1.95 3334 -0.59 1.72 0.59 0.02
Bao et al. [10] R123 1.95 3334 -0.55 1.73 0.44 0.08
Baird et al. [33] CO2 1.95 3323 -0.55 1.7 0.34 0.15
Lin et al. [20] R141b 1.1 3324 -0.6 1.79 0.86 0.58
Agostini [34] R-134a 0.77 3272 -0.05 1.18 1.23 1.82
Agostini [34] R-134a 2.01 3452 -0.24 1.7 0.8 3
General
Model 3328 -0.5 1.74 0.29 0.3
Parameter range (3319, (-0.6, (1.27, (0.29, (0.01,
3551) -0.24) 1.8) 1.23) 3)
The parameter ranges were added by the authors.
3.1 Governing Equations

The governing equations for two phase flow without phase change have been presented by Gupta
et al. [35]. For boiling flow, due to the interphase mass transfer, additional source terms are needed
in the conservation equations. The governing equations for flow boiling and heat transfer are shown
in eqns (1)–(3). There are some general assumptions made in the derivation of these equations:
1. All the fluids are Newtonian.
2. Viscous heating/dissipation is neglected.
3. The flow velocity is sufficiently low that a thermal energy, as opposed to a total energy,
equation can be used.
Continuity:
∂(α i ρi )
+ ∇ ⋅ (α i ρi vi ) = Sρi (1)
∂t
Momentum:
∂(α i ρi vi )
+ ∇ ⋅ (α i ρi vi ⊗ vi ) = −α i ∇Pi + ∇ ⋅ (α i τ i ) + α i ρi g + Fij + Fσ + Svi (2a)
∂t
where τ i = μi ( ∇vi + ∇vi )

T
(2b)
Energy:
∂(α i ρi ei )
+ ∇ ⋅ (α i ρi vi hi ) = ∇ ⋅ (α i ki ∇Ti ) + Shi (3)
∂t
where i = g, l denotes the phases, g is the gas phase and l is the liquid phase. ai is the volume
fraction, ri is the density, mi is the dynamic viscosity and ki is the thermal conductivity of the ith
phase. vi, Pi, Ti, ei and hi are the velocity, pressure, temperature, internal energy and enthalpy of the
ith phase. g is the acceleration due to gravity, t is time and the Fij denotes the momentum transfer

between the two phases due to forces such as drag, added mass, lift etc. The force Fs represents the
interfacial force due to surface tension. Sri is the phase change mass source, Svi is the momentum
transfer due to the mass source and Shi is the phase change energy source. The form taken by these
source terms is discussed later.
The implementation of the phase change mass source in the governing equations is dependent
on the interface capturing method used in the simulation. In the volume of fluid (VOF) method,
where a continuous interface is assumed to exist between the fluids, the mass source in eqn (1) can
be expressed as:
Sρ g = m ∇α i (4a)
Sρl = −m ∇α i
(4b)
.
where m is the phase change mass flux which is positive in evaporation/boiling processes and is
negative in condensation processes. |—ai| is the liquid-gas interfacial area density (the choice of
fluid used in the evaluation does not affect the value of |—ai|). The liquid-gas interfacial area
density is characterized by the interfacial area per unit volume between the two phases. The
integration of the interfacial area density over the computational domain gives the interfacial area
between the two phases, as shown in eqn (5):
∫ ∇α i dV = Aint (5)
v
where Aint is the interfacial area between the two phases [36]. A variety of functions can be used to
determine the interfacial area per unit volume depending on the degree of sharpness required. Any
function of the form f(a)|—ai| where f(a) satisfies the equation
1
(6)
∫ f (α ) dα = 1
0
provides an estimate of the interfacial area per unit volume [37]. Candidates for f(a) in increasing
order of sharpness are 1, 2a, 6a(1- a), etc. The inclusion of this source term over the computational
domain simulates the interphase mass transfer between the two phases.
The term Svi represents the momentum carried by a fluid as it changes phase. The term
represents the phase change momentum flux but in some literature is known as the evaporation
reaction force. The obvious treatment is to define the term as follows
Svg = Sρ g vg Sρ g < 0 (7a)
Svg = Sρ g vl Sρ g > 0 (7b)
Svl = −Svg
(7c)
so that each phase carries its own momentum with it during phase change.
For the energy conservation equations, the energy source terms for the gas and liquid phases are
calculated based on the mass source terms, the specific enthalpies of the fluids and the latent heat
of vaporisation as shown below.
In evaporation/boiling processes, Srg > 0. As liquid is removed from the liquid phase and the
same mass of vapour is added to the gas phase, the energy that needs to be removed from the liquid
phase is equal to the enthalpy of the liquid removed plus the latent heat of vaporisation at the liquid
temperature. To maintain energy conservation, the same amount of energy is added to the gas
phase. The energy source terms for the two phases in evaporation/boiling processes are expressed
as shown in eqn (8):

Shl = −Sρ g ( hl + hlv ) (8a)
Shg = Sρ g ( hl + hlv ) (8b)
where hlv is the latent heat of vaporisation.

In condensation processes, Srg < 0. As vapour is removed from the gas phase and the same mass
of liquid is added to the liquid phase, the energy that needs to be removed from the gas phase is
equal to the enthalpy of the vapour removed. To maintain energy conservation, the same amount of
energy is added to the liquid phase. The energy source terms for the two phases in condensation
processes are expressed as shown in eqn (9):
Shg = Sρ g hg (9a)
Shl = −Sρ g hg (9b)
3.1.1. Equations of state

Most simulations assume constant densities for the gas and liquid phases. This assumption limits
the CFD work to flows with relatively low velocities, as compressibility effects are not taken into
account. The speed of sound in two phase mixtures has been shown to be much lower than the
single phase sound speeds of the vapour and liquid separately [38] and are therefore potentially
important. Compressibility effects may become important at high mass fluxes, especially for flows
with relatively higher vapour quality (x > 0.1).
In order to take the effect of compressibility into account, equations of state for the gas and
liquid phases need to be specified. For water at low pressure and temperature, direct use of steam
table data is possible. However, for other materials or very high pressure and temperature
conditions where no tabular data are available, suitable equations of state need to be chosen.
3.2 Boundary conditions

Typically, at the inlet boundaries, the velocity, volume fraction and temperatures are specified. At
the outlet, an average relative pressure of zero is widely used [39-41]. An outflow boundary
condition [42] which gives a zero gradient boundary condition for the velocity and temperature is
used in some simulations [43, 44].
3.2.1. Wall conditions

Due to the presence of fluid viscosity, the no-slip condition is applied at the wall:
n × (v − vw ) = 0 (10)
where n is the unit normal to the wall and nw is the velocity of the wall.
The wall is impermeable so there is no flow into the wall, giving:
n ⋅ v = n ⋅ vw
(11)
The thermal boundary condition at the wall can be adiabatic (zero heat flux), a specified heat
flux, a specified temperature or a heat transfer coefficient with a specified ambient temperature.
3.2.2. Interface conditions

For two fluids separated by an interface, the conditions on the normal and tangential velocity fields
at the interface are as follows:
ρ g n ⋅ v g = ρ l n ⋅ vl (12)

n × ( vg − vl ) = 0 (13)
where vg and vl are the velocities of gas and liquid at the interface and n is vector normal to the
interface. Note that whilst the tangential velocities remain equal there must be a jump in the normal
velocity to account for any density change upon phase change.
The normal and tangential forces must balance at the interface giving eqns (14) and (15) [45].
The variation of surface tension coefficient with temperature needs to be modelled if the
temperature changes significantly along the interface [46].
−Pl + Pg + n ⋅ (τ l − τ g ) = σκ (14)
n × (τ l − τ g ) = ∇σ (15)
here s is the surface tension coefficient, k is the curvature of the interface and ti are the interfacial
shear stresses.
It is important to note that these interfacial conditions are not implemented directly in an
interface capturing model, as no explicit interface is present. They have to be built into the
numerical method.
3.3 Numerical Methods for Flow Boiling Simulation

Many numerical methods have been developed for multiphase flow modelling, such as boundary
integral methods, finite element methods, interface capturing methods, etc. Detailed descriptions of
these methods can be found in Gupta et al. [35]. This article is focused on flow boiling simulations
which require integration of phase change with boundary determination. Therefore, this section
focusses on the interface capturing methods that have been used for this purpose. However, an
interface tracking method is also mentioned because of its potential to improve the interface
modelling accuracy.
In interface capturing methods, the homogeneous assumption can be made for the velocity field,
temperature field or both of them. If the homogeneous model is applied to the velocity and
temperature, all the fluids are assumed to have the same velocity and temperature at any given
point. The bulk properties of the fluids are calculated as volume fraction weighted averages of the
properties of the two fluids. Based on these assumptions, the governing equations of the two fluids
can be added to form a single-set of equations, as given in eqns (16)-(18), using the bulk properties
of the fluids.
∂ρ
+ ∇ ⋅ ( ρ v ) = 0 (16)
∂t
∂( ρ v )
+ ∇ ⋅ ( ρ v ⊗ v ) = −∇P + ∇ ⋅ μ ( ∇v + ∇v T ) + Fσ + ρ g
( ) (17)
∂t
∂E  = ∇ ⋅ k ∇T
∂t
( )
+ ∇ ⋅ vH ( ) (18)
The homogenous equations are readily derived from the transport equations for the separate
phases by summing them. As noted above properties are then calculated based on the phase volume
~
fraction weighted values, e.g. r = al rl + ag rg. The homogeneous energy equation is somewhat
~
more complicated as the internal energy and enthalpy are defined via E = al rl el + ag rg hg and
~
H = al rl hl + ag rg hg and the temperature is obtained from the caloric equation under the constraint
that the temperature is the same for both phases.
Additionally, a transport equation for a colour function is solved to identify the interface location
and to calculate the volume fractions, as given below:
∂ ρi C
+ v ⋅ ∇ρi C = Sρi (19)
∂t

where C is the colour function. The definition of the colour function is different in different
interface capturing methods. The volume of fluid method and the level-set method are briefly
reviewed in the following sections.
As the mass transfer rate is high in flow boiling conditions, numerical instabilities may occur at
the interface when the above equations are solved numerically. This is particularly so if a high
evaporation flux is calculated in a cell that contains relatively little liquid. A numerical smoothing
procedure proposed by Hardt and Wondra [36] can be used to enhance the numerical stability by
means of shifting such a source to the liquid side of the interface and smearing the mass source term
(Sri) over a few cells near the interface.
The surface tension force (Fs), which is the embodiment of the interface condition (eqn(14))
within the numerical method, is normally modelled by the Continuum Surface Force (CSF) model
proposed by Brackbill et al. [47]. In this model, the surface tension is calculated as follows:
Fσ = σκδ ( r − rint ) n (20)

∇C
n= ; κ = ∇⋅n
∇C (21)
here d(r - rint) is the Dirac delta function. Lafaurie et al. [48] and Harvie et al. [49] pointed out that
this implementation of the surface tension force can induce unphysical velocities near the interface,
known as ‘spurious or parasitic currents’. In the review of CFD modelling of microchannel flows
performed by Fletcher et al. [50], they pointed out that the main challenges in the modelling of
multiphase flow in microchannels are the correct identification of the interface location and the
minimisation of the parasitic currents near the interface that arise due to surface tension modelling.
Magnini et al. [43] proposed a height function algorithm to replace the curvature calculation
method in eqn (21) to improve the performance of surface tension modelling, with reference to
Malik et al. [51] and Hernandez et al. [52]. Magnini et al. [43, 44] successfully performed flow
boiling simulations by implementing this height function algorithm in ANSYS Fluent.
In addition, both the interfacial momentum transfer terms and the source terms due to phase
change cancel out in the summed momentum equations yielding great simplification. Similarly, the
addition of the two energy equations for the two fluids (eqn (3)) has cancelled out the energy source
terms (eqns (8) or (9)) and therefore yields eqn (18). This is correct as the summed energy source
needs to be zero in order to maintain the energy conservation of the system. It is important that
when these energy sources are implemented in any CFD code that care is taken to ensure that the
different fluids have the same reference conditions for the enthalpy and that the exact manner in
which source terms are implemented is understood. It is very easy to introduce artificial source
terms if great care is not taken with respect to these points. It should be noted that if the
homogeneous (summed) energy equation is used the liquid and vapour are assumed to be in local
thermal equilibrium at every point in the flow.
3.3.1. Volume of fluid (VOF) method

The volume of fluid (VOF) method proposed by Hirt and Nichols [53] is a widely used technique
for interface capturing in multiphase flow simulations. In this method, the presence of the two
phases (which may in fact be liquid, gaseous or a multi-component mixture of either gas or liquid
but not a mixture of both) is tracked using a volume fraction variable (a) to be the colour function
in eqn (19). For the ith phase in a specific cell, ai = 1 denotes that the cell is full of the ith phase,
whereas ai = 0 indicates no ith phase in that cell. Therefore, the interface can be captured at the cells
with value 0 < ai < 1, using numerical methods that use geometrical reconstruction schemes
[54-56] or high order difference schemes [57] to keep track of the interface. As the sum of the
volume fractions of the two fluids must be unity everywhere, the volume fraction for one fluid is
solved using eqn (19) and the other is obtained from ai = 1 - aj.
The VOF method is easy to implement into the simulations of complex problems and has been
successfully used in many simulations of flow boiling in microchannels [43, 44, 58-60]. However,
one thing that needs to be noted is that the performance of different numerical schemes used in
interface capturing can be different. Typically compressive differencing schemes, in which
numerical methods that include anti-diffusion are used for convection of the colour function,
provide an interface with diffusion of over 3-4 cells, whereas the interface is located inside one cell

using geometric reconstruction schemes. This fact causes a high dependency of the applied surface
tension force on the grid size and therefore results in differences in performance of simulations
using different interface tracking schemes with surface tension force being considered [35].
3.3.2. Level-set method

This method was first proposed by Osher and Sethian [61] and then improved by Sussman
et al. [62]. Different from the VOF method, the level-set method uses a level-set function (f) as the
colour function in eqn (19) to define the interface. The value of the level-set function (f) is positive
in one phase and negative in the other phase, where a zero value of f indicates an interface. The
absolute value of the level-set function (f) at a given location represents the shortest distance from
that location to the interface. As volume fractions are required in the solution of the governing
equations, a smoothed Heaviside function is used to compute the volume fractions of the fluids at
the interface. Simulations of flow boiling have been successfully performed using the level-set
method [39-41, 63-66].
The level-set method has been shown to be excellent for the calculation of interface curvature
[35]. However, Rider and Kothe [54] pointed out a disadvantage of the level-set method that loss
of mass conservation could happen in simulations using this method. Recently, Nourgaliev et al.
[67] suggested that by using a high order discretization scheme to solve the advection equations
and applying finer mesh near the interface, the mass conservation in level-set methods can be
improved.
3.3.3. Marker point method

Unverdi and Tryggvason [68] developed a marker point method, which was later improved by
Tryggvason et al. [69]. This technique is also known as a front tracking method because it directly
tracks the interface with a Lagrangian grid. In this method, the centre of the smeared interface is
marked by many points and the complete interface is defined by connecting these marker points by
lines (2-D) or triangular surfaces (3-D). The marker points are advected so that the interface can be
explicitly tracked. The material properties at the interface are smoothed and need to be reset at
every time step.
In this method the surface tension force is obtained from the tensile forces. The surface tension
force on a two-dimensional interface can be expressed as given in eqn (22):
∂t
Fσ = ∫ σκ n ds = σ ∫ ∂S ds = σ (t 2 − t1 ) (22)
ΔS ΔS
where ti denotes the tangents of the marker points. Therefore, instead of calculating the curvature,
only the calculation of the tangents of the interface marker points is needed and this is
straightforward in the front tracking method because the interface is tracked explicitly with the
Lagrangian marker points. The main disadvantage of this method is that when the interface is
stretched or compressed the marker points must be redistributed. In addition, this method is
computationally more expensive than the interface capturing methods, such as the VOF and level-
set methods.
Tryggvason et al. [69] have used this method to perform simple film boiling simulations.
Baltussen et al. [70] extracted the tensile force model of this method and applied it to the VOF
model to compare with other surface tension models commonly used in the VOF method. The
results showed that the tensile force model can be applied successfully using the reconstructed
interface in the VOF model, and the accuracy of surface tension calculation was improved
significantly compared with the CSF model. However, there has not been any numerical work on
microchannel flow boiling that uses this approach and no commercial code has this approach
embedded in it to date.
3.3.4. Discussion
There are other interface capturing techniques, such as phase field methods [71], lattice Boltzmann
methods (LBM), etc. Some of these methods are still in development. More discussion of these
methods can be found in the review by Gupta et al. [35]. In summary, the VOF and level-set
methods are the two most frequently used interface capturing methods in complex multiphase flow
modelling, including the modelling of flow boiling in microchannels, because of their high level of

development and the balance between computational cost and calculation accuracy. Currently, the
VOF method is the most widely implemented interface capturing method in commercial CFD
codes. The VOF method is computationally cheap relative to the level-set method, it conserves
mass but has no explicit distance function so it can be hard to determine the interface curvature
accurately. The level-set method gives a distance function which can be used for the interface
curvature calculation but suffers mass conservation issues. It is computationally more expensive
than the VOF method because of a redistancing step required after each time step and may in fact
not necessarily give a better calculation of curvature. A possible strategy is to couple the level-set
and VOF methods together in the calculation to combine their advantageous parts [72]. In addition,
the front tracking method also has potential to be used in numerical simulations of microchannel
flow boiling.
3.4 Phase Change Models

The phase change models used in the simulations are critical because they decide the way in which
the interfacial mass transfer is handled. Various approaches have been used in dealing with phase
change in the studied literature and can be generally split into three classes.
3.4.1. Phase change models based on experimental data and theoretical models
The first class of phase change models is based on experimental data or theoretical models. In these
methods, a bubble growth rate is calculated from either a theoretical model [73] or an experimental
dataset [74]. Then this bubble growth rate is imposed upon the simulation by means of applying a
heat flux at the liquid-gas interface from which a mass flux is calculated [58, 60] or injecting a
vapour mass flow into the tube [59]. The applied heat flux and the injected vapour mass flow rate
are back-calculated from the bubble growth rate. Typically, once the bubble diameter reaches the
tube diameter, the imposition of a bubble growth rate is stopped and the phase change is driven by
the heat conduction across the liquid film between the interface and the tube wall.
Zhuan and Wang [58, 60] and Zu et al. [59] used this type of phase change model. The
imposition of the bubble growth rate provides simulations with bubble growth processes similar to
the experimental observations or theoretical calculations. Then the simulations can concentrate on
the modelling of the processes that follow, such as film evaporation. However, the requirement of
the experimental data and theoretical models limits the utility of this type of model.
3.4.2. Phase change models based on heat conduction

The second class of phase change models is based on the assumption that the phase change at the
liquid-gas interface is driven by the difference between the energy conducted to and from the
liquid-gas interface. As shown in Figure 3, the heat energy available for phase change is
Wall
ql = n.kl ∇Tl
Liquid
Interface
Gas
qg = n. kg ∇Tg m = (ql −qg )/hlv
Figure 3. Phase change is driven by the heat conduction across the interface.

.
proportional to the heat conducted from the liquid phase to the interface (q l) minus the heat
.
conducted from the interface to the gas phase (q g).
.
Therefore, the mass flux due to phase change (m) at the interface can be calculated as shown in
eqn (23):
n ⋅ (kl ∇Tl - kg ∇Tg )

m = (23)
hlv
.
Again, m is positive in evaporation/boiling process and is negative in condensation process.
In some numerical work [75], the heat conduction from the interface to the gas phase was
assumed to be negligible. Therefore, if the heat conduction from the interface to the gas phase is
set to zero, the phase change mass flux can be calculated as shown in eqn (24). The work performed
by Mukherjee and co-workers [63, 64, 66, 76] and Son’s group [39-41] used this approach to model
the phase change.
n ⋅ kl ∇Tl (24)
m =
hlv
Note that in this formulation the mass source term is calculated on the basis of local energy
fluxes and contains no information about the behavior of the interface. No account is taken of any
possible rate limiting processes arising from kinetic effects at the interface nor is any temperature
specified at the gas-liquid interface.
3.4.3. Kinetic-based phase change models

Kinetic-based models form the third class of phase change models. These models are derived from
the kinetic theories of liquid and gas. Schrage [77] firstly analysed the mass transfer process at the
liquid-gas interface using the Maxwellian distribution-based kinetic theory of gases. Tanasawa [78]
reviewed the work of Schrage [77] and other workers [79], and provided an interphase mass
transfer model for interphase condensation as shown below.
m = ϕ (Tint − Tv ) (25)
j, defined in eqn (26), is the kinetic mobility which was derived from the kinetic theory of gases,
assuming the temperature jump across the liquid-vapour interface is very small relative to the absolute
vapour temperature. The Clausius–Clapeyron relation was utilized in the derivation to change the
pressure difference term into a temperature difference term. The kinetic mobility is given by
1/2
2γ ⎛ 1 ⎞ ρv hlv
ϕ= ⎜ ⎟ (26)
2 − γ ⎝ 2 π R ⎠ Tv3/2
where R is the gas constant, Tint is the interface temperature, rv is the vapour density and g is the
phase change or accommodation coefficient. Tv is the vapour temperature and is normally equal to
the gas phase temperature (Tg) in a multicomponent gas. Tanasawa [78] pointed out that the
measurement of g is difficult and the values of g range from 0.04 to 1 depending on the fluid.
The Tanasawa [78] model has been used in some flow boiling and evaporation simulations
[43, 44, 65]. However, this model is only applicable to the simulations in which the phase change
happens at conditions close to saturation. This is because the Clausius–Clapeyron relation, which
is applicable to saturation conditions, was used in the derivation of this model. Actually, a more
general model without using the Clausius–Clapeyron relation can be obtained from Tanasawa [78],
as shown in eqn (27):
2γ Pv − Psat
m = (27)
2 − γ 2 π RTv
where Pv and Psat are the vapour pressure and the saturation pressure at the interface, respectively.
Note, if the gas phase is a multicomponent mixture, Pv denotes the partial pressure of the vapour.
This model can be applied to evaporation as well as boiling and can account for the presence of
permanent gases.

The assumption that the vapour-liquid interface is at the saturation temperature is valid over a
wide range of conditions. Tanasawa [78] shows that for condensation of water at 1 atm the
interfacial temperature drop is 0.007 K which is clearly insignificant, although for pressures of
0.001 atm and liquid metals it can be ~10 K. Hardt and Wondra [36] have combined eqns (24) and
(25) to show that for very thin liquid films the interface temperature can be elevated, corresponding
to a situation in which the very high heat flux across the liquid film cannot be carried away by the
phase change process sufficiently quickly for equilibrium to be maintained.
When implementing eqn (25) into a CFD model of boiling the interface temperature must either
be specified as Tsat or it can be set equal to Tl so that some superheating of the liquid must be
present for phase change to occur.
3.4.4. Discussion
The above three categories of phase change models include most of those used for microchannel
flow boiling simulations. The models based on experimental data and theoretical considerations are
not easy to implement because they require experimental datasets and analytical correlations for
each specific case being simulated. For any experimental dataset, the construction of a theoretical
model from the experimental data is needed before it can be implemented in the simulation. The
heat conduction models are relatively simple and easy to implement but they are potentially
inaccurate because the phase change rate only depends on the heat conduction in the phases,
therefore no interface kinetic limitations are captured. The kinetic-based models are the most
accurate amongst the three types of models as they calculate the phase change flux based on the
physics at the interface and therefore capture the kinetic limitations. However, when the kinetic-
based models are chosen, great care needs to be taken to ensure the appropriateness of the model
for the particular simulation case. For example, if eqn (26) is to be used in a simulation, then the
interface condition needs to be always very close to saturation, as the Clausius–Clapeyron relation
has been included in the model.
4. NUMERICAL WORK FOR MICROCHANNEL FLOW BOILING

Although local data from experiments are relatively sparse, plenty of global data obtained from
experiments and analytical solutions offer an abundant database to be compared with simulation
results. Many numerical simulations using various numerical methods have been performed to
study important aspects of boiling flow: bubble growth [63, 76], bubble formation and the
subsequent bubble dynamics [41, 58-60, 64, 65], bubble merging [66], evaporation of elongated
bubbles [43, 44, 59], flow boiling in finned channels [39] and parallel channels [40]. Some of these
studies are reviewed in the next sections.
4.1 Bubble Dynamics Simulations by Kandlikar and co-workers

Mukherjee and Kandlikar [63] performed 3-D transient simulations of vapour bubble growth
during flow boiling of water in a microchannel. The software used for the simulations is not
specified by the authors but the implemented numerical methods are explained in detail. The
Navier-Stokes equations were solved using the SIMPLER algorithm of Patankar [42]. The gas-
liquid interface was captured using the level-set method developed by Sussman et al. [62]. A fifth-
order scheme proposed by Fedkiw et al. [80] was used for the discretization of the level-set
function f. The transient timestep was determined from eqn (28), where Dx is the grid spacing and
u0 is the inlet liquid velocity, essentially enforcing a pseudo-Courant number of 0.5.
Δx
Δt = (28)
2u 0
The vapour temperature in the simulation was assumed to be constant at the saturation
temperature. The phase change model presented in Section 3.4.2 was used in this simulation,
assuming the heat conduction from the interface to the gas phase is negligible. The homogeneous
temperature assumption was used in this simulation. Therefore the phase change mass flux at the
interface was calculated as shown in eqn (24) using the bulk temperature.
The computational domain used is shown in Figure 4. Flow is in the x direction. A spherical
vapour bubble was placed at x = 0.99, y = 0, and z = 0, with radius 0.1 in the reference frame shown

X
Z
0.5
0
5
4
−0.5 3
−0.5 2
0 1
0.5 0 0.000 ms
Figure 4. Computational domain and initial bubble position. Taken from Mukherjee
and Kandlikar [63].
in Figure 4. Note, the coordinates in the figure are made non-dimensional using a length factor of
200 mm. The bubble is equidistant from the walls in the y and the z directions. 480, 96 and 48 cells
are used in the x, y and z directions, respectively. The simulations were terminated when the vapour
bubble nose reached an axial position of four-fifths of the channel length.
The channel walls were assumed to be no-slip. All the velocities in the domain were set to zero
at t = 0 s. The temperatures of the walls were set to 107 °C. The vapour temperature in the bubble
was assumed to remain at the saturation temperature of 100 °C. The liquid temperature in the
domain was set equal to the inlet liquid temperature. The local axial velocity at the inlet, u0, was
calculated from a specified Reynolds number. Inlet liquid temperatures of 102 °C, 104 °C and 107
°C were studied in different runs, which represent different levels of liquid superheating. The
Reynolds numbers used were 50, 100 and 150. If the bubble touched the wall and formed dry-
patches at the wall, the contact angle between the liquid-vapour interface and the wall was set to
50°. The gravitational force was neglected in all but the last simulation in which the authors tested
the effects of gravitational force on the bubble growth rate. Both fluids were assumed to be
incompressible.
In the first run, the inlet liquid temperature was set to 102 °C and the Reynolds number of the
inlet liquid was 100. At the beginning of the simulation, the vapour bubble was observed to grow
spherically with a linear time-law for the radius. As the bubble grew bigger, its shape changed from
spherical to being elongated in the axial direction. Subsequently, the bubble stretched and generated
a thin liquid film between the gas-liquid interface and the channel wall, as shown in Figure 5. This
observation of bubble shape evolution agrees qualitatively with the experimental data obtained by
Kandlikar and Balasubramanian [81] and Xu et al. [82]. A quantitative comparison was not possible
because the physical parameters used in the experiments were different from those used in the
simulations.
Figure 6 shows the change of bubble growth rate by plotting the bubble diameter and the axial
bubble length against time. The equivalent diameter was calculated from the volume of the bubble,
assuming the bubble remained spherical. The bubble length shown in the figure has a dramatic
increase after around t = 0.2 ms. Figure 6 shows that after around t = 0.232 ms, the bubble is
elongated and appears to grow in the axial direction. A thin liquid layer is formed between the walls
and the liquid-vapour interface. Mukherjee and Kandlikar [63] believe that the dramatic increase
of bubble length after around t = 0.2 ms is due to this thin liquid layer between the walls and the
liquid-vapour interface, where a high rate of evaporation takes place. In fact, the equivalent
diameter of the bubble shown in Figure 6 shows an approximately linear growth law during the
entire simulation period from 0 ms to around 0.28 ms, with only a slight increase after t = 0.2 ms.
This indicates that the bubble growth rate does not increase very much when the bubble starts to
become elongated compared with the growth rate when the bubble shape is spherical. Thus the
evaporation rate of the film cannot be shown to be very high using these data.
A parametric study of the effects of liquid superheat and inlet Reynolds number on the bubble
growth rate was performed. The results showed that higher superheats significantly increased the

X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.000 ms
X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.082 ms
X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.158 ms
X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.232 ms
X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.256 ms
X
Z
0.5
0
5
−0.5 4
3
−0.5 2
0 1
0.5 0 0.278 ms
Figure 5. Growth of a vapour bubble in a microchannel for Tin = 102 °C and Re = 100.
Taken from Mukherjee and Kandlikar [63].
bubble growth rate, whereas the inlet liquid Reynolds number did not affect the growth rate much.
The effect of gravity on the bubble growth rate was also evaluated. Under the specific conditions
of an inlet liquid temperature of 107 °C and an inlet Reynolds number of 100, the effect of gravity
was found to be negligible on the bubble growth rate, as would be expected.

Equivalent diameter
1 Bubble length
Upstream end location
0.9 Downstream end location
0.8
Dimensions (mm) 0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4
Time (ms)
Figure 6. Bubble growth with time. Taken from Mukherjee and Kandlikar [63].
(a) 0.4 (b) 0.4

∆T − 5K Re − 50
Equivalent diameter (mm)
∆T − 8K Re − 100
0.3 0.3
∆T − 10K Re − 200
0.2 0.2
0.1 0.1
0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Time (ms) Time (ms)
(c) (d) 0.4

0.4
σ − 0.03 N/M ϕ − 20
σ − 0.04 N/M ϕ − 40
0.3 σ − 0.05 N/M 0.3 ϕ − 60
ϕ − 80
σ − 0.0589 N/M
0.2 0.2
0.1 0.1
0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Time (ms) Time (ms)
Figure 7. Bubble growth rate with varied (a) wall superheats; (b) liquid Reynolds
number; (c) surface tension; (d) contact angle. Taken from Mukherjee et al. [64].
Mukherjee et al. [64] conducted a detailed parametric study of the flow boiling of a water
vapour bubble in a square microchannel in 2011. This work used the same numerical model as
developed by Mukherjee and Kandlikar [63], described above. There is one difference in the new
work in that the vapour bubble was placed at the centre of a wall and touched the wall at a point,
representing nucleation at the wall.

Figure 7 shows the bubble equivalent diameter as a function of time for different values of four
parameters: wall superheat, liquid Reynolds number, surface tension and contact angle between the
liquid-vapor interface and the wall. The bubble growth rate is shown to increase with the wall
superheat. This is in accord with the author’s previous simulation results [63]. An increase in liquid
Reynolds number is shown to have a slightly negative effect on the bubble growth rate. The effect
of surface tension on the bubble growth is found to be negligible in the simulations. The contact
angle between the liquid-vapour interface and the superheated wall shows a remarkable effect on
evaporative heat transfer because the bubble growth rate yields much higher values for the
condition of a contact angle of 20° compared with those obtained from the simulations with contact
angles of 40°, 60° and 80°. This is because in the simulation period dryout on all the walls can be
seen for all the cases except for the 20° case in which only a small vapour patch occurs on the wall
touched by the bubble initially. Thus a higher contact angle is supposed to result in a larger wall
area exposed to the vapour and therefore leads to a lower wall heat transfer rate, but the dryout in
this simulation may be caused by numerical effects, for example the failure in capturing the thin
liquid films with the adopted mesh resolution.
The authors also evaluated the effect of the above parameters on the wall Nusselt number,
calculated based on the area-averaged heat transfer coefficient. The local heat transfer coefficient
is calculated from one-dimensional heat conduction across the liquid film between the wall and the
liquid-vapour interface, assuming that the wall temperature and the gas-vapour interfacial
temperature are constant. Figure 8 shows the variation of Nusselt number on the north side wall
(the wall opposite to the wall touched by the bubble initially) with time obtained from simulations
performed for different values of the four parameters. The variations of the Nusselt number are
similar to those observed for the bubble growth rate: high liquid superheats and small contact
angles give high Nusselt numbers, whereas the effects of liquid Reynolds number and surface
tension on the heat transfer are negligible. The effect of Reynolds number on Nusselt number is
qualitatively in agreement with some experimental observations that the heat transfer coefficient is
weakly dependent on mass flow rate [11-13].
The above work performed by Kandlikar and co-workers showed that bubble growth could be
(a) 50 (b) 50
∆T − 5K Re − 50
40 ∆T − 8K 40 Re − 100
∆T − 10K Re − 200
Nu north
Nu north
30 30
20 20
10 10
0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Time (ms) Time (ms)
(c) 50 (d) 50
σ − 0.03 N/M ϕ − 20
σ − 0.04 N/M ϕ − 40
40 40 ϕ − 60
σ − 0.05 N/M
σ − 0.0589 N/M ϕ − 80
Nu north
Nu north
30 30
20 20
10 10
0 0
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
Time (ms) Time (ms)
Figure 8. Nusselt numbers of the north side wall for various (a) wall superheats; (b) liquid Reynolds
numbers; (c) surface tension coefficients; (d) contact angles. Taken from Mukherjee et al. [64].

simulated successfully during flow boiling in a rectangular microchannel. However, two

assumptions used in their work limit the validity of the CFD model. First, the velocities inside the
domain were set to zero initially and the liquid local velocity was constant across the inlet. This
assumption makes the liquid flow transient when the simulation begins. Thus the bubble formation
is not happening in a steady laminar flow. Hence, the results obtained from the simulations with
this assumption are not able to be compared with the experimental data obtained for bubble growth
in steady laminar flow. Second, the evaporation rate in eqn (24) is determined from the heat
conductivity of the liquid and the interfacial temperature gradient. This evaporation model ignores
any limitations from kinetic effects and therefore is not able to capture any kinetic limitations, such
as variation in evaporation rate caused by local saturation pressure/temperature changes, although
these are likely to be small.
4.2 Bubble Growth and Flow Pattern Analysis from Shanghai Jiaotong University
Group
Zhuan and Wang [58] studied bubble formation, coalescence and flow pattern transitions within
circular microchannels in flow boiling conditions by performing 3-D transient simulations. The
gas-liquid interface was captured by use of the Volume of Fluid algorithm (VOF). The
computational domain was a horizontal circular tube 0.5 mm in diameter and 70.7 mm in length.
Based on a grid independence analysis, the grid spacing was determined as h = 20 mm for all three
directions. The timesteps chosen for calculations were from 10-6 s to 10-7 s. The authors proposed
a theoretical model which predicted the bubble lift-off size to be around 100 mm so that they
believed this grid size would successfully capture the bubble interface. However, the simulation
pictures showed inaccuracy in interface capturing which will be discussed later.
The evaporation was handled in two different ways for nucleation bubble growth and slug
bubble growth, respectively. First, the occurrence of bubble nucleation was controlled by a wall
superheat equation [83]. Then a theoretical model for bubble growth rate was used to simulate the
bubble growth [73]. The evaporation heat transfer rate at the bubble surface was back-calculated
from the bubble growth rate, assuming the bubble shape was spherical, as shown in eqn (29), where
.
m denotes the rate of bubble mass increase. The rate of mass increase dmb/dt can be derived from
the bubble radius growth rate given by the above model [73].
dmb
qb = hlv m = hlv (29)
dt
Second, once the bubble diameter reached the tube diameter, the heat flux at the bubble interface
was calculated assuming steady-state heat conduction from the wall to the slug surface via a liquid film.
The fluid used in the simulation was R-134a, with a saturation temperature set at 30 ºC and an
inlet liquid temperature of 27 ºC. The simulations covered inlet liquid mass velocities from 350 to
2000 kg m-2 s-1, wall heat fluxes of 4 to 129 kW m-2 and vapor qualities ranges from 0.004 to 0.85.
Figure 9 depicts a horizontal section of bubbly flow from their simulation. Note the bubbles
shown in the picture have diffused boundaries which correspond to gas-liquid mixture regions.
These regions suggest that the gas-liquid interfaces were not successfully captured because the grid
spacing used in this simulation was not fine enough.
The monitored flow transition processes in flow boiling was qualitatively in agreement with
previous experimental observations [17, 84]. The simulation results also showed that for a high
wall heat flux, the bubble lift-off size increased and transition of the flow regime from bubbly flow
to slug-annular and annular flow was faster than the transition happening at a low wall heat flux.
The effects of mass flow rate on bubble growth and flow pattern transitions were not clearly shown
in this work. One thing should be noted is that during the bubble growth process, the bubble
Figure 9. A horizontal section of bubbly flow for G = 500 kg m-2 s-1. The gas volume
fractions are equal to 0 in dark blue regions and 1 in red regions. Taken from Zhuan and
Wang [58].

expansion rate was controlled by a theoretical model in this simulation. This choice limits the
applicability of the model because the theoretical model used to control the bubble growth was
constructed from the experimental data. Also many of the simulations appear to be under resolved
with very diffuse gas-liquid interfaces and direct contact of vapour with the walls.
4.3 Simulations by the Brunel University Group

A 3-D, transient numerical model of nucleation and growth of a bubble in a confined microchannel
was proposed by Zu et al. [59]. The simulations were performed using ANSYS Fluent. The VOF
method was used to capture the gas-liquid interface.
The boiling model formulated in these simulations involved some approximations based on a
concept of “pseudo-nucleate boiling” which separated the bubble evolution into two stages. The
first stage was the formation and growth of the bubble on the wall before lift-off. During this
period, the volumetric growth rate of the bubble was controlled to match the experimental
observations or results from a theoretical model. This was achieved by injecting vapour through a
small inlet located on the tube wall, as shown in Figure 10. The injection rate was calculated from
the author’s previous experimental observations [74]. When the lift-off time given by the
experiments was reached in the simulation, the injection stopped and stage one finished. During the
second stage (after the bubble detached from the wall), the simulated bubble growth was changed
to evaporation driven by a constant heat flux equal to the average wall heat flux, imposed on the
contact area between the bubble and the walls, assuming a thin film existed there to prevent dryout
from happening. However, the authors did not explain how they implemented the heat source at the
interface equal to the average wall heat flux with the use of the VOF method in such simulations.
The computational domain is shown in Figure 16. The channel dimensions are 40 mm, 1.6 mm
and 0.38 mm in the x, y and z directions, respectively. The number of cells used in each of the three
directions was 225, 50 and 20, giving cell sizes of 178 mm, 32 mm and 19 mm, respectively. The
transient timestep used in the simulation was 10-6 s.
For boundary conditions, a constant mass flux of 747 kg m-2 s-1 of hot water at T = 370 K entered at
the inlet. This mass flux gives an average inlet velocity of 0.747 m s-1 which corresponds to a Reynolds
number of 428. Therefore the flow is laminar. However, the standard k-e turbulence model was used in
the simulation by the authors because they expected the flow around the bubble would become turbulent
during the bubble growth. The accuracy of this assumption is discussed later. A constant heat flux of
210 kW m-2 was applied at the wall. Both fluids were assumed incompressible.
A preliminary steady-state liquid single phase simulation was performed to prepare the initial
condition for the transient simulations. This simulation avoided the inaccuracy of modelling bubble
nucleation by injecting vapour at the wall. The authors observed a liquid film between the bubble
and the tube after the bubble radius approached the channel wall. This is in accord with some
experimental observations [74, 82]. The pressure distribution along the domain obtained in the
simulation was found to match the 1-D model proposed in part I of the authors’ paper [74].
However, the bubble growth stage in this simulation is a reconstruction of the author’s
experimental observations so that this reproduction is not surprising.
Additionally, two assumptions made by the authors limit the validity of the CFD model. First,
the model assumes that the heat flux at the contact area between the bubble and the wall is equal
Outlet
Vapour inlet
z
x
g
y
Liquid inlet
Figure 10. Vapour inlet shown in the computational domain. Taken from Zu et al. [59].

to the wall heat flux and all the heat is used to drive the evaporation process. This assumption may
not be true and the method of implementation of the vaporization source term in the simulation is
not provided. Second, the standard k-e turbulence model was used to account for turbulence caused
by the high liquid and vapour velocities following bubble growth that the authors expected.
However, Figure 11 shows that the liquid velocities are less than 1 m s-1 over the flow field which
corresponds to a Reynolds number of 573 under the conditions specified in the simulation, making
the assumption of turbulent flow incorrect.
4.4 Simulations of Bubble Growth in Flow Channels by Son and co-workers
Lee and Son [41] performed 3-D transient analyses of the growth of a water vapour bubble in
rectangular microchannels using a numerical approach developed by Son et al. [85]. This approach
used the level-set method to capture the liquid-vapour interface.
The phase change model based on heat conduction (Section 3.4.2) was used in the simulations,
assuming that the vapour temperature was constant at the saturation temperature. The homogenous
temperature assumption was applied. Thus the mass flux due to phase change was calculated as
expressed in eqn (24) using the bulk temperature. In addition, an additional mass source term,
which was not included by eqn (24), is introduced into the mass conservation equation, to compute
1.00
0.80
0.60
0.40
0.20
0.00
t = 3.9 ms
Figure 11. Pathlines colored with velocity magnitude (m s-1). Taken from Zu et al. [59].
z = –W/2
y=H
g y
z x
x=L
z = W/2 Bubble
z = −W/2
Top
view z x
z = W/2
y=H y
Front Bubble Side
y view
view ϕ
x x=L z
Microlayer
h/2
δ
Figure 12. Computational domain and a detailed picture of the microlayer between the
bubble and the bottom wall of the channel. The distance h denotes the grid spacing. Taken
from Lee and Son [41].

the evaporation mass transfer in a microlayer between the bubble and the bottom wall near the
bubble-wall contact location, as shown in Figure 12. Lee and Son [41] considered that the
evaporation in the microlayer which had a thickness d less than h/2 was not computed in the heat
conduction model because the microlayer thickness was smaller than half of the grid spacing. They
used a theoretical model from Son et al. [85] to calculate this additional term. However, in the later
work performed by Lee et al. [39], this additional term was abandoned and the phase change was
handled by the heat conduction model (Section 3.4.2) directly. In fact, the accuracy of the
simulation is highly dependent on the grid spacing. If the grid spacing is made small enough to
successfully capture the liquid film between the bubble and the wall, the heat conduction model
should be able to handle the phase change in the microlayer presented in this model. The grid
spacing of this simulation will be discussed later when looking at the computational results.
Figure 12 shows the computational domain of the simulation work. The height and width of the
channel were made equal (H = W), ranging from 0.4 mm to 3 mm. The length of the channel (L)
varied from 10 mm to 22.5 mm. The grid spacing was chosen to be 12.5 mm, giving H/h = 80 for
the smallest channel and 240 for the largest channel. The timestep used in the simulation was not
specified.
A truncated spherical bubble that satisfies a specified contact angle (j) was placed at the bottom
wall of the channel near the inlet. The initial diameter of the bubble was 5 mm. The static contact
angle of the bubble was varied from 30° to 50° with a dynamic contact line model used so that the
contact angle depended on whether the bubble was advancing or receding and on the contact line
speed. The flow was in x direction. Saturated water at 1 atm entered the inlet with a specified
velocity u0 = 0.12 m s-1. The no-slip condition was applied at the walls. A constant temperature
boundary condition was applied at the bottom wall with a superheating level ranging from 5 K to
20 K. The top and side walls were adiabatic. The fluids were assumed incompressible. A liquid-
only steady-state simulation was performed to provide the initial condition for the transient
analysis. The transient runs were terminated when the bubble touched the outlet boundary.
(a)
0.0
2.0
4.0
t
6.0
8.0
0.4
8.6
y
0
0 10
X
(b)
0
25
50
t
62
75
3
100
y
0
0 22.5
X
Figure 13. Bubble growth in microchannels at DT = 5 K and j = 40°. (a) H = 0.4 mm;
(b) H = 3 mm. The black solid lines are the isotherms. Taken from Lee and Son [41].

25
H = 0.4
H = 0.8
20 H = 3.0
15
q (W/cm2)
10
0
0 5 10 15 20 25
t
Figure 14. Effect of channel size on area-average wall heat flux at DT = 5 K and j = 40°.
Taken from Lee and Son [41].
The bubble growth process is displayed in Figure 13, and shows that the bubble in the smaller
channel (H = 0.4 mm) grew and touched the top wall very quickly and then transformed into an
elongated bubble, forming a thin liquid film between the bubble and the walls and also a large area
of dryout, which was in accord with the simulation results of Mukherjee and Kandlikar [63].
However, the bubble in the relatively larger channel (H = 3 mm) detached from the bottom wall
before it touched the top wall of the channel and hence no thin liquid film was formed. The
numerical results showed that the heat transfer in the smaller channels was much better than that in
the larger channels. As shown in Figure 14, the average wall heat flux obtained for the 0.4 mm
channel increases very fast from the start of the simulation, whereas the heat flux was
approximately constant during the simulation for the 3 mm channel. The dense isotherms shown in
Figure 13(a) suggested high heat fluxes across the thin liquid films between the bubble and the wall
in the 0.4 mm channel. These observations indicated that the formation of a thin liquid film might
play an important role in heat transfer of flow boiling in microchannels.
The effect of contact angle on the heat transfer was also investigated by the authors. Their results
showed that the average wall heat flux decreased with increasing contact angle and this is in
agreement with the results of Mukherjee et al. [64]. Similar to the discussion in Section 4.1, it is
suspected that the occurrence of the large area of dryout might be caused by poor mesh resolution.
In fact, in the simulation work on Taylor bubbles evaporation by Magnini
et al. [43], where the channel size was 0.5 mm and the wall superheats were close to this work, a
permanent liquid film was observed between the bubble and the tube wall during the simulation
period.
Lee et al. [39] extended this numerical model to simulate flow boiling in finned channels. The
above channel was adapted by placing parallel fins at the bottom wall of the channel. As shown in
Figure 15, Lfin and Hfin are the length and height of the fins and Sfin denotes the fin spacing. The
size of the domain was fixed at 4 mm, 0.32 mm and 0.48 mm in x, y and z directions, respectively.
Note, the domain includes the bottom and side walls of Hb = 0.12 mm and Ws = 0.04 mm. The grid
spacing was chosen to be 0.01 mm, giving 400, 32, 48 cells in x, y and z directions, respectively.
The boundary conditions were fixed in order to perform the parametric analysis for the fin sizes.
The bottom wall temperature was specified at Tw = 110 °C. Note, the authors included the heat
conduction across the wall in the simulation, using a wall material heat conductivity of kw = 230kl.
The inlet liquid was saturated water at 1atm with a specified velocity u0 = 0.3 m s-1. The contact

y
x
z
Side wall
Ws
Top x Fin
view Wc
z
Ws
Sfin Lfin
Front Bubble Hc
view Hfin
θ
y Hb
x x=L
Bottom wall
Figure 15. Computational domain of a finned channel. The actual numbers of fins
depend on the fin length and spacing. q is the contact angle. Taken from Lee et al. [39].
angle at the liquid-vapour-solid contact line was set to 40°. A bubble with a radius of 0.03 mm was
initiated at the bottom wall near the inlet. The initial condition for the transient analysis was
obtained from a liquid-only steady-state simulation.
The first analysis was the effect of fin height on the heat transfer of flow boiling in
microchannels. The fin spacing and the length was fixed as Lfin = Sfin = 0.1 mm. The simulations
were performed in channels with Hfin = 0.05 mm and 0.03 mm. To compare with the heat transfer
performance in plain channels, a simulation with no fins (Hfin= 0 mm) was also performed.
Figure 16 shows the flow maps at the symmetry plane (z = 0) for channels with different fin
sizes. It is clear that the bottom wall surface undergoes complete dryout in the channel with no fins
and therefore led to a significant decrease in heat transfer. Liquid films were formed in the gaps
between the fins in both the channels with fins. The liquid films in the channel with smaller fin
height (0.3 mm) were thinner than those in the channel with larger fin height (0.5 mm). As a
constant temperature wall condition was used, the thinner liquid film resulted in a larger
temperature gradient and hence led to a higher heat conduction rate through the film, as shown by
the dense isotherms in Figure 16(c). However, the dryout regions were still observed at the top
surfaces of the fins. The authors evaluated the heat transfer enhancement over the area of the
bottom wall surface and the total finned surface. Their results showed that the channel with a fin
height of 0.3 mm offered the best heat transfer enhancement.
Parametric studies for fin spacing and length were also performed. The results showed that the
heat transfer enhancement was weakly dependent on the fin spacing but increased significantly
with decreasing fin length. The authors concluded that the dryout areas at the top surfaces of the
fins, which had a negative effect on the heat transfer, were diminished by decreasing fin length.
The above work performed by Son and co-workers successfully simulated bubble growth in
conventional and finned microchannels. Their results showed that a suitable fin size could improve
the heat transfer of flow boiling in microchannels. Dryout areas were observed in all the cases
performed in their work, whereas a permanent liquid film was found to exist between the bubble
and the tube wall in the simulation work by Magnini et al. [43], where the wall superheats and the

(a) z=0 (b) z=0

0.32 0.32
0.16 0.16
y
y
109.9
0.00 0.00
2.5 3.0 3.5 2.5 3.0 3.5
x x
(c) z=0
0.32
0.16
y 0.00
2.5 3.0 3.5
x
Figure 16. Flow maps at the symmetry plane at t = 3 ms for various fin heights. (a) No
fins; (b) Hfin= 0.05 mm; (c) Hfin= 0.03 mm. Yellow represents the gas phase whereas
white denotes the liquid phase. The black solid lines in the liquid phase are isotherms.
The dotted lines are isotherms in the wall. Taken from Lee et al. [39].
tube sizes were close to those of Lee et al. [39]. The grid spacing used in Magnini et al. [43] was
less than one fifth of that used in this work. This indicates that the large dryout areas in this work
may be caused by poor mesh resolution so that the heat transfer performance of flow boiling in
finned channels still needs to be analysed.
4.5. Simulation Work on Taylor Bubble Evaporation by Thome’s Group

Magnini et al. [43] investigated the hydrodynamics and heat transfer of an elongated Taylor bubble
during flow boiling in microchannels by performing 2-D transient simulations in ANSYS Fluent
12. The gas-liquid interface was captured using the VOF method.
With reference to dimensional analyses performed by Morini [86] and Ong and Thome [87],
viscous heating and gravitational effects were neglected in this simulation work. The fluids were
assumed incompressible. Therefore the mass conservation equation can be expressed as in eqn (30),
.
where m is the evaporation mass flux of vapour across the liquid-vapour interface.
⎛1 1⎞
∇ ⋅ ν = ⎜⎜ − ⎟⎟ m ∇α (30)
⎝ ρv ρ l ⎠
The surface tension term in the momentum conservation equation was not modelled using the
default ANSYS Fluent default scheme to calculate interface curvature but was replaced by a height
function algorithm proposed by the authors with reference to Malik et al. [51] and Hernandez et al.
[52] in order to improve the estimation of curvature. The algorithm was implemented using user-
defined functions (UDF).
The evaporation model implemented in the simulations is a modified version of that of
Tanasawa [78] (Section 3.4.3) assuming that at the interface Tv ≈ Tsat(p∞), with p∞ being the system
pressure. Therefore, the evaporation mass flux was calculated as shown in eqns (25) and (26) with
the interface temperature replaced by the saturation temperature corresponding to the system
pressure. Note, the saturation temperature varies throughout the solution domain because of the
pressure drop along the tube and the Laplacian pressure jumps across the liquid-vapour interface.
However, the authors neglected the effect of these two terms based on estimates of their size. In
addition, the authors modified the source terms by means of a smoothing procedure proposed by
Hardt and Wondra [36] which was mentioned in Section 3.3.
Five simulation runs were performed, including for four different fluids and different conditions.
Figure 17 shows a schematic view of the computational domain and the initial and boundary
conditions for the first run.
A 2-D axisymmetric channel was modelled, as shown in Figure 17 above, with a tube diameter
D = 0.5 mm. The length of the domain is 20D. The domain is split to an adiabatic region of length
8D from the inlet and a heated region of length 12D to the outlet. With reference to Gupta et al.

8 4000
h (W/m2K)
Tw-Tsat (K)
6 3000
4 2000
2 1000
q=0 q = 9 kW/m2
R113 liquid
Vapor
G = 600 kg/m2s,
T = 323.15 K
0 2 4 6 8 10 12 14 16 18 20
z/D
0 1 2 3 4 5 6 7
T-Tsat (K)
Figure 17. Schematic view of simulation case 1. Heat flux is applied from z/D = 8. Taken
from Magnini et al. [43].
4000
Two-phase
3000 Single-phase
h (W/m2K)
2000
1000
0.45
r/D
0.4
0.35
0.3
8 10 12 14 16 18 20
z/D
0 1 2 3 4 5 6 7
T-Tsat (K)
Figure 18. Heat transfer coefficient and temperature field. The black solid line on the
temperature contour identifies the interface. Taken from Magnini et al. [43].
[88] and Han and Shikazono [89, 90], the maximum grid space (D) was determined as D/D = 100
in order to successfully capture the liquid film between the bubble and the tube wall. However, the
authors thought finer grids were needed here because with evaporation the temperature in the liquid
film could change faster than that without phase change, which was the case in Gupta et al. [88].
Therefore they performed a grid independence analysis and determined the grid size required in this
simulation to be D/D = 300. The timesteps were determined from the grid sizes, using the default
Courant number of 0.25.
As shown in Figure 17, a constant mass flux of 600 kg m-2 s-1 of R113 liquid is imposed at the
inlet. A vapor bubble of length 3D at saturation conditions is initialized as a cylinder with spherical
ends near the inlet. Zero gradient boundary condition for velocity and temperature are applied to
the outlet. A constant heat flux of 9 kW m-2 is applied at the heated walls.
A prior liquid-only steady state simulation was performed to obtain the initial condition for the
transient analysis. When the bubble enters the heated region, the liquid-vapour interface enters the
superheated thermal boundary layer created by the wall heat flux. Then the liquid at the interface
evaporates and the removal of latent heat cools down the interfacial cells to the saturation temperature.
Therefore the bubble grows and accelerates. The adiabatic region in the entry section of the tube is used

0.3
0.25
(htp-hsp)/hsp
0.2
0.15
0.1
0.05
0
8 10 12 14 16 18 20
z/D
Figure 19. Enhancement of heat transfer coefficient relative to the single phase value. The
dashed lines indicate the bubble nose and tail positions. Taken from Magnini et al. [43].
to allow the bubble to become a steady Taylor bubble before it enters the heated region.
The simulation results showed that the local heat transfer coefficient was increased when the
bubble passed through the heated region. Figure 18 depicts the heat transfer coefficients and
the temperature field around the bubble. Small waves were observed on the gas-liquid interface at
the bubble tail, that the author concluded to be capillary waves, referring to the experimental work
performed by Liberzon et al. [91]. It is shown in both Figure 18 and Figure 19 that the
enhancement of heat transfer coefficient relative to single phase flow increases with decreasing
film thickness, and this is in agreement with the observation of Gupta et al. [75] for Taylor flow
without phase change.
Magnini et al. [44] extended the above work to simulate two consecutive elongated bubbles and
investigated the influence of bubble interaction on flow boiling heat transfer in microchannels. In
this case, all the previous assumptions were maintained. The boundary conditions were not changed
except the adiabatic channel length was doubled in order to allow both of the bubbles to reach a
steady-state upstream of the heated section. The bubbles were initialized at the inlet with an
arbitrarily selected separation distance of 6D.
The results showed that the leading bubble had a larger acceleration and higher bubble speed.
This may be because the leading bubble reduces the temperature of the super heated liquid region
and there is not enough time for the system to restore the steady temperature field when the two
bubbles are close to each other (6D separation). Meanwhile, the leading bubble head was observed
to be less blunt than that of the trailing bubble. These two observations indicate that the simulation
is far from fully-developed. Actually, the heat transfer coefficients observed at a fixed axial position
are higher when the second bubble is passing than that observed when the leading bubble is
passing. Therefore, more bubbles are needed to obtain a fully-developed flow field and temperature
field so that the flow boiling mechanism in micro channels can be explicitly examined.
The velocity field in the liquid slug between the two bubbles showed a variation of axial velocity
in the radial direction which gave rise to an internal recirculating flow in the slug, as shown in
Figure 20. In this recirculation zone, the heated liquid near the wall was brought to the centre of
the slug and therefore the heat transfer was enhanced compared with that of the liquid-only steady
state flow. Gupta et al. [75] obtained a similar recirculation zone in Taylor flow heat transfer
simulations without phase change and their calculations showed that the Nusselt number in the slug
was ~2.5 times higher than that for the fully-developed single phase laminar flow with the same
flow conditions. Qualitatively similar temperature and velocity fields can also be found in the work
of Fukagata et al. [92] and Lakehal et al. [93].

(a)
0.5
0.4
r/D 0.3
0.2
0.1
0
24 25 26 27 28 29 30 31 32 33 34
z/D
(b)
0.5
0.4
0.3
r/D
0.2
0.1
0
24 25 26 27 28 29 30 31 32 33 34
z/D
Figure 20. Streamlines of the (a) velocity field and (b) relative velocity field (u-Ur,N)
where Ur,N is the velocity of the nose of the trailing bubble. The blue lines identify the
bubble locations. Taken from Magnini et al. [44].
For comparison, Magnini et al. [43, 44] firstly used Thome’s three zone model [27] (Section 2.3)
to predict the above simulation results for both a single bubble and two bubbles, however the
numerical results were very different from the predictions. The authors believed that this
disagreement was because Thome’s model was based on the assumption of a quasi-steady-state
condition which was not the case in this simulation. Hence, the authors introduced a new non-CFD
model considering the transient heat conduction in the film to predict the results of the current
work. They showed that the new model predicted the numerical results with a positive error of up
to 18% [44]. However, the input parameters used in the theoretical model, such as the film
thickness, were obtained from the simulation results so this transient model is not an independent
model.
The above work performed by Magnini et al. [43, 44] simulated flow boiling of elongated Taylor
bubbles, utilizing a kinetic theory based evaporation model and a smoothing procedure proposed
by Hardt and Wondra [36] to handle the evaporation. A height function for curvature estimation was
introduced to replace ANSYS Fluent default scheme. For the simulation itself, one important fact
is that the flow fields obtained in the simulations are not fully-developed and therefore the boiling
heat transfer mechanism cannot be fully elucidated.
5. CONCLUSIONS AND RECOMMENDATIONS

A large number of experimental and theoretical studies of microchannel flow boiling have been
performed for decades and an increasing number of numerical studies have been reported in recent
years. Although some agreement between experimental data, the theoretical calculations and the
computational results has been obtained, debate continues and gaps in microchannel flow boiling
research remain. Nevertheless, some conclusions can be made from the above review.
• The available experimental data and analysis have not always provided a clear picture from
which modelling studies can be started. However, there now seems to be a consensus that
nucleation leads to the formation of large bubbles that rapidly fill the tube, leading to Taylor
flow. It seems likely under certain conditions this elongated bubble regime can also undergo
transition to annular flow but no such studies have been performed.
• A limited number of analytical models have been developed for the prediction of micro channel
flow boiling. Amongst them, Thome’s three-zone model [27] is considered to be the most
developed so far. However, this model achieved good agreements with experimental data from

different research groups by applying different parameters obtained by least-squares fitting and
therefore it is not able to provide general predictions for microchannel flow boiling. It takes no
account of inlet compressibility, even though some datasets have subsequently been shown to
have had highly compressible inlet conditions.
• An increasing number of CFD studies of microchannel flow boiling have been performed in
the past few years. These solve the fundamental equations for conservation of mass,
momentum and energy using either commercial software or purposely written codes.
• As flow boiling simulations require the integration of a phase change model with liquid-
vapour boundary determination, interface capturing methods have been used in the CFD
studies for this purpose. The volume of fluid (VOF) method is the most popular among the
interface capturing methods used because of its simplicity and wide availability in
commercial CFD codes, although several groups have used the level-set method in an attempt
to get better interface capturing and curvature for use in the surface tension force calculation.
• In the reported numerical studies homogenous velocity and temperature fields have been
assumed. The effect of allowing local dis-equilibrium has not been studied. However, if a
sharp interface is postulated to exist, as opposed to a dispersion of droplets or bubbles, the
homogeneous approach would seem to be the most suitable, especially as it is the least
computationally expensive approach.
• The treatment of phase change in the available computational studies varies considerably,
ranging from the use of experimental data based models to kinetic rate limited models. The
most recent CFD work uses kinetic rate limited models as they better reflect the physics
happening during the flow boiling processes.
• The mass transfer rate at the liquid-gas interface is high in flow boiling simulations.
Smoothing and source shifting procedures for interface mass transfer need to be considered
in flow boiling simulations if numerical instabilities are to be avoided.
• Most of the current simulations are performed using two-dimensional computational
domains. This is almost certainly due to the combination of small cell sizes and short
timesteps needed in the simulations leading to very long run times. Whilst authors are not
reporting information on computational requirements, the relatively low number of
simulations reported in a study is evidence of the huge computational costs.
• Implicit in the use of 2D simulations is the assumption that gravitational effects are not
important. One set of simulations has explicitly shown this to be the case. This highlights the
fact that microchannel boiling correlations in which gravitational effects are included, are
probably not correlating the correct physics.
• The modelling strategy varies in the various studies. Some of the simulations started with a
pre-defined bubble and investigated the growth of the bubble and the following film
evaporation process. Other studies began with bubble formation and investigated the bubble
coalescence and the flow regime transition process. This is again a reflection of numerical
limitations. Mesh sizes cannot be fine enough to track bubble nucleation events and
simulations are so costly that initial conditions must be chosen to give results of interest in
a practical computational time. However, it would seem sensible to adopt such approaches
as the details of the nucleation process are not important in the simulation of the bubble
dynamics and heat transfer once Taylor bubbles are formed.
ACKNOWLEDGEMENTS
The work was supported under Australian Research Council Discovery Grant DP120103235.
Z. Guo also acknowledges the China Scholarship Council and the University of Sydney.
REFERENCES
[1] C.B. Sobhan, S.V. Garimella, A comparative analysis of studies on heat transfer and fluid flow in microchannels,
Microscale Thermophysical Engineering, 2001, 5(4), 293–311.
[2] B. Agostini, M. Fabbri, J.E. Park, L. Wojtan, J.R. Thome, B. Michel, State of the art of high heat flux cooling
technologies, Heat Transfer Engineering, 2007, 28(4), 258–281.
[3] S.G. Kandlikar, Fundamental issues related to flow boiling in minichannels and microchannels, Experimental
Thermal and Fluid Science, 2002, 26(2–4), 389–407.

[4] S.V. Garimella, Advances in mesoscale thermal management technologies for microelectronics, Microelectronics
Journal, 2006, 37(11), 1165–1185.
[5] S.S. Bertsch, E.A. Groll, S.V. Garimella, Review and comparative analysis of studies on saturated flow boiling in
small channels, Nanoscale and Microscale Thermophysical Engineering, 2008, 12(3), 187–227.
[6] J.R. Thome, Boiling in micro channels: a review of experiment and theory, International Journal of Heat and Fluid
Flow, 2004, 25(2), 128–139.
[7] S.V. Garimella, C.B. Sobhan, Transport in microchannels - a critical review, Annual Review of Heat Transfer, 2003,
13(13), 1–50.
[8] S.G. Kandlikar, History, advances, and challenges in liquid flow and flow boiling heat transfer in microchannels: a
critical review, Journal of Heat Transfer, 2012, 134(3), 034001.
[9] C. Baldassari, M. Marengo, Flow boiling in microchannels and microgravity, Progress in Energy and Combustion
Science, 2013, 39(1), 1–36.
[10] Z.Y. Bao, D.F. Fletcher, B.S. Haynes, Flow boiling heat transfer of Freon R11 and HCFC123 in narrow passages,
International Journal of Heat and Mass Transfer, 2000, 43(18), 3347–3358.
[11] T.N. Tran, M.W. Wambsganss, D.M. France, Small circular- and rectangular-channel boiling with two refrigerants,
International Journal of Multiphase Flow, 1996, 22(3), 485–498.
[12] W. Yu, D.M. France, M.W. Wambsganss, J.R. Hull, Two-phase pressure drop, boiling heat transfer, and critical
heat flux to water in a small-diameter horizontal tube, International Journal of Multiphase Flow, 2002, 28(6),
927–941.
[13] R. Mertz, A. Wein, M. Groll, Experimental investigation of flow boiling heat transfer in narrow channels, Heat and
Technology, 1996, 14(2), 47–54.
[14] R.J. Armstrong, The temperature difference in nucleate boiling, International Journal of Heat and Mass Transfer,
1966, 9(10), 1148–1149.
[15] M.G. Cooper, Heat flow rates in saturated nucleate pool boiling-a wide-ranging examination using reduced
properties, in: P.H. James, F.I. Thomas, eds. Advances in Heat Transfer, Elsevier, 1984, 157–239.
[16] M.K. Akbar, D.A. Plummer, S.M. Ghiaasiaan, On gas–liquid two-phase flow regimes in microchannels,
International Journal of Multiphase Flow, 2003, 29(5), 855–865.
[17] J. Lee, I. Mudawar, Two-phase flow in high-heat-flux micro-channel heat sink for refrigeration cooling applications:
Part II—heat transfer characteristics, International Journal of Heat and Mass Transfer, 2005, 48(5), 941–955.
[18] X.F. Peng, H.Y. Hu, B.X. Wang, Flow boiling through V-shape microchannels, Experimental Heat Transfer, 1998,
11(1), 87–100.
[19] S.S. Bertsch, E.A. Groll, S.V. Garimella, Refrigerant flow boiling heat transfer in parallel microchannels as a
function of local vapor quality, International Journal of Heat and Mass Transfer, 2008, 51(19), 4775–4787.
[20] S. Lin, P.A. Kew, K. Cornwell, Two-phase heat transfer to a refrigerant in a 1 mm diameter tube, International
Journal of Refrigeration, 2001, 24(1), 51–56.
[21] R. Yun, J. Hyeok Heo, Y. Kim, Evaporative heat transfer and pressure drop of R410A in microchannels, International
Journal of Refrigeration, 2006, 29(1), 92–100.
[22] H. J. Lee, S.Y. Lee, Heat transfer correlation for boiling flows in small rectangular horizontal channels with low
aspect ratios, International Journal of Multiphase Flow, 2001, 27(12), 2043–2062.
[23] W. Qu, I. Mudawar, Flow boiling heat transfer in two-phase micro-channel heat sinks––II. Annular two-phase flow
model, International Journal of Heat and Mass Transfer, 2003, 46(15), 2773–2784.
[24] W. Qu, I. Mudawar, Flow boiling heat transfer in two-phase micro-channel heat sinks––I. Experimental investigation
and assessment of correlation methods, International Journal of Heat and Mass Transfer, 2003, 46(15), 2755–2771.
[25] Y. Liu, D.F. Fletcher, B.S. Haynes, On the importance of upstream compressibility in microchannel boiling heat
transfer, International Journal of Heat and Mass Transfer, 2013, 58(1-2), 503–512.
[26] A.M. Jacobi, J.R. Thome, Heat transfer model for evaporation of elongated bubble flows in microchannels, Journal
of Heat Transfer, 2002, 124(6), 1131–1136.
[27] J.R. Thome, V. Dupont, A.M. Jacobi, Heat transfer model for evaporation in microchannels. Part I: presentation of
the model, International Journal of Heat and Mass Transfer, 2004, 47(14–16), 3375–3385.
[28] V. Dupont, J.R. Thome, A.M. Jacobi, Heat transfer model for evaporation in microchannels. Part II: comparison with
the database, International Journal of Heat and Mass Transfer, 2004, 47(14–16), 3387–3401.

[29] D. Brutin, F. Topin, L. Tadrist, Experimental study of unsteady convective boiling in heated minichannels,
International Journal of Heat and Mass Transfer, 2003, 46(16), 2957–2965.
[30] G.M. Lazarek, S.H. Black, Evaporative heat-transfer, pressure-drop and critical heat-flux in a small vertical tube with
R-113, International Journal of Heat and Mass Transfer, 1982, 25(7), 945–960.
[31] M.W. Wambsganss, J.A. Jendrzejczyk, T.N. Tran, D.M. France, Boiling heat transfer in a horizontal small-diameter
tube, Journal of Heat Transfer (Transactions of the ASME (American Society of Mechanical Engineers), Series C);,
1993, 115(4), 963–972.
[32] Y.-Y. Yan, T.-F. Lin, Evaporation heat transfer and pressure drop of refrigerant R-134a in a small pipe, International
Journal of Heat and Mass Transfer, 1998, 41(24), 4183–4194.
[33] J.R. Baird, Z.Y. Bao, D.F. Fletcher, B.S. Haynes, Local flow boiling heat transfer coefficients in narrow conduits,
Multiphase Science and Technology, 2000, 12(3&4), 16.
[34] B. Agostini, Etude expérimentale de l’ébullition de fluide réfrigérant en convection forcée dans des mini-canaux,
Ph.D Thesis, Universite Joseph Fourier, 2002.
[35] R. Gupta, D.F. Fletcher, B.S. Haynes, Taylor flow in microchannels: a review of experimental and computational
work, The Journal of Computational Multiphase Flows, 2010, 2(1), 1–32.
[36] S. Hardt, F. Wondra, Evaporation model for interfacial flows based on a continuum-field representation of the source
terms, Journal of Computational Physics, 2008, 227(11), 5871–5895.
[37] Y. Egorov, ANSYS Germany, Private communication, 2013.
[38] M.A. Grolmes, H.K. Fauske, Propagation characteristics of compression and rarefaction pressure pulses in one-
component vapor-liquid mixtures, Nuclear Engineering and Design, 1970, 11(1), 137–142.
[39] W. Lee, G. Son, H.Y. Yoon, Direct numerical simulation of flow boiling in a finned microchannel, International
Communications in Heat and Mass Transfer, 2012, 39(9), 1460–1466.
[40] Y. Suh, W. Lee, G. Son, Bubble dynamics, flow, and heat transfer during flow boiling in parallel microchannels,
Numerical Heat Transfer, Part A: Applications, 2008, 54(4), 390–405.
[41] W. Lee, G. Son, Bubble dynamics and heat transfer during nucleate boiling in a microchannel, Numerical Heat
Transfer, Part A: Applications, 2008, 53(10), 1074–1090.
[42] S.V. Patankar, Numerical heat transfer and fluid flow, Taylor & Francis, 1980.
[43] M. Magnini, B. Pulvirenti, J.R. Thome, Numerical investigation of hydrodynamics and heat transfer of elongated
bubbles during flow boiling in a microchannel, International Journal of Heat and Mass Transfer, 2013, 59, 451–471.
[44] M. Magnini, B. Pulvirenti, J.R. Thome, Numerical investigation of the influence of leading and sequential bubbles
on slug flow boiling within a microchannel, International Journal of Thermal Sciences, 2013, 71(0), 36–52.
[45] L.D. Landau, E.M. Lifshitz, Fluid mechanics, 2nd, in, London: Elsevier, 1987.
[46] S. Sugden, VI.—The variation of surface tension with temperature and some related functions, Journal of the
Chemical Society, Transactions, 1924, 125, 32–41.
[47] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, Journal of Computational
Physics, 1992, 100(2), 335–354.
[48] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, G. Zanetti, Modelling merging and fragmentation in multiphase
flows with SURFER, Journal of Computational Physics, 1994, 113(1), 134–147.
[49] D.J.E. Harvie, M.R. Davidson, M. Rudman, An analysis of parasitic current generation in Volume of Fluid
simulations, Applied Mathematical Modelling, 2006, 30(10), 1056–1066.
[50] D.F. Fletcher, B.S. Haynes, J. Aubin, C. Xuereb, Modelling of microfluidic devices, in: V. Hessel, J.C. Schouten, A.
Renken, J.-I. Yoshida, eds. Handbook of Micro Reactors. Fundamentals, Operations and catalysts, Wiley-VCH
2009, 117–144.
[51] M. Malik, E.S.-C. Fan, M. Bussmann, Adaptive VOF with curvature-based refinement, International Journal for
Numerical Methods in Fluids, 2007, 55(7), 693–712.
[52] J. Hernández, J. López, P. Gómez, C. Zanzi, F. Faura, A new volume of fluid method in three dimensions—Part I:
Multidimensional advection method with face-matched flux polyhedra, International Journal for Numerical
Methods in Fluids, 2008, 58(8), 897–921.
[53] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of
Computational Physics, 1981, 39(1), 201–225.
[54] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, Journal of Computational Physics, 1998, 141(2), 112–152.

[55] M. Rudman, Volume-tracking methods for interfacial flow calculations, International Journal for Numerical
Methods in Fluids, 1997, 24(7), 671–691.
[56] D. Youngs, Time-dependent multi-material flow with large fluid distortion, Numerical Methods for Fluid Dynamics,
1982, 24, 273–285.
[57] O. Ubbink, R.I. Issa, A method for capturing sharp fluid interfaces on arbitrary meshes, Journal of Computational
Physics, 1999, 153(1), 26–50.
[58] R. Zhuan, W. Wang, Flow pattern of boiling in micro-channel by numerical simulation, International Journal of Heat
and Mass Transfer, 2012, 55(5–6), 1741–1753.
[59] Y.Q. Zu, Y.Y. Yan, S. Gedupudi, T.G. Karayiannis, D.B.R. Kenning, Confined bubble growth during flow boiling in
a mini-/micro-channel of rectangular cross-section part II: Approximate 3-D numerical simulation, International
Journal of Thermal Sciences, 2011, 50(3), 267–273.
[60] R. Zhuan, W. Wang, Simulation on nucleate boiling in micro-channel, International Journal of Heat and Mass
Transfer, 2010, 53(1–3), 502–512.
[61] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi
formulations, Journal of Computational Physics, 1988, 79(1), 12–49.
[62] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow,
Journal of Computational Physics, 1994, 114(1), 146–159.
[63] A. Mukherjee, S.G. Kandlikar, Numerical simulation of growth of a vapor bubble during flow boiling of water in a
microchannel, Microfluidics and Nanofluidics, 2005, 1(2), 137–145.
[64] A. Mukherjee, S.G. Kandlikar, Z.J. Edel, Numerical study of bubble growth and wall heat transfer during flow
boiling in a microchannel, International Journal of Heat and Mass Transfer, 2011, 54(15–16), 3702–3718.
[65] S. Zhou, X. Xu, B.G. Sammakia, Modeling of boiling flow in microchannels for nucleation characteristics and
performance optimization, International Journal of Heat and Mass Transfer, 2013, 64(0), 706–718.
[66] A. Mukherjee, V.K. Dhir, Study of lateral merger of vapor during nucleate pool boiling, Journal of Heat Transfer-
Transactions of the ASME, 2004, 126(6), 1023–1039.
[67] R.R. Nourgaliev, S. Wiri, N.T. Dinh, T.G. Theofanous, On improving mass conservation of level set by reducing
spatial discretization errors, International Journal of Multiphase Flow, 2005, 31(12), 1329–1336.
[68] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, Journal of
Computational Physics, 1992, 100(1), 25–37.
[69] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.J. Jan, A front-tracking
method for the computations of multiphase flow, Journal of Computational Physics, 2001, 169(2), 708–759.
[70] M.W. Baltussen, J.A.M. Kuipers, N.G. Deen, A critical comparison of surface tension models for the volume of fluid
method, Chemical Engineering Science, 2014, 109, 65-74.
[71] D.M. Anderson, G.B. McFadden, A.A. Wheeler, Diffuse-interface methods in fluid mechanics, Annual Review of
Fluid Mechanics, 1998, 30(1), 139–165.
[72] M. Sussman, E.G. Puckett, A coupled level set and volume-of-fluid method for computing 3D and axisymmetric
incompressible two-phase flows, Journal of Computational Physics, 2000, 162(2), 301–337.
[73] N. Zhang, D.F. Chao, Models for enhanced boiling heat transfer by unusual Marangoni effects under microgravity
conditions, International Communications in Heat and Mass Transfer, 1999, 26(8), 1081–1090.
[74] S. Gedupudi, Y.Q. Zu, T.G. Karayiannis, D.B.R. Kenning, Y.Y. Yan, Confined bubble growth during flow boiling in
a mini/micro-channel of rectangular cross-section Part I: Experiments and 1-D modelling, International Journal of
Thermal Sciences, 2011, 50(3), 250–266.
[75] R. Gupta, D.F. Fletcher, B.S. Haynes, CFD modelling of flow and heat transfer in the Taylor flow regime, Chemical
Engineering Science, 2010, 65(6), 2094–2107.
[76] A. Mukherjee, Contribution of thin-film evaporation during flow boiling inside microchannels, International Journal
of Thermal Sciences, 2009, 48(11), 2025–2035.
[77] R.W. Schrage, A theoretical study of interphase mass transfer, Columbia University Press, 1953.
[78] I. Tanasawa, Advances in condensation heat transfer, Advances in Heat Transfer, 1991, 21, 55–139.
[79] W.M. Rohsenow, Film condensation of liquid metals, ICHMT Digital Library Online, 1988.
[80] R.P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach to interfaces in multimaterial
flows (the ghost fluid method), Journal of Computational Physics, 1999, 152(2), 457–492.

[81] S.G. Kandlikar, P. Balasubramanian, Effect of gravitational orientation on flow boiling of water in 1054 ¥ 197 μm
parallel minichannels, Journal of Heat Transfer, 2004, 127, 820–829.
[82] J. Xu, W. Zhang, G. Liu, Seed bubble guided heat transfer in a single microchannel, Heat Transfer Engineering,
2011, 32(11–12), 1031–1036.
[83] D. Liu, P.-S. Lee, S.V. Garimella, Prediction of the onset of nucleate boiling in microchannel flow, International
Journal of Heat and Mass Transfer, 2005, 48(25–26), 5134–5149.
[84] L. Consolini, Convective boiling heat transfer in a single micro-channel, M. Sc. Thesis, University of Illinois, 2008.
[85] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and heat transfer associated with a single bubble during nucleate
boiling on a horizontal surface, Journal of Heat Transfer, 1999, 121(3), 623–631.
[86] G.L. Morini, Viscous heating in liquid flows in micro-channels, International Journal of Heat and Mass Transfer,
2005, 48(17), 3637–3647.
[87] C.L. Ong, J.R. Thome, Macro-to-microchannel transition in two-phase flow: Part 1 – Two-phase flow patterns and
film thickness measurements, Experimental Thermal and Fluid Science, 2011, 35(1), 37–47.
[88] R. Gupta, D.F. Fletcher, B.S. Haynes, On the CFD modelling of Taylor flow in microchannels, Chemical Engineering
Science, 2009, 64(12), 2941–2950.
[89] Y. Han, N. Shikazono, Measurement of the liquid film thickness in micro tube slug flow, International Journal of
Heat and Fluid Flow, 2009, 30(5), 842–853.
[90] Y. Han, N. Shikazono, The effect of bubble acceleration on the liquid film thickness in micro tubes, International
Journal of Heat and Fluid Flow, 2010, 31(4), 630–639.
[91] D. Liberzon, L. Shemer, D. Barnea, Upward-propagating capillary waves on the surface of short Taylor bubbles,
Physics of Fluids, 2006, 18(4), 048103.
[92] K. Fukagata, N. Kasagi, P. Ua-arayaporn, T. Himeno, Numerical simulation of gas–liquid two-phase flow and
convective heat transfer in a micro tube, International Journal of Heat and Fluid Flow, 2007, 28(1), 72–82.
[93] D. Lakehal, G. Larrignon, C. Narayanan, Computational heat transfer and two-phase flow topology in miniature
tubes, Microfluidics and Nanofluidics, 2008, 4(4), 261–271.

Guo Et Al 2014 A Review of Computational Modelling of Flow Boiling in Microchannels

Uploaded by

Copyright:

Available Formats

Guo Et Al 2014 A Review of Computational Modelling of Flow Boiling in Microchannels

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Guo Et Al 2014 A Review of Computational Modelling of Flow Boiling in Microchannels

Uploaded by

Copyright:

Available Formats

79

A Review of Computational Modelling of

Recieved: 3 December 2013; Accepted: 1 March 2014

2. EXPERIMENTAL AND ANALYTICAL WORK ON MICROCHANNEL FLOW BOILING

*Corresponding author: E-mail: david.fletcher@sydney.edu.au

2.1 Discussion of the dominant mechanism of flow boiling heat transfer

Journal of Computational Multiphase Flows

2.2 Effects of upstream compressibility

2.3 Thome’s three-zone model

Volume 6 · Number 2 · 2014

3. CFD FRAMEWORK OF FLOW BOILING

Journal of Computational Multiphase Flows

Table 1. Values of parameters in optimization function for Thome’s three-zone model

Study Fluid Diameter Specific parameters for each study

The parameter ranges were added by the authors.

3.1 Governing Equations

where τ i = μi ( ∇vi + ∇vi )

Volume 6 · Number 2 · 2014

Svg = Sρ g vg Sρ g < 0 (7a)

Svg = Sρ g vl Sρ g > 0 (7b)

Journal of Computational Multiphase Flows

Shl = −Sρ g ( hl + hlv ) (8a)

Shg = Sρ g ( hl + hlv ) (8b)

where hlv is the latent heat of vaporisation.

Shl = −Sρ g hg (9b)

3.1.1. Equations of state

3.2 Boundary conditions

3.2.1. Wall conditions

3.2.2. Interface conditions

Volume 6 · Number 2 · 2014

3.3 Numerical Methods for Flow Boiling Simulation

Journal of Computational Multiphase Flows

Fσ = σκδ ( r − rint ) n (20)

3.3.1. Volume of fluid (VOF) method

Volume 6 · Number 2 · 2014

3.3.2. Level-set method

3.3.3. Marker point method

Journal of Computational Multiphase Flows

3.4 Phase Change Models

3.4.2. Phase change models based on heat conduction

qg = n. kg ∇Tg m = (ql −qg )/hlv

Volume 6 · Number 2 · 2014

n ⋅ (kl ∇Tl - kg ∇Tg )

3.4.3. Kinetic-based phase change models

Journal of Computational Multiphase Flows

4. NUMERICAL WORK FOR MICROCHANNEL FLOW BOILING

4.1 Bubble Dynamics Simulations by Kandlikar and co-workers

Volume 6 · Number 2 · 2014

Journal of Computational Multiphase Flows

Volume 6 · Number 2 · 2014

Dimensions (mm) 0.7

(a) 0.4 (b) 0.4

Equivalent diameter (mm)

(c) (d) 0.4

Journal of Computational Multiphase Flows

Volume 6 · Number 2 · 2014

simulated successfully during flow boiling in a rectangular microchannel. However, two

Journal of Computational Multiphase Flows

4.3 Simulations by the Brunel University Group