Numerical study of bubble growth and wall heat transfer during flow boiling in a microchannel

International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

Contents lists available at ScienceDirect

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

Numerical study of bubble growth and wall heat transfer during flow boiling in a microchannel A. Mukherjee a,⇑, S.G. Kandlikar b, Z.J. Edel a a b

Department of Mechanical Engineering and Engineering Mechanics, Michigan Technological University, 1400 Townsend Drive, Houghton, MI 49931, USA Department of Mechanical Engineering, Rochester Institute of Technology, 76, Lomb Memorial Drive, Rochester, NY 14623, USA

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 5 February 2010 Received in revised form 27 January 2011 Accepted 27 January 2011 Available online 12 April 2011

A numerical study has been performed to analyze the wall heat transfer mechanisms during growth of a vapor bubble inside a microchannel. The microchannel is of 200 lm square cross section and a vapor bubble begins to grow at one of the walls, with liquid coming in through the channel inlet. The complete Navier-Stokes equations along with continuity and energy equations are solved using the SIMPLER method. The liquid vapor interface is captured using the level set technique. Experiments have been conducted to validate the numerical model. The bubble growth rate and shapes show good agreement between numerical and experimental results. The numerical results show that the wall heat transfer increases with wall superheat but stays almost unaffected by the liquid flow rate. The liquid vapor surface tension value has little influence on bubble growth and wall heat transfer. However, the bubble with the lowest contact angle resulted in the highest wall heat transfer. Ó 2011 Elsevier Ltd. All rights reserved.

Keywords: Flow boiling Microchannels Bubbles

1. Introduction Flow through microchannels is a subject of extensive study due to wide ranging applications in engineering and biological sciences. Microchannel heat sinks with liquid cooling are extensively used in various applications such as electronic chip cooling. Bubble formation inside microchannels can take place if the fluid is a mixture of a gas and a liquid or the temperature of the wall reaches above the local liquid saturation temperature. During flow boiling, bubbles nucleate on the microchannel walls and may grow big enough to fill up the entire channel cross-section. When the bubbles are of the same size as the microchannel, they regulate the flow characteristics and if applicable the wall heat transfer. The wall heat transfer from the channel wall to the liquid is affected by the bubble nucleation and growth inside the channels. At microscale the surface tension forces are expected to dominate the gravitational forces and control the bubble dynamics. Water has been the preferred coolant for flow boiling for its excellent thermal properties. Direct application of boiling water on a chip surface may not be desirable due to its poor dielectric properties and high boiling temperature. Fluorochemicals with excellent dielectric properties are also often used for electronics cooling. The heat transfer mechanism in microchannels is affected not only by the thermal properties but also by the contact angle between the bubbles and the microchannel walls. The bubble nucleation tem⇑ Corresponding author. Tel.: +1 906 487 1174; fax: +1 906 487 2822. E-mail addresses: [email protected] (S.G. Kandlikar), [email protected] (Z.J. Edel).

(A.

Mukherjee),

[email protected]

0017-9310/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2011.01.030

perature is also known to be dependent on the contact angle. Dielectric liquids typically are highly wetting in nature with much lower contact angle and surface tension as compared to water.

2. Literature review Excellent reviews are available on flow boiling in microchannels, e.g. Bergles et al. [1], Garimella and Sobhan [2] and Thome [3]. Kandlikar [4] carried out a critical literature review on flow boiling through minichannels and microchannels. He identified the effect of surface tension to be significant in microchannels causing the liquid to form small uniformly spaced slugs that fill the channel. He pointed out that one of the main reasons for the boiling instability in microchannels is the explosive growth of a vapor bubble after it nucleates. It is therefore important to understand the dynamic growth characteristics of a bubble upon its nucleation during flow boiling in microchannels. Yen et al. [5] studied convective flow boiling in a circular pyrex glass microtube and a square pyrex glass microchannel. Higher heat transfer coefficient was observed in the square microchannel as compared to the circular crossectional microtube because of square corners acting as active nucleation sites. Agostini et al. [6–8] studied high heat flux flow boiling of refrigerants in microchannel heat sinks in a three-part paper. The authors compared the heat transfer results with their own theoretical models. Lee and Pan [9] studied eruptive boiling in silicon based microchannels and argued that it may be caused by the nano-sized cavities present at the channel walls.

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3703

Nomenclature A Cp D D G G H h hfg k K1 K2 L L1 L2 l0 m ms Nu p q Re T DT t t0 u u0 v w

wall area specific heat at constant pressure characteristic dimension grid spacing mass flux gravity vector heaviside function heat transfer coefficient latent heat of evaporation thermal conductivity 2 qL Kandlikar number 1 Ghq qG fg  2 Kandlikar number 2 hq rDq fg G length of bubble upstream bubble cap location downstream bubble cap location length scale mass transfer rate at interface milliseconds Nusselt number pressure heat flux Reynolds number temperature superheat, Tw  Tsat time time scale x direction velocity velocity scale y direction velocity z direction velocity

Bertsch et al. [10] measured local flow boiling heat transfer coefficient in a microchannel-based cold plate evaporator using HFC-134a. The heat transfer coefficient showed a peak value at 0.2 vapor quality in all of the experiments. Lee and Garimella [11] investigated flow boiling of water in a microchannel array and presented new correlations for predicting two-phase pressure drop and local saturated boiling heat transfer coefficient. Wang and Cheng [12] investigated subcooled flow boiling and microbubble emission boiling (MEB) of water in a microchannel. The occurrence of MEB resulted in removal of high heat flux at moderate rise in wall temperatures. Geisler and Bar-Cohen [13] studied saturated flow boiling CHF (Critical Heat Flux) in narrow vertical microchannels and observed that bubble confinement led to heat transfer enhancement in the low heat flux region of the nucleate boiling curve. Harirchian and Garimella [14] defined a new transition criterion to qualify microscale two-phase flow. They termed this number as the ‘convective confinement number’ that incorporated mass flux, fluid properties and channel cross-sectional area. Yang et al. [15] simulated bubbly two phase flow in a narrow channel using a numerical code FlowLab based on the LatticeBoltzmann method. Single or multiple two-dimensional Taylor bubbles were placed in a vertical channel and their behavior was studied for different values of surface tension and body forces. No heat transfer or phase change was considered between the two phases. The authors found little effect of surface tension on the movement of the bubbles or the flow regime transition. Jacobi and Thome [16] developed an analytical model of elongated bubble flows in microchannels and compared the results with experimental data. The central idea to this model was the thin film evaporation around the elongated bubbles. The model correctly predicted the heat transfer coefficient to be dependent on

x y z

distance in x direction distance in y direction distance in z direction

Greek symbols bT coefficient of thermal expansion j interfacial curvature l dynamic viscosity m kinematic viscosity q density r surface tension s time period / level set function u contact angle Subscripts evp evaporation G gas or vapor in inlet L, l liquid sat saturation v vapor w wall x o/ox y o/oy z o/oz Superscripts ⁄ non-dimensional quantity ? vector quantity

heat flux but insensitive to mass flux. The model did not include the possibility of formation of vapor patches at the walls and thus completely excluded the effect of contact angle. The success of this model depended exclusively on initial guesses of the critical nucleation radius and initial liquid film thickness. The authors argued that the experimental studies that show heat-flux dependence of the convection coefficient along with relative independence from quality and mass flux cannot be ascribed only to the nucleate boiling mechanism. They developed a hypothesis that microchannel evaporation is thin-film dominated. There are a few studies that have appeared in the literature on the effect of contact angle on vapor bubble growth inside microchannels. Mukherjee and Mudawar [17] studied microchannel electronics cooling using dielectric coolant FC-72 and water which has large differences in surface tension and contact angle. The results showed opposite trends in critical heat flux (CHF) with decrease in channel size for the two liquids. Water produced large diameter bubbles that blocked liquid flow through the channel cross-section and hence CHF decreased with decrease in channel gap below the departure diameter of the bubbles. In the case of the dielectric liquid, the departure bubble diameters were much smaller as compared to water and these bubbles blocked the liquid flow only below a channel gap of 0.13 mm. When the channel gap varied between 3.56 mm and 0.13 mm the CHF increased with a decrease in channel diameter because partial blockage caused by the tiny bubbles increased the two-phase mixture velocity and improved heat transfer. Thome et al. [18] and Dupont et al. [19] developed a three-zone flow boiling heat transfer model to describe evaporation of elongated bubbles in a microchannel and compared the time-averaged local heat transfer coefficients from several independent experimental studies. Their numerical model consisted of sequential

3704

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

and cyclic passage of (a) a liquid slug, (b) an evaporating elongated bubble with a thin liquid film around it and (c) a vapor slug, through a microchannel. The model had three adjustable parameters, the initial thickness of the liquid film, the minimum thickness of the liquid film at dry-out, and the bubble departure frequency. The comparison of the results with experimental data indicated limited success of the model. Kandlikar [20] observed that at low Reynolds number, the convective boiling is diminished in microchannels and nucleate boiling plays a major role with periodic flow of liquid and vapor slugs in rapid succession. He compared the transient conduction under the approaching rewetting liquid slug and the heat transfer in the evaporating meniscus region of the liquid–vapor interfaces in the contact line region to the nucleate boiling phenomenon. The author also listed several non-dimensional groups relevant to study of two-phase flow in microchannels. He developed a mechanistic model of the flow boiling phenomena based on the different forces acting on a growing vapor bubble. He introduced two new non-dimensional groups K1 and K2 relevant to flow boil 2 qL ing in microchannels. K 1 ¼ Ghq q represents the ratio between fg

G

the evaporation momentum force and the inertia force whereas  2 K 2 ¼ hq rDq represents the ratio between evaporation momenfg

G

tum force and the surface tension force. K1 is important for overall heat transfer considerations since the boiling number embedded in it is used extensively in flow boiling correlations. K2 is relevant to the interface movement during bubble growth phenomena and is more pertinent to the present study. When K2 is divided by the cosine of the contact angle it takes into account the actual surface tension force acting at the bubble base. The author identified transient heat conduction under the approaching rewetting film and thin film evaporation around the bubbles as the main mechanisms of heat transfer during flow boiling inside microchannels. He concluded that the wetting characteristic of the liquid which is indicated by the contact angle plays a significant role during dry out and rewetting and also during the events leading to CHF. Lee et al. [21] experimentally studied bubble dynamics in trapezoidal microchannels with hydraulic diameter of 41.3 lm and recorded bubble departure size and frequency using high speed digital camera. The bubble departure radius was found to decrease with heat flux whereas there was a mixed effect of mass flux on the bubble departure radius. The authors concluded that the bubble departure radius was primarily influenced by surface tension forces and drag due to bulk flow. Steinke and Kandlikar [22] experimentally studied flow boiling of water in microchannels. The most common flow boiling regime observed was annular-slug flow where a vapor bubble filled the entire channel length with a thin liquid film around it. Flow reversal in parallel microchannels occurred during explosive growth of nucleating bubbles. They recorded visual observations of the liquid vapor interface near the wall during dry-out conditions. They observed that the liquid vapor surface contact line shifted from an advancing contact position to receding contact orientation at the onset of dry out due to rapid evaporation of the liquid at the contact line region. Kandlikar et al. [23] used inlet restrictor and fabricated nucleation sites to to improve flow boiling instability in microchannels. The fabricated nucleation sites in absence of pressure drop elements increased flow boiling instability. Kosar et al. [24] used inlet orifices to suppress flow boiling instability in parallel microchannels. They defined a non-dimensional pressure drop multiplier and correlated it with the suppression of flow boiling oscillations. Kuo and Peles [25] studied the effect of reentrant cavities on suppression of flow boiling instabilities inside microchannels. Kuo and Peles [26] reported that an increase in

system pressure prevented rapid bubble growth thereby delaying flow boiling instability. Wang et al. [27] tested the effect of different inlet/outlet configurations on flow boiling instability in parallel microchannels. They found that inlet throttling led to stable flow boiling with no oscillations of temperature and pressure. Mukherjee and Kandlikar [28] numerically computed bubble growth in a microchannel with inlet constriction. The results showed presence of higher inlet velocities due to the restriction which suppressed the bubble growth and flow reversal. The authors suggested parallel microchannel designs with increasing cross-sectional area to prevent flow instability. Recently Zhang et al. [29] derived a lumped oscillator model to predict and control dynamic pressure drop instability in flow boiling microchannel systems. Several analytical and numerical models of bubble growth in microchannels have appeared in the literature. However, most of them include many simplifying assumptions that do not represent the complete physics of dynamic bubble growth. Ajaev and Homsey [30] mathematically modeled a three dimensional steady vapor bubble inside a microchannel with a square cross-section. The bottom wall of the channel was heated where evaporation took place while condensation occurred near the top wall. The inertia terms in the momentum equation and the convective heat transfer terms in the energy equation were neglected. Stress contributions from vapor recoil and thermocapillarity at the vaporliquid interface were assumed to be insignificant. The problem was solved in the limits of small capillary numbers and temperature gradients. A local solution at the contact line was obtained based on the lubrication theory. Results indicated that highest evaporative mass flux occurred in the macroscopic regions of the thin liquid film at the contact line where the disjoining pressure was not important. The overall bubble length was found to increase with the heater temperature. Mukherjee and Dhir [31] developed a three dimensional numerical model using the level-set method to study lateral merger of vapor bubbles during nucleate pool boiling. Mukherjee and Kandlikar [32] developed a complete numerical model of a growing vapor bubble during flow boiling inside a microchannel and studied the bubble growth rate for different liquid superheats and flow velocities. However, in that model the contact angle at the contact region of the liquid–vapor interface with the wall as well as the liquid and vapor thermal properties were maintained constant. The results indicated that bubble growth rate increased with increase in liquid superheat but decreased with incoming liquid flow rate. Effect of gravity was found to be negligible on bubble growth. Lee and Son [33] simulated bubble dynamics and heat transfer during nucleate boiling in a microchannel using level-set method. Later Suh et al. [34] used a similar model to simulate flow boiling in parallel microchannels. 3. Objective The various numerical and experimental studies indicate that wall heat transfer during flow boiling inside microchannels depends strongly on the wall heat flux but weakly on the mass flux. However, there is no general agreement on whether the dominant wall heat transfer mechanism is nucleate boiling or thin film evaporation. The wide difference in thermal properties and contact angle of water and dielectric liquids which are often employed for electronics cooling makes it difficult to properly ascertain their influence on the wall heat transfer. It is important to distinguish between the effects of surface wettability and thermal properties on the bubble growth and wall heat transfer in order to understand the underlying physics that governs the wall heat transfer

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

mechanisms. As a first step to that effect a vapor bubble growing inside a microchannel during flow boiling of water is simulated in the present study. All liquid and vapor properties are kept constant except the value of surface tension and contact angle which are systematically varied. The objective is to analyze and explain the effect of wall superheat, liquid mass flux, surface tension and contact angle on bubble growth and the corresponding wall heat transfer. 4. Numerical model 4.1. Computational domain Fig. 1 shows the typical computational domain. The domain is 3.96  0.99  0.99 non-dimensional units in size. Cartesian coordinates are used with uniform grid. The liquid enters the domain at x⁄ = 0 and leaves the domain at x⁄ = 3.96. To take advantage of symmetry and reduce computation time, a nucleating cavity is placed equidistant from the walls in the x–y plane. The two horizontal walls in the x–z planes are named as South Wall (y⁄ = 0) and North Wall (y⁄ = 0.99). The vertical wall in the x–y plane is named the Top Wall (z⁄ = 0.495).

3705

function. The level set function was typically a smooth function, denoted as /. This formulation eliminated the problems of adding/subtracting points to a moving grid and automatically took care of merging and breaking of the interface. Furthermore, the level set formulation generalized easily to three dimensions. The present analysis is done using this level set technique. The liquid vapor interface is identified as the zero level set of a smooth distance function /. The level set function / is negative inside the bubble and positive outside the bubble. The interface is located by solving the level set equation. A fifth order WENO (Weighted, Essentially Non-Oscillatory) scheme is used for left sided and right sided discretization of / [38]. While / is initially a distance function, it will not remain so after solving the level set equation. Maintaining / as a distance function is essential for providing the interface with a width fixed in time. This is achieved by reinitialization of /. A modification of Godunov’s method is used to determine the upwind directions. The reinitialization equation is solved in fictitious time after each fully complete time step. With Ds ¼ 2ud0 ; ten s steps are taken with a third order total variation diminishing (TVD) Runge Kutta method. 4.4. Governing equations Momentum equation:

4.2. Grid independence check The number of computational cells in the domain is 320  80  40, i.e. 80 grids are used per 0.99l0. This grid size is chosen from previous work of Mukherjee and Kandlikar [32] to minimize numerical error and optimize computation time. Variable time step is used which varied typically between 1e4 and 1e5. Negligible change in the results is observed when calculations are repeated with half the time step, which ensured that calculations are time step independent. 4.3. Method

  o~ u u ¼ rp þ q~ g  qbT ðT  T sat Þ~ þ~ u  r~ g  rjrH þ r ot

q

 lr~ u þ r  lr~ uT :

ð1Þ

Energy equation:

qC p

  oT þ~ u  rT ¼ r  krT ot

for u > 0;

ð2Þ

for u 6 0:

T ¼ T sat

Continuity equation:

The complete incompressible Navier–Stokes equations are solved using the SIMPLER method [35], which stands for Semi-Implicit Method for Pressure-Linked Equations Revised. The continuity equation is turned into an equation for the pressure correction. A pressure field is extracted from the given velocity field. At each iteration, the velocities are corrected using velocity-correction formulas. The computations proceed to convergence via a series of continuity satisfying velocity fields. The algebraic equations are solved using the line-by-line technique, which uses tri-diagonal matrix algorithm (TDMA) as the basic unit. The speed of convergence of the line-by-line technique is further increased by supplementing it with the block-correction procedure [36]. The multigrid technique is employed to solve the pressure equations. Sussman et al. [37] developed a level set approach where the interface was captured implicitly as the zero level set of a smooth

r ~ u¼

~ m

q2

 rq:

ð3Þ

The curvature of the interface:

jðuÞ ¼ r 





ru : jruj

ð4Þ

The mass flux of liquid evaporating at the interface:

~¼ m

kl rT : hfg

ð5Þ

The vapor velocity at the interface due to evaporation:

~ uevp ¼

~ m

qv

¼

Fig. 1. Computational domain.

kl rT

qv hfg

:

ð6Þ

3706

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

To prevent instabilities at the interface, the density and viscosity are defined as:





q ¼ qv þ ql  qv H;   l ¼ lv þ ll  lv H:

ð7Þ ð8Þ

100 °C (T⁄ = 0). The initial liquid temperature inside the domain is set equal to the inlet liquid temperature of 102 °C. All physical properties are taken at 100 °C. The contact angle at the walls is set as 40° unless otherwise specified, which is obtained from the experimental data of Balasubramanian and Kandlikar [39].

H is the Heaviside function given by [37]:

H ¼ 1 if u P þ1:5d; H ¼ 0 if u 6 1:5d;

4.7. Boundary conditions

ð9Þ

H ¼ 0:5 þ u=ð3dÞ þ sin ½2pu=ð3dÞ=ð2pÞ if juj 6 1:5d:

At the inlet (x⁄ = 0):

where d is the grid spacing. Since the vapor is assumed to remain at saturation temperature, the thermal conductivity is given by

k ¼ kl H1 :

ð10Þ

The level set equation is solved as:

 ou  þ ~ u þ~ uevp  ru ¼ 0: ot

ð11Þ

After every time step the level-set function /, is reinitialized as:

ou ¼ Sðu0 Þð1  jrujÞu0 ; ot uðx; 0Þ ¼ u0 ðxÞ:

u

ð13Þ

1 A

Z

ð14Þ

A

hdA;

ð15Þ

0

where A is the wall area and h is obtained from:

h¼

kl oT jwall oy

for horizontal walls; T w  T sat kl oT jwall oz for vertical walls: and h ¼ T w  T sat

ux ¼ 0:

At the outlet (x⁄ = 3.96):

ux ¼ 0:

0 hl : kl

ð19Þ

⁄

At the plane of symmetry (z = 0):

uz ¼ 0:

ð20Þ

⁄

At the walls (y = 0, y = 0.99):

u ¼ v ¼ w ¼ 0;

T ¼ Tw;

uy ¼  cos /;

ð21Þ

u ¼ v ¼ w ¼ 0;

T ¼ Tw;

uz ¼  cos /:

ð22Þ

4.8. Experimental validation Experimental studies were conducted to obtain single bubble growth data during flow boiling inside a microchannel. The data is used to validate the numerical model. The experimental flow loop consists of a syringe pump, a constant temperature bath, an in-line heater, a particle filter, and a microchannel test section (see Fig. 2). The syringe pump provides a constant flow rate of de-gassed, de-ionized water. After the syringe pump, the water is passed through a constant temperature bath consisting of a metal tube running through a beaker of continuously boiling water. The flow is then passed through an insulated in-line heater to adjust the inlet temperature of the water. A 7-lm

ð16aÞ ð16bÞ

The wall Nusselt number is defined as:

Nu ¼

ð18Þ

where u is the contact angle. At the wall (z⁄ = 0.495):

The Nusselt number (Nu) is calculated based on the area-averaged heat transfer coefficient (h) at the wall given by

h¼

T ¼ T in ;

Constant inlet flow velocity has been specified in the numerical calculations. In parallel microchannel heat exchangers constant inlet flow velocity is necessary to maintain stable operating conditions, which can be achieved using flow restrictions at the inlet, [23]:

⁄

The governing equations are made non-dimensional using a length scale and a time scale. The length scale l0 given by the channel width/height and is equal to 200 lm. Thus for water at 100 °C, and Re = 100, the velocity scale u0 is calculated as 0.146 m/s. The corresponding time scale t0 is 1.373 ms. The non-dimensional temperature is defined as:

T  T sat : T w  T sat

v ¼ w ¼ 0;

uz ¼ v z ¼ w ¼ T z ¼ 0;

4.5. Scaling factors

T ¼

u ¼ u0 ;

ux ¼ v x ¼ wx ¼ T x ¼ 0; ð12Þ

S is the sign function which is calculated as: 0 : Sðu0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u20 þ d2

The boundary conditions are as following:

ð17Þ

4.6. Initial conditions The bubble is placed at x⁄ = 0.99, y⁄ = 0 and z⁄ = 0, with 0.1l0 radius in the domain shown in Fig. 1. All velocities in the internal grid points are set to zero. The liquid inlet temperature is set to 102 °C and the wall temperature is set to the specified superheat (T⁄ = 1). The vapor inside the bubble is set to saturation temperature of

Fig. 2. Schematic of experimental microchannel flow loop.

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

Fig. 3. Comparison of bubble growth rates between experimental data and numerical simulation.

particle filter is located after the in-line heater, followed by a tee with a 0.25 mm type T thermocouple submersed in the flow to measure the fluid inlet temperature. After the microchannel test section, the water is discharged into a beaker. The microchannel test section consists of a single microchannel cut into a 400 micron-thick piece of brass using a micro-mill-

3707

ing machine. The cross-section is slightly trapezoidal, with the bottom of the trench about 40 lm smaller than the top of the trench. The nominal rectangular cross-section dimensions are 266 lm deep by 201 lm wide, giving a hydraulic diameter of 229 lm. The channel is 25.4 mm long with 7 mm of uncut brass surrounding it on all sides. A microheater is placed behind the brass with three 0.25 mm type T thermocouples located directly underneath the channel on the brass surface, between the microheater and the brass. Water is introduced into the microchannel through a polycarbonate face plate which is bolted in-place between a steel cover and a block of Teflon with a pocket milled in it. The steel cover has a through-slot cut in the center for flow visualization. The syringe pump was operated at a constant flow rate of 0.41 ml/min, which corresponds to a liquid Reynolds number of 100 using the saturation temperature of water. A fiber-optic cold light source illuminated the microchannel and high speed video was recorded at 10,000 fps over the downstream end of the channel. Temperature readings were taken every 100 ms using a laboratory data acquisition system. Bubble growth rates were measured from the high speed video images using the Phototron PFV video analysis software package. Fig. 3 shows the comparison of bubble growth rates between the simulation results and the experimental data. Numerical simulation was conducted in a square channel with 229 lm hydraulic diameter assuming saturated flow at the inlet. The wall superheat was set to the experimental value of 102.1 °C with an initial linear temperature profile in the liquid between North and South walls. The North wall boundary condition was defined as adiabatic assuming insignificant heat loss from the Plexi-glass cover used

Fig. 4. Comparison of bubble shapes: (a) numerical simulation, (b) experimental data.

3708

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

in the experiments. The wall contact angle was set to 30° as obtained from the experiments. The bubble growth rates show good agreement between the numerical simulation and experimental data. The simulation results, however, overpredicts the data with increase in time, the largest difference being around 16%. This is probably due to the assumption of constant wall temperature in the numerical simulation, whereas, local variation in wall temperatures around the bubble occurs in the experiments. There is a slight decrease in the bubble growth rate observed around 1.0 ms in both results. This is due to the fact that as the bubble expands, it touches the side walls of the microchannel leading to the formation of vapor patches that decreases heat transfer into the bubble.

Fig. 4 compares the bubble shapes between the numerical simulation and experimental data. Fig. 4(a) shows the numerically obtained bubble shapes. Fig. 4(b) shows the experimental data at similar times with the bubble outline indicated in each frame. In both figures, the bubbles are seen to move downstream along with the flow. The bubbles are seen to touch the side walls around 1.0 ms. Vapor patches can be observed on the plexi-glass surface between the bubble end caps in the experimental data at 1.4 and 1.8 ms. Similar vapor patch formation on the North wall is observed in the numerical simulation at 1.39 and 1.8 ms. Comparison of bubble growth rate and shapes indicate good agreement between the numerical results and experimental data

Fig. 5. Bubble shapes (DT = 8 K, Re = 100, r = 0.589 N/m, / ¼ 40 ).

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3709

5. Results

Fig. 6. Bubble cap locations inside the computational domain (DT = 8 K, Re = 100, r = 0.589 N/m, / ¼ 40 ).

and provides validation for the numerical model used in the following parametric study.

A parametric study has been carried out to analyze the effect of wall superheat, liquid flow rate, surface tension and contact angle on bubble growth inside a microchannel during flow boiling of water with constant incoming liquid superheat of 2 K. The wall superheat is varied as 5, 8 and 10 K. The Reynolds no. is varied as 50, 100 and 200, which correspond to inlet liquid velocities of 0.073, 0.146 and 0.29 m/s, respectively. The values of contact angle used are 20°, 40°, 60° and 80°. The surface tension values studied are 0.03, 0.04, 0.05 and 0.0589 N/m. Fig. 5 shows the growth of a vapor bubble inside the microchannel with Re = 100, wall superheat of 8 K, surface tension value of 0.0589 N/m and contact angle of 40°. The bubble grows predominantly in the direction of flow, due to heat transfer from the wall and the surrounding superheated liquid. The time in each frame is indicated at the upper right corner. The bubble is also seen to form vapor patches on the vertical walls at 0.311 ms. Fig. 6 plots the bubble cap locations from the microchannel inlet as a function of time, corresponding to Fig. 5. The upstream and downstream bubble cap locations, L1 and L2 respectively, are measured from the channel inlet as indicated in the third frame in Fig. 5. The upstream cap distance L1 stays almost constant initially and increases slightly as the bubble fills up the channel length. The

Fig. 7. Bubble growth rates: (a) effect of wall superheat, (b) effect of inlet Reynolds number, (c) effect of surface tension, (d) effect of contact angle.

3710

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

downstream cap distance L2 increases linearly at the beginning of the bubble growth but afterwards increases exponentially as the bubble turns into a vapor plug. The results indicate that the bubble grows faster with time as the spherical bubble turns into an elongated bubble. Fig. 7 plots the bubble equivalent diameters as a function of time for the different values of wall superheat, Reynolds number, surface tension and contact angle. The bubble equivalent diameter increases as the bubble grows due to heat transfer. The bubble equivalent diameter is calculated assuming a sphere of equal volume. The bubble equivalent diameter is 0.04 mm at 0 ms as the initial radius was set to 0.1l0. Fig. 7(a) shows that the bubble growth rate increases with increase in wall superheat. However, in Fig. 7(b) it can be seen that bubble growth rate decreases slightly with increase in the Reynolds number. Comparing Fig. 7(a) and (b), it can also be inferred that the effect of Reynolds number on the bubble growth is less as compared to the effect of wall superheat for the range of parameters used in the calculations. Fig. 7(c) compares the bubble equivalent diameters against time for different values of surface tension. There is little difference between the equivalent diameters for different values of surface tension indicating that the bubble growth rate is unaffected by the changes in the values of surface tension. Fig. 7(d) plots the bubble equivalent diameters against time for the cases with different con-

tact angle. The case with 20° contact angle clearly shows the highest growth rate compared to the other three values of contact angle. In this case, the presence of liquid film around and below the downstream cap of the bubble contributed to the increased wall heat transfer which will be explained further in the subsequent sections of this paper. Fig. 8 shows the effect of the wall superheat on the bubble size and shape after 0.3 ms of growth. The bubble is seen to grow comparatively slowly at 5 K wall superheat as it elongates along the channel in the direction of flow. The bubble at 10 K wall superheat is the biggest in size and has filled the entire computational domain along the length of the channel. Both the bubbles at 8 and 10 K wall superheat are seen to have developed vapor patches on the vertical walls along the x–y planes. Fig. 9 compares the bubble shapes for the two limiting values of surface tension used in this study, i.e. 0.0589 and 0.03 N/m. The frames show the bubbles when it has grown large enough to almost fill the entire channel length. The time taken for the bubbles to grow is 0.3 ms which is shown in the lower right corner of each frame. The two bubbles show very similar shapes and sizes in spite of widely varying surface tension values. One noted difference is that the bubble with higher surface tension has formed vapor patches at the vertical walls. This is because, higher surface tension causes the bubble to try to maintain its spherical shape and hence it grows comparatively more in the lateral direction. The bubble

Fig. 8. Effect of wall superheat on bubble shapes at 0.3 ms.

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3711

Fig. 9. Effect of surface tension on bubble shapes.

with lower value of surface tension does not show vapor patches at the walls in the x–y plane indicating higher wall heat transfer locally around the bubble. Fig. 10 shows the effect of contact angle on bubble shapes inside the microchannel after they grew and elongated along the channel length. Calculations were stopped just before the bubbles reached the channel outlet. The contact angle and the time corresponding to each case have been indicated at the lower right corner of each frame. The figure shows that for the case with 20° contact angle, the bubble has taken the least time to fill the channel length. Vapor patch formation on the vertical walls in the x–y plane can be seen for all the cases except for the 20° case. The extent of the vapor patch formation is found to increase with the contact angle indicating that the surface wettability decreases with increase in contact angle. The bubble base and the vapor patches at the vertical walls have merged for the case with 80° contact angle exposing large area of the walls to the vapor thereby decreasing wall heat transfer. The case with 20° contact angle shows a thin layer of liquid trapped between the wall and the bubble below the downstream end (starting from around x⁄ = 2.6). The extent of this liquid layer between the bubble base and the downstream bubble cap location can be seen to decrease with increase in contact angle. Figs. 11–13 compare the area averaged heat transfer at the North, South and Top Wall (refer Fig. 1) respectively for all the cases as a function of time. The Nu is calculated using Eqs. (16) and (17) as explained earlier. The wall heat transfer in each case is very high initially (t  0) as the incoming liquid contacts the heated wall. The wall heat transfer decreases with time as the thermal boundary begins to develop. At all the microchannel walls, it can be seen that the heat transfer improves with increase in the wall superheat. This is because the bubble grows faster with increase in the wall superheat, pushing the liquid towards the opposite wall, thereby inhibiting the thermal boundary layer development. Fig. 11 plots the effect of the different parameters on wall heat transfer at the North Wall. It is seen in Fig. 11(a) that the wall heat transfer decreases initially but can be seen to increase at the later

stages of bubble growth, as the bubble changes to vapor plug and elongates. The vapor plug elongates and in a sweeping action presses the already developed thermal boundary layer at its downstream end towards the wall, thereby causing increased wall heat transfer. Fig. 11(b) shows the effect of Reynolds number on the heat transfer at the North Wall. Compared to Fig. 11(a) it can be seen that Reynolds number, or the mass flux has little effect on the wall heat transfer at the North Wall. This is because the velocities associated with bubble growth are much higher compared to the incoming liquid velocity. From Fig. 5 it can be seen that the rate of change of L2 is initially about 1 m/s which increases to 2 m/s at later stages of bubble growth, whereas the specified incoming liquid velocity is only around 0.1 m/s. Thus, the boundary layer development and the wall heat transfer are primarily influenced by the bubble growth and not by the incoming liquid mass flux. Fig. 11(b) also shows that the heat transfer at the North Wall increases slightly after 0.3 ms with a decrease in Re. This is because the bubble growth rate increases with a decrease in Re, which opposes the expansion of the thermal boundary layer. Fig. 11(c) compares the wall averaged Nusselt numbers for all the different values of surface tension as function of time at the North wall. The heat transfer decreases initially as the flow develops and the thermal boundary layer thickens, but as the bubble grows and approaches the wall it pushes the thermal boundary layer against the wall causing the heat transfer to increase after 0.25 ms. There is little difference noted between the heat transfers at the North wall for different values of surface tension. Fig. 11(d) shows the effect of contact angle on the wall averaged Nusselt number at the North Wall as a function of time. As the thermal boundary layer thickens with time, the Nusselt number decreases until the bubble growth on the opposite South wall prevents it from thickening any further. As the bubble growing on the South wall elongates along the channel it eventually starts to push back the thermal boundary layer at the North Wall increasing the heat transfer as seen around 0.2 ms for the 20° contact angle case. This case has the highest bubble growth rate as was explained

3712

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

Fig. 10. Effect of contact angle on bubble shapes.

earlier with reference to Fig. 7. The larger bubble influences the thermal boundary layer development over a wider area and hence it results in the highest wall heat transfer compared to the other three cases. Fig. 12 shows the area averaged heat transfer at the South Wall for all the cases. The heat transfer decreases continuously from time 0 ms as the bubble base expands. At the bubble base, vapor is in contact with the wall and since the vapor has much lower thermal conductivity as compared to liquid, the wall heat transfer decreases with increase in the bubble base diameter. Fig. 12(a) shows the effect of wall superheat. At a time of 0.1 ms, the wall heat transfer is seen to increase with the increase in the

wall superheat. This is because the bubble grows faster at a higher wall superheat, which causes the liquid to move faster downstream in the channel resulting in a thinner thermal boundary layer at the South Wall. However, as the bubble base area increases with the bubble growth, the wall heat transfer decreases with time for all the cases due to the vapor contact. At around 0.3 ms, the wall heat transfer at the South Wall for the bubble at 10 K superheat has decreased considerably and is close to that of the bubble at 5 K superheat, due to its larger base area. Fig. 12(b) shows the effect of the inlet liquid flow rates. The bubble growth rate is suppressed by higher flow rates and hence the bubble base also expands at a slower rate. Thus, the higher flow

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3713

Fig. 11. Heat transfer at the north wall: (a) effect of wall superheat, (b) effect of inlet Reynolds number, (c) effect of surface tension, (d) effect of contact angle.

rate not only improves the convective heat transfer at the walls but also ensures that a larger area at the wall is in contact with the liquid. Fig. 12(c) shows the effect of surface tension on the heat transfer at the South Wall. The bubble with least surface tension 0.03 N/ m has expanded most along the channel length and has a larger base area and hence it shows least wall heat transfer at around 0.3 ms. Fig. 12(d) shows the effect of contact angle on the wall heat transfer at the South wall as a function of time. The increase in bubble base area depends on surface tension forces as well as evaporation momentum forces. Decrease in contact angle increases the surface tension forces at the bubble base resulting in decreased bubble base diameter. At 0.2 ms, the Nu at the South wall increases with decrease in contact angle due to smaller bubble base diameters except for the case with 20° contact angle. The evaporation momentum force in this case is high (due to the thin liquid layer around the bubble base) increasing the bubble base diameter and decreasing wall heat transfer at the South wall as compared to the 40° contact angle case. Fig. 13 shows the effect of the different parameters at the Top Wall (refer Fig. 1). The heat transfer is high initially (t  0) which decreases due to the thermal boundary layer development. The wall heat transfer decreases more rapidly at the later stages of bubble growth due to formation of vapor patches at the Top Wall as was seen earlier in Fig. 5.

Fig. 13(a) shows that heat transfer improves with increase in wall superheat due to increased bubble growth rates. Fig. 13(b) indicates little effect of the liquid flow rate on the wall heat transfer until the vapor patch forms at the Top Wall which is around 0.25 ms. After initiation of the vapor patch formation, heat transfer improves with Re. This is because higher liquid flow rates inhibit further expansion of the vapor patch. Fig. 13(c) shows that surface tension values too have little effect on the wall heat transfer at the Top Wall until the initiation of a vapor patch at around 0.27 ms. However, the vapor patch formation is slightly delayed by the lower surface tension values causing a small improvement in the wall heat transfer. Fig. 13(d) shows the effect of the contact angle values. The bubble with the 20° contact angle shows the highest wall heat transfer at around 0.2 ms due to its higher growth rate. However, after the initiation of the vapor patch the wall heat transfer decreases rapidly for all the cases. Fig. 14 compares the temperature field around the vapor bubbles at the central x–y plane inside the microchannel, for wall superheats of 5 and 10 K, after 0.3 ms of bubble growth. Nondimensional temperature contours are plotted between 0 and 1 with intervals of 0.1. The plot shows the thermal boundary layers at the North and South Walls. The wall heat transfer increases with wall superheat at the North Wall due to two reasons. The faster growing bubble pushes the thermal boundary layer thinner and at the same time the bigger bubble affects a larger area along the bubble length. At the South wall, larger temperature gradient is

3714

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

Fig. 12. Heat transfer at the south wall: (a) effect of wall superheat, (b) effect of inlet Reynolds number, (c) effect of surface tension, (d) effect of contact angle.

observed in the liquid at the upstream end of the bubble at 10 K superheat. At the downstream end of the bubble, higher temperature gradient is also present in the liquid below the bubble nose (x⁄  3.5) for 10 K wall superheat. However, at 10 K wall superheat, a significant portion of the South Wall is in contact with the vapor which has much lower thermal conductivity as compared to the liquid. Hence, the net heat transfer at the South Wall is found to be almost equal at 0.3 ms for the two wall superheats as was seen earlier in Fig. 12(a). Fig. 15 compares the effect of the wall superheat on the velocity field at the center of the computational domain as in Fig. 14. The reference vector has been indicated in each frame and is equal to 10 units. A significantly higher rate of evaporation is observed near the upstream and downstream interfaces of the bubble at the wall for the case with 10 K superheat. The rate of evaporation is also found to be significantly higher along the bubble interface close to the North Wall at 10 K superheat. The bubble at 5 K wall superheat is growing at a significantly slower rate which is indicated by the smaller velocity vectors at the downstream of the microchannel as compared to the case with 10 K wall superheat. Fig. 16 compares the temperature field in the central vertical x– y plane in the computational domain for the cases with 20° and 80° contact angle at the same time of 0.232 ms. Isotherms are plotted of non-dimensional temperature T⁄ from 0 to 1 with intervals of 0.1. T⁄ = 0 inside the bubble in the vapor and T⁄ = 1 at the wall.

The bubble with 20° contact angle has grown much bigger compared to the bubble with 80° contact angle. Thus, the former bubble has pushed the thermal boundary layer at the North Wall more effectively increasing the local wall heat transfer as mentioned earlier with respect to Fig. 11(d). A thin layer of liquid is clearly seen below the bubble downstream cap for the case with 20° contact angle. The crowding of the isotherms between the bubble interface and the wall indicates an area of increased wall heat transfer. Fig. 17 compares the velocity field in the central vertical x–y plane in the computational domain for the cases with 20° and 80° contact angle at the same time of 0.232 ms. The reference velocity vector equal to 20 units is indicated in each frame. High rate of evaporation near the North wall can be seen in the upper frame indicated by larger velocity vectors in the vapor near the interface inside the bubble. Large velocity vectors can also be observed near the thin liquid film at the downstream end below the bubble cap indicating intense evaporation. This explains the higher rate of bubble growth for the case with 20° contact angle. The velocity vectors near the channel inlet are much smaller compared to the channel outlet in both the frames indicating that liquid is being pushed out at a faster rate due to bubble growth. The liquid flow rate indicated by the velocity vectors at the channel outlet for the case with 20° contact angle is much greater as compared to the case with 80° contact angle, indicating faster bubble growth rate in the former case.

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3715

Fig. 13. Heat transfer at the top wall: (a) effect of wall superheat, (b) effect of inlet Reynolds number, (c) effect of surface tension, (d) effect of contact angle.

6. Discussion The above results indicate that the wall heat transfer inside the microchannel is primarily influenced by the wall superheat and the contact angle made by the moving liquid–vapor interfaces at the walls. With increase in wall superheat, bubble growth rate increases pushing the liquid more rapidly against the wall and thereby increasing the wall heat transfer. This explains the experimental findings that the heat transfer coefficient during flow boiling inside microchannels is directly dependent on the wall heat flux. However, increase in bubble growth rate at higher heat fluxes also accelerates vapor patch formation at the walls which is detrimental to wall heat transfer and may lead to early CHF. Comparison of heat transfer at the three different wall shows that the wall heat transfer decreases rapidly due to growth of bubble base and vapor patch at the South and Top walls respectively. However, at the North wall, heat transfer improvement is observed as the bubble growth prevents thickening of the thermal boundary layer at that wall. The incoming liquid mass flux suppresses bubble growth inside the microchannel thereby limiting the enhancement in the wall heat transfer due to the liquid motion generated by the bubble growth. On the other hand an increase in incoming liquid mass flux delays vapor patch formation thereby increasing the wall heat

transfer. The liquid velocities generated from the bubble growth are found to be much higher in this study as compared to the incoming liquid velocities in the microchannel. Thus as a net effect, the incoming liquid mass flux does not significantly affect the wall heat transfer, as seen in Figs. 11–13. The contact angle also has significant influence on the bubble growth and wall heat transfer inside microchannels. Contact angle not only affects the surface tension forces acting at the bubble base but also influences formation of liquid layer between the bubble and the walls. Small values of contact angle increases surface tension force parallel to the heater surface that opposes bubble base expansion. At the same time smaller contact angle leads to formation of thin liquid layer between the bubble and the walls where high rate of evaporation leads to large evaporation momentum force that assists expansion of bubble base. These opposing forces affect the bubble growth and movement of the upstream and downstream interfaces which in turn affects the wall heat transfer. The bubble with 20° contact angle is seen to have the highest growth rate and wall heat transfer as compared to all other cases presented in this study. The effect of surface tension on bubble growth and corresponding wall heat transfer is found to be insignificant for the cases studied. The net surface tension force acting on a bubble at the walls depends on the surface tension values as well as the contact angle.

3716

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

Fig. 14. Effect of wall superheat on the thermal field at 0.3 ms.

Fig. 15. Effect of wall superheat on the velocity field at 0.3 ms.

However, we have observed earlier that there is a significant influence of the contact angle on bubble growth and wall heat transfer for similar values of surface tension. This indicates that surface wettability rather than the surface tension plays a more important role in determining the wall heat transfer inside microchannels.

7. Conclusions Numerical simulation of a vapor bubble growing on a heated wall inside a microchannel during flow boiling of water has been performed. The inlet liquid flow rates, the wall superheat, the

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

3717

Fig. 16. Effect of contact angle on the thermal field at 0.232 ms.

Fig. 17. Effect of contact angle on the velocity Field at 0.232 ms.

surface tension values and the contact angle at the walls have been systematically varied keeping all other properties constant. Experiments were conducted to validate the numerical model. The bubble growth rate and shapes show good agreement between numerical and experimental results.

The following observations are made from the parametric numerical study: (1) The wall heat transfer is found to improve with increase in wall superheat and the bubble growth rates.

3718

A. Mukherjee et al. / International Journal of Heat and Mass Transfer 54 (2011) 3702–3718

(2) The wall heat transfer is found to be unaffected by the changes in the incoming liquid flow rates. (3) There is little effect of the liquid surface tension on bubble growth. The bubble shapes are affected by the liquid surface tension values with lower surface tension producing longer and thinner bubbles. The effect of surface tension on wall heat transfer is found to be negligible. (4) The movements of the upstream and downstream interfaces of the bubble as well as formation of vapor patches at the walls are found to be dependent on the contact angle at the walls. A decrease in contact angle causes formation of a liquid layer between the bubble downstream interface and the wall that has significant influence on bubble growth and wall heat transfer. The bubble with the lowest contact angle exhibited the highest growth rate and also the highest wall heat transfer. (5) The wall heat transfer increases primarily due to the motion of the evaporating liquid–vapor interface. The bubble growth is found to push the liquid against the microchannel walls, thus preventing the growth of the thermal boundary layers.

Acknowledgments The work was conducted in the Thermal Analysis and Microfluidics Laboratory at RIT and in the Advanced Energy Systems and Microfluidics Laboratory at MTU. References [1] A.E. Bergles et al., Boiling and evaporation in small diameter channels, Heat Transfer Eng. 24 (2003) 18–40. [2] S.V. Garimella, C.B. Sobhan, Transport in microchannels – a critical review, Annu. Rev. Heat Transfer 13 (2003) 1–50. [3] J.R. Thome, Boiling in microchannels: a review of experiment and theory, Int. J. Heat Fluid Flow 25 (2) (2004) 128–139. [4] S.G. Kandlikar, Fundamental issues related to flow boiling in minichannels and microchannels, Exp. Therm. Fluid Sci. 26 (2–4) (2002) 389–407. [5] T.-H. Yen et al., Visualization of convective boiling heat transfer in single microchannels with different shaped cross-sections, Int. J. Heat Mass Transfer 49 (21–22) (2006) 3884–3894. [6] B. Agostini et al., High heat flux flow boiling in silicon multi-microchannels – Part III: saturated critical heat flux of R236fa and two-phase pressure drops, Int. J. Heat Mass Transfer 51 (21–22) (2008) 5426–5442. [7] B. Agostini et al., High heat flux flow boiling in silicon multi-microchannels – Part I: heat transfer characteristics of refrigerant R236fa, Int. J. Heat Mass Transfer 51 (21–22) (2008) 5400–5414. [8] B. Agostini et al., High heat flux flow boiling in silicon multi-microchannels – Part II: heat transfer characteristics of refrigerant R245fa, Int. J. Heat Mass Transfer 51 (21–22) (2008) 5415–5425. [9] P.C. Lee, C. Pan, On the eruptive boiling in silicon-based microchannels, Int. J. Heat Mass Transfer 51 (19–20) (2008) 4841–4849. [10] S.S. Bertsch, E.A. Groll, S.V. Garimella, Refrigerant flow boiling heat transfer in parallel microchannels as a function of local vapor quality, Int. J. Heat Mass Transfer 51 (19–20) (2008) 4775–4787. [11] P.-S. Lee, S.V. Garimella, Saturated flow boiling heat transfer and pressure drop in silicon microchannel arrays, Int. J. Heat Mass Transfer 51 (3–4) (2008) 789– 806. [12] G. Wang, P. Cheng, Subcooled flow boiling and microbubble emission boiling phenomena in a partially heated microchannel, Int. J. Heat Mass Transfer 52 (1–2) (2009) 79–91.

[13] K.J.L. Geisler, A. Bar-Cohen, Confinement effects on nucleate boiling and critical heat flux in buoyancy-driven microchannels, Int. J. Heat Mass Transfer 52 (11– 12) (2009) 2427–2436. [14] T. Harirchian, S.V. Garimella, A comprehensive flow regime map for microchannel flow boiling with quantitative transition criteria, Int. J. Heat Mass Transfer 53 (13–14) (2010) 2694–2702. [15] Z.L. Yang, B. Palm, B.R. Sehgal, Numerical simulation of bubbly two-phase flow in a narrow channel, Int. J. Heat Mass Transfer 45 (3) (2002) 631–639. [16] A.M. Jacobi, J.R. Thome, Heat transfer model for evaporation of elongated bubble flows in microchannels, J. Heat Transfer 124 (6) (2002) 1131–1136. [17] S. Mukherjee, I. Mudawar, Smart pumpless loop for micro-channel electronic cooling using flat and enhanced surfaces, IEEE Trans. Compon. Pack. Technol. 26 (1) (2003) 99–109. [18] J.R. Thome, V. Dupont, A.M. Jacobi, Heat transfer model for evaporation in microchannels – Part I: presentation of the model, Int. J. Heat Mass Transfer 47 (14–16) (2004) 3375–3385. [19] V. Dupont, J.R. Thome, A.M. Jacobi, Heat transfer model for evaporation in microchannels – Part II: comparison with the database, Int. J. Heat Mass Transfer 47 (14–16) (2004) 3387–3401. [20] S.G. Kandlikar, Heat transfer mechanisms during flow boiling in microchannels, J. Heat Transfer 126 (1) (2004) 8–16. [21] P.C. Lee, F.G. Tseng, C. Pan, Bubble dynamics in microchannels – Part I: single microchannel, Int. J. Heat Mass Transfer 47 (25) (2004) 5575–5589. [22] M.E. Steinke, S.G. Kandlikar, An experimental investigation of flow boiling characteristics of water in parallel microchannels, J. Heat Transfer 126 (4) (2004) 518–526. [23] S.G. Kandlikar et al., Stabilization of flow boiling in microchannels using pressure drop elements and fabricated nucleation sites, J. Heat Transfer 128 (4) (2006) 389–396. [24] A. Kosar, C.-J. Kuo, Y. Peles, Suppression of boiling flow oscillations in parallel microchannels by inlet restrictors, J. Heat Transfer 128 (3) (2006) 251–260. [25] C.J. Kuo, Y. Peles, Flow boiling instabilities in microchannels and means for mitigation by reentrant cavities, J. Heat Transfer 130 (7) (2008) 072402– 072410. [26] C.J. Kuo, Y. Peles, Pressure effects on flow boiling instabilities in parallel microchannels, Int. J. Heat Mass Transfer 52 (1–2) (2009) 271–280. [27] G. Wang, P. Cheng, A.E. Bergles, Effects of inlet/outlet configurations on flow boiling instability in parallel microchannels, Int. J. Heat Mass Transfer 51 (9– 10) (2008) 2267–2281. [28] A. Mukherjee, S.G. Kandlikar, The effect of inlet constriction on bubble growth during flow boiling in microchannels, Int. J. Heat Mass Transfer 52 (21–22) (2009) 5204–5212. [29] T. Zhang et al., Analysis and active control of pressure-drop flow instabilities in boiling microchannel systems, Int. J. Heat Mass Transfer 53 (11–12) (2010) 2347–2360. [30] V.S. Ajaev, G.M. Homsy, Three-dimensional steady vapor bubbles in rectangular microchannels, J. Colloid Interf. Sci. 244 (1) (2001) 180–189. [31] A. Mukherjee, V.K. Dhir, Study of lateral merger of vapor bubbles during nucleate pool boiling, J. Heat Transfer 126 (6) (2004) 1023–1039. [32] A. Mukherjee, S.G. Kandlikar, Numerical simulation of growth of a vapor bubble during flow boiling of water in a microchannel, J. Microfluidics Nanofluidics 1 (2) (2005) 137–145. [33] W. Lee, G. Son, Bubble dynamics and heat transfer during nucleate boiling in a microchannel, Numer. Heat Transfer A: Appl. Int. J. Comput. Meth. 53 (10) (2008) 1074–1090. [34] Y. Suh, W. Lee, G. Son, Bubble dynamics, flow, and heat transfer during flow boiling in parallel microchannels, Numer. Heat Transfer A: Appl. Int. J. Comput. Meth. 54 (4) (2008) 390–405. [35] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, Washington, DC, 1980. [36] S.V. Patankar, A calculation procedure for two-dimensional elliptic situations, Numer. Heat Transfer 4 (4) (1981) 409–425. [37] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [38] R.P. Fedkiw, et al., A Non-Oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (The Ghost Fluid Method), Department of Mathematics, UCLA, 1999. [39] P. Balasubramanian, S.G. Kandlikar, Experimental study of flow patterns, pressure drop, and flow instabilities in parallel rectangular minichannels, Heat Transfer Eng. 26 (3) (2005) 20–27.