Numerical modelling of bubble growth in microchannel using Level Set Method

International Journal of Heat and Mass Transfer 101 (2016) 719–732

Contents lists available at ScienceDirect

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

Numerical modelling of bubble growth in microchannel using Level Set Method Gaurav Katiyar, Shyamprasad Karagadde, Sandip K. Saha ⇑, Atul Sharma Department of Mechanical Engineering, Indian Institute of Technology, Bombay, Mumbai 400076, India

a r t i c l e

i n f o

Article history: Received 13 November 2015 Received in revised form 21 January 2016 Accepted 6 May 2016

Keywords: Two phase flow Microchannel Numerical Level-Set Method

a b s t r a c t Development of more efficient thermal management systems is of prime importance not only in the context of environmental and energy concerns, but also due to ever-increasing demands of computational power. Flow boiling in microchannels holds a lot of promise and is capable of removing high heat fluxes. However, the physics behind the heat transfer and fluid flow during flow boiling at micro scales is not completely understood. Various studies have been performed to classify the flow regimes and identify the dominant mode of heat transfer in two phase flow through microchannels. In the present work, a numerical study is performed to investigate the bubble dynamics in a confined microchannel. A DGLSM (Dual-Grid Level Set Method) based numerical model is used to capture the unsteady bubble interface dynamics. The Navier–Stokes equation is being solved using Finite Volume Method (FVM) based Semi-Explicit Pressure Projection Method. The effect of parameters namely contact angle, surface tension, wall superheat, Reynolds number and system pressure on the bubble dynamics and bubble growth rates is investigated. Three distinct stages of heat transfer corresponding to the rapid reduction, stabilization and enhancement of evaluated Nusselt number are identified from the parametric investigation. The results show that the system pressure plays a vital role in controlling the bubble shape, as compared to remaining parameters. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Thermal management is increasingly becoming a bottleneck for a variety of applications such as integrated circuits, solar cells, microprocessors, and energy conversion devices. The reliability of an electronic device is defined as its ability to perform a required function under stated conditions for a given time. An electronic device fails to fulfil its intended function when its application or environmental condition exceeds its tolerable limit. Theoretically, electronic components are said to be reliable, if they are operated for a long duration of time at recommended operating temperatures, known as critical temperature. However, adverse environment and unusual operation reduces the effective operating time. It has been found that a 1 °C decrease in a component temperature may lower its failure rate by as much as 4% and 10–20 °C increase in component temperature can increase its failure rate by 100% [1]. According to a survey by the US Air Force, the percentage of temperature related failures in electronics exceeded 55% [1]. International Technology Roadmap for Semiconductors (ITRS 2005, [2]) ITRS predicted that the chip power and heat flux of the future ⇑ Corresponding author. Tel.: +91 22 2576 7392; fax: +91 22 2572 6875. E-mail address: [email protected] (S.K. Saha). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.05.029 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

generation electronic devices would reach 350 W and 200 W/cm2 in 2018. Hence, there is a tremendous need for innovative cooling technologies. The new cooling technologies should also comply with the major design goals of thermal management viz. performance, cost, physical size and reliability. In the recent years, the research on boiling and two-phase flow in microchannel has gained momentum due to its ability to meet the requirement of dissipating high heat fluxes in a relatively small space and volume. Several studies have reported the experimental observations of the flow boiling phenomena in microchannels [3–9]. Authors captured the liquid–vapour interface during dry-out conditions at the contact line region and observed flow reversal in some cases indicating a high evaporation rate. Kandlikar et al. [8] studied the effect of pressure drop and artificial nucleation sites on the stability of flow boiling and concluded that the application of pressure drop alone decreases the instabilities whereas the application of artificial nucleation sites alone increases them for the same mass flow and heat flux conditions. Kandlikar [9] proposed a possible mechanism of heat transfer during flow boiling in microchannels. It was observed that convective boiling is diminished in the microchannels and the major mode responsible for heat transfer is nucleate boiling. Two non-dimensional groups K 1 and K 2 were identified, where K 1 represents the ratio of evaporative momentum force

720

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

Nomenclature cp d Fr g h h12 Hð/Þ Ja k l lc _ m _ M ~ M; ~ ^ n; n Nu p P Pr rc r c

specific heat (J/kg K) bubble diameter (mm) Froude number gravitational acceleration (m/s2) heat transfer coefficient (W/m2 K) latent heat (J/kg) Heaviside function Jacob number thermal conductivity (W/m K) length of channel (m) characteristic length (m) interfacial mass flux (kg/m2 s) non-dimensional interfacial mass flux normal and unit vector to the interface Nusselt number pressure (N/m2) non-dimensional pressure Prandtl number Radius of bubble (m)
 non-dimensional radius of bubble ¼ rlcc

Re t T u; ~ u uc U; V x; L; l X, Y X⁄ Y⁄ Xc, Yc W W⁄ We

Reynolds number based on channel width time (s) temperature (K) velocity and velocity vector (m/s) characteristic velocity (m/s) non-dimensional velocities in X–Y coordinate length (m) coordinate axes non-dimensional length non-dimensional width non-dimensional location of bubble in X–Y coordinate width of channel (m)
 non-dimensional width of channel ¼ W lc Weber number

dð/Þ

g



c j l

r /

q r h

s u 1

Dirac delta function viscosity ratio (l2 =l1 ) diffused half interface thickness specific heat ratio (cp;2 =cp;1 ) curvature of the interface dynamic viscosity (Pa s) gradient operator level set function density (kg/m3) coefficient of surface tension contact angle non-dimensional time improved reinitialization parameter thermal conductivity ratio (k2 =k1 )

Subscripts 1; 2 liquid (fluid 1) and vapour (fluid 2) DC downstream cap cap bubble cap eq equivalent i index in inlet int interface m mean max maximum s pseudo p phase PC phase change sat saturation UC upstream cap w wall Superscripts ⁄ non-dimensional quantity ? vector quantity

Greek symbols v density ratio (q2 =q1 )

and inertia force and K 2 represents the ratio of evaporative momentum force and surface tension force. The effect of surface tension was found to be significant in microchannels, which leads to slug-annular flow regime, when the vapour bubble fills the entire channel allowing the liquid to flow around the bubble in thin films. Moreover, the primary reason for flow boiling instability was attributed to the explosive growth of bubble after its nucleation, resulting in reversal of flow in parallel microchannels. Numerical studies have also been performed by various researchers for both pool boiling and flow boiling [10–17]. Tomar et al. [12] performed a study on bubble formation in film boiling using Coupled Level Set and Volume of Fluid (CLSVOF) method on water and a refrigerant, R134a at near and far critical pressures. They concluded that an increase in wall superheat initially reduces instability and later enhances it and results in higher frequency of bubble formation. A similar numerical study of film boiling was performed by Gada and Sharma [13], where Dual Grid Level Set Method (DGLSM) was used to study the ebullition cycle of water vapour bubbles at near critical pressures for both normal and reduced gravity conditions. They reported periodic release of bubbles at nodes and antinodes alternatively. Son et al. [14] simulated nucleate boiling on a horizontal surface using Level Set Method (LSM) and studied the effect of static contact angle and wall superheat on bubble dynamics. They included the effect of microlayer

evaporation in their analysis. They reported outwards and inwards motion of the vapour contact line with the wall during the growth and departure of the bubble respectively. Dhir et al. [15] employed LSM to study the effect of wall superheat, liquid subcooling, contact angle, gravity level, noncondensables and conjugate heat transfer on the bubble dynamics of single as well as multiple bubbles during pool boiling. Pan and Chen [16] performed a twodimensional numerical study of bubble dynamics in a microchannel using a front tracking method based on immersed-boundary method. They obtained a regime map for bubble slug passing through a contraction section using Weber number (We) and Reynolds number (Re). They also studied the interactions between multiple bubbles rising through a cavity under gravity. Mukherjee et al. [17] performed a three-dimensional numerical study of the bubble growth and heat transfer where-in they studied the effects of static contact angle, Reynolds number, wall superheat and surface tension on bubble dynamics and heat transfer. However, the mechanisms of bubble growth and heat transfer are not fully understood. Besides, several aspects such as influence of pressure and channel dimensions on the onset of dry out and boiling instabilities are yet to be systematically investigated. Therefore, in this paper, a numerical study is performed on the heat transfer during flow boiling in microchannels using different parameters such as wall superheat, surface tension, contact angle,

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

Reynolds number and system pressure. Further the influence of these parameters on bubble shape, growth rate and Nusselt number is investigated, which provide critical understanding of the mechanism of heat transfer. For this purpose, a pre-defined vapour bubble located at the bottom wall of a two dimensional microchannel is analysed using Computational Multi-Fluid Dynamics (CMFD) methodology. 2. Description of the physical problem The physical problem consists of a vapour bubble located at a particular location at the bottom wall of a two dimensional channel as shown in Fig. 1. Both the bottom and top walls are assigned a constant wall temperature higher than the liquid inlet temperature. The inlet liquid temperature is kept slightly higher than the saturation temperature of the liquid, in order to sustain the bubble growth after nucleation, Hsu [18]. Constant velocity is employed at the inlet which leads to hydrodynamically and thermally developing flow. At the outlet, the insulated boundary wall condition is assigned for the temperature and outflow boundary conditions are employed for velocity and pressure. The two-dimensional computational domain (dimensionless, sized 4  1) with the nomenclature and boundary conditions is shown in Fig. 1. The boundaries are named as top wall (Y⁄ = 1), bottom wall (Y⁄ = 0), inlet (X⁄ = 0) and outlet (X⁄ = 4). The nucleating cavity is located at X⁄ = 1, Y⁄ = 0. 3. Numerical modelling and CMFD methodology A two-dimensional transient two-phase laminar flow numerical model is developed for studying the bubble dynamics in a microchannel heated from top and bottom walls. Continuity, momentum and energy equations together with the level set (LS) equations form the governing equations. The continuity, momentum and energy equations are used to calculate the temporal evolution of the flow variables such as velocity, pressure and temperature, whereas the LS equation is used for calculating the temporal evaluation of the interface. Scaling factors used for nondimensionalizing the governing equations are given as,

u x tuc p T  T sat ~¼ ~ _ U ; X ¼ ; s ¼ ; P ¼ ; T ¼ ; M uc lc lc q1 u2c T w  T sat _ q l k2 cp;2 m ; v ¼ 2 ; g ¼ 2 ;1 ¼ ; c ¼ ¼ q1 uc q1 l1 k1 cp;1 where lc and uc are characteristic length and velocity scales, respectively. The non-dimensional governing parameters such as Reynolds number, Prandtl number, Weber number, Froude number and Jacob number are defined as,

Rei ¼ Ja ¼

721

qi lc uc l cp;i uc q u2 lc ; Pri ¼ i ; Fr ¼ pffiffiffiffiffiffi ; We ¼ 1 c ; li ki r glc

cp;2 ðT w  T sat Þ h12

(a) Length scale (lc Þ is given by the channel width and is equal to 0.2 mm. (b) Velocity scale uc = 0.146 m/s based on water at Relc ¼ 100 and at 100 °C. (c) Time scale (tc ) is given by lc =uc , which results in t c ¼ 1:373 ms. Non-dimensional continuity, momentum and energy equations can be written using non-dimensional parameters from Gada and Sharma [13,19] as, 1. Continuity equation:

~¼ rU

  1v _ Md2 ð/Þ

v

ð1Þ

_ is the interfacial mass flux calculated as, where M

_ ¼ M

Ja

cRe1 Pr1

^ ½ðrT  Þ1  1ðrT  Þ2   n

ð2Þ

2. Momentum equation:

 ~ @ðqm UÞ ~UÞ ~ ¼ rP þ 1 r  ð2l DÞ  qm ^j þ r  ðqm U m @s Re1 Fr 2 1 þ jn^d2 ð/Þ We

ð3Þ

~ þ ðrUÞ ~ T Þ, ^J is the unit vector in Y-direction. qm where D ¼ 0:5 ðrU and lm are non-dimensional density and viscosity given as,

qm ¼ H ð/Þ þ v½1  H2 ð/Þ lm ¼ H ð/Þ þ g½1  H2 ð/Þ / ^ ^ ¼ jr n r/j and j ¼ rn

ð4Þ

3. Energy conservation equation:

@T  ~  Þ ¼ 1 r2 T  þ r  ðUT Rei Pr i @s

ð5Þ

where i ¼ 1 if / < 0; i ¼ 2 if / > 0 and T ¼ T sat if / ¼ 0. 4. Level-Set advection equation:

@/ ~ þ U a  r/ ¼ 0 @s

ð6Þ

~a is the sum of bulk velocity U ~ and phase change velocity where U ~PC which is defined as, U

  ~PC ¼ 0:5 1  v M _n ^ U

v

Fig. 1. Schematic of the computational domain with boundary conditions and terminology.

ð7Þ

722

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

4.1. Initial and boundary conditions

5. Reinitialization equation:

@/ ¼ Lð/0 ; /Þ þ uf ð/Þ @ ss

ð8Þ

ss is the pseudo time and Lð/0 ; /Þ ¼ Se ð/0 Þð1  jr/jÞ; f ð/Þ ¼ d2 ð/Þjr/j and reinitialization parameter is defined as, R  d2 ð/ÞLð/0 ; /Þ u¼ R d2 ð/Þf ð/Þ

ð9Þ

6. Extension equation: The normal temperature gradients in fluids are extrapolated inside the interface using the PDE described by Son and Dhir [20] as,

    @ @T  @T ^r þ Sð/Þn ¼0 @ ss @n i @n i

ð10Þ

where i ¼ 1 for extrapolation inside fluid 1 and i ¼ 2 for extrapolation inside fluid 2, ss is the pseudo time and Sð/Þ is the sharp sign function calculated as,

8 > < 1 if Sð/Þ ¼ 0 if > : 1 if

9 / > 0> = /¼0 > ; /<0

The bubble radius is taken to be 0.1 lc and placed at X  ¼ 1; Y  ¼ 0 in the computational domain as shown in Fig. 1. Velocities at all the internal grid points are set to zero initially. Liquid inlet temperature is taken to be 102 °C. The initial temperature of the entire domain is set to this temperature. Wall temperature is set as a specified degree of superheat, i.e. T  ¼ 1. The saturation temperature T sat is taken as 100 °C. All the thermophysical properties are considered at T sat . Contact angle at the walls is set at a constant value of 40° based on the experimental work of Balasubramanium and Kandlikar [21]. The boundary condition applied on the boundaries can be written as, Inlet boundary conditions: at X  ¼ 0 :

U ¼ 1; V ¼ 0; T  ¼ ðT  102Þ=ðT w  T sat Þ;

@/ @X 

¼0

Outlet boundary conditions: at X  ¼ 4 :

@U @V @T  @/ ¼ 0; ¼0  ¼  ¼ @X @X @X  @X 

ð11Þ

U ¼ V ¼ 0; T  ¼ 1;

@/ ¼ cos h @Y 

@T  Nu ¼   @y

4.2. Verification of the present numerical model

ð12Þ

The average Nusselt number is obtained as follows,

1 L

L

Nu dx

ð13Þ

0

4. Solution procedure The solution methodology for the present work is adopted from Gada and Sharma [13]. The Navier–Stokes equations along with Level-Set (LS) equations are solved on a Cartesian staggered dual grid arrangement using an in-house code. The method is called as Dual-Grid Level-Set (LS) Method (DGLSM) [13], where the Level-Set and energy equations are solved on a uniform grid which is twice finer than that of the continuity and momentum equations. This novel grid arrangement is aimed at improving interface resolution with only a slight increase in computational time [13]. In the proposed dual-resolution grid structure, the LS grid points are available at all locations where the thermo-physical properties are needed, thus avoiding the spatial averaging of properties (needed in the traditional LSM). A Finite Volume (FVM) based Semi-Explicit Pressure Projection Method is used for solution of both continuity and momentum equations. The advection terms in the predictor step for velocity are obtained by using 3rd order Lin– Lin total variation diminishing (TVD) scheme and the diffusion terms are obtained by using central difference scheme. For advection equation, 3rd order Runge–Kutta method is applied for temporal discretization and 5th order upwind WENO (weighed essentially non-oscillatory) is used for spatial discretization. For reinitialization equation, 3rd order ENO (essentially nonoscillatory) scheme is used and pseudo time step is taken one tenth of the LS node spacing. The reinitialization equation is solved till steady state is achieved.

ð16Þ

At the bottom wall: at Y  ¼ 0:

U ¼ V ¼ 0; T  ¼ 1;

Nu ¼

ð15Þ

At the top wall: at Y  ¼ 1:

Eq. (10) is the advection equation of ð@T  =@nÞp with the advection velocity as Sð/Þ. Nusselt number: Nusselt number is calculated as,

Z

ð14Þ

@/ ¼  cos h @Y 

ð17Þ

where h is the contact angle.

For verification without phase case, the physical domain consists of a bubble as dispersed phase 1 suspended in continuous phase 2, inside a channel having a length l = 2 and width W = 1. The density as well as viscosity is taken as 2 for the dispersed and 1 for the continuous phase. The radius of the bubble is rc ¼ 0:15 and its non-dimensional location is X c ¼ 0:4 and Y c ¼ 0:5. Phase 2 is entering the channel at a uniform velocity of U in ¼ 0:01. With the passage of time, the bubble advects through the channel and as a result, its shape changes due to the action of hydrodynamic and surface tension forces. For further verification, two different cases of surface tension, r ¼ 0 and r ¼ 0:1 are considered from Yap et al. [22] and the bubble shapes at different points in time are compared. A grid size of 100  50 and a time step of Dt ¼ 0:1 are used same as that used by Yap et al. [22]. It can be observed from Fig. 2 that the change from slug like shapes of the bubble at r ¼ 0 to more circular shapes at r ¼ 0:1. The figure shows a good agreement between the present model and the published literature in predicting the bubble shape for different surface tensions. For the phase change case, the verification of the present code was reported by Gada [23] for a two-dimensional Stefan problem and by Gada and Sharma [13] for a two-dimensional film boiling. In this study, further verification is presented below by comparing the present two-dimensional model with the published three numerical as well as experimental results [17]. For verification in the case of present problem, the physical domain along with the boundary conditions is shown in Fig. 1. A channel of width 0.229 mm, static contact angle of 30° and wall superheat of 2 °C is considered for simulation. For verification, the rate of bubble growth is compared with the experimental as well as three numerical results of Mukherjee et al. [17] and is shown in Fig. 3. Upto 0.5 ms, the present two-dimensional model over-predicts the bubble growth and after that it under-predicts.

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

Fig. 2. Comparison of bubble shape for surface tension (a)

723

r ¼ 0 and (b) r ¼ 0:1 with Yap et al. [22]. LSM and VOF results correspond to Yap et al. [22]

An additional comparison is performed with the analytical results presented in Mikic et al. [24]. Their analysis considers unconstrained bubble growth in a large pool of superheated liquid, nucleating on a superheated wall. In the present study, the bubble growth is influenced by the presence of walls as well as flow, which contribute to faster growth rate compared to that of pool boiling. 4.3. Grid independence study

Fig. 3. Comparison of the present numerical results with numerical and experimental results of Mukherjee et al. [17], for temporal evolution of equivalent bubble diameter.

Overall the maximum error remains about 18% with respect to the experimental data. This deviation may be due to the fact that the present simulations are two-dimensional as compared to the three simulations in the published result [17].

A grid independence study is performed by varying the number of grid points per unit length/height of the channel as 33, 55, 80 and 120, with a uniform square control volume, considering Re = 100, wall superheat of 5 K, contact angle of 40°, and surface tension of 0.0589 N/m at atmospheric pressure. The time required for the bubble to reach to the end of the channel for the above grid sizes is found to be 2.4, 1.9, 1.75, 1.65 ms respectively, which correspond to a reduction of 26.3%, 8.6% and 6.1% (for a coarser as compared to the finest grid). Thus, the grid size of 80 per unit length, resulting in a grid size of 320  80, is considered sufficient for the grid independent results and is used in all the present simulations. 5. Results and discussion 5.1. Bubble dynamics in microchannel The dynamics of bubble in microchannel is studied for a liquid inlet Re = 100, wall superheat of 5 K, contact angle of 40°, surface

724

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

tension of 0.0589 N/m and channel width of W at atmospheric pressure. Fig. 4 shows the growth of a bubble with time along with the non-dimensional velocity vectors and temperature contours. The bubble expands due to evaporation because of heat input from bottom and top walls. The bubble also advects forward due to fluid flow from the inlet. The velocity downstream to the bubble is higher as compared to the velocity upstream at any time and also, as time progresses, due to the increase in evaporative flux at the interface, the velocity increases near the interface downstream to the bubble. Upon observing the temperature contours with the passage of time, the thermal boundary layer downstream to the bubble is disrupted due to the motion of the bubble interface which results in increase in heat transfer. Also, with the increase in bubble size with time, the top interface of the bubble reaches near the top wall which compresses and disrupts the thermal boundary layer resulting in an increase in heat transfer there. 5.2. Effect of contact angle Contact angle is one of the most important parameters involved in the equilibrium of the interface. The effect of contact angle is studied on bubble dynamics and heat transfer considering static contact angles of 20°, 40°, 60° and 80°. Fig. 5(a) shows the equivalent bubble diameters with respect to time. The equivalent diameter of bubble is obtained by calculating the diameter of a circle occupying the same area of the bubble. The figure shows that the rate of bubble growth decreases with an increase in contact angle. In the figure, points marked as (a), (b), (c) and (d) on the curves represent the inflection points from where sudden growth in the vapour bubble begins. This marks the time when the bubble reaches near the top wall and due to increase in heat influx from the top wall, the bubble size starts increasing significantly. The size

of vapour bubble is the largest for 80° as indicated by the inflection points and it decreases with decreasing contact angle. The portion of channel length occupied by the vapour bubble is also highest for 80° and decreases with decrease in contact angle. Fig. 5(b) shows the non-dimensional upstream and downstream cap locations of the vapour bubble against time. The dotted lines in the figure represent the position of downstream bubble cap, whereas the solid lines represent the location of the upstream cap. The velocity of downstream cap of the bubble is greater than that of the upstream cap. This is because liquid flow from the inlet of channel opposes the velocity of the upstream bubble cap (velocity obtained by the interface due to evaporation), whereas it aids the velocity of the downstream bubble cap. Also there is a little difference in the downstream cap positions for contact angles of 40°, 60° and 80°, however large difference is observed in case of the upstream cap. Thus it can be concluded that the difference in bubble growth rates between different contact angles arises mainly due to the upstream interface movement. Fig. 5(c) shows the bubble interfaces at the inflection points marked as (a), (b), (c) and (d) in Fig. 5(a). Starting from these points, there is an increase in the slopes indicating increase in bubble growth rates. This is because of the proximity of the top interface to the top wall, due to which heat influx from the top wall increases leading to an increase in bubble growth rate. The inflection points occur at different times for different contact angles. Fig. 6 shows the bubble interface, non-dimensional velocity vectors and temperature contours throughout the domain at 0.6 ms. The bubble size is maximum for 20° and decreases with an increase in contact angle, but the length occupied by the base of the bubble increases with increasing contact angle. As a result, heat influx from the bottom wall decreases with an increase in contact angle leading to a lower growth rate. Fig. 6 also depicts

Fig. 4. Bubble shapes, non-dimensional velocity and temperature contours, for a wall superheat of 5 K, at (a) 0.1, (b) 0.4 and (c) 0.8 ms.

725

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

4

0.35 o

=20 o =40 o =60 o =80

0.3

(b)

0.25

(c)

3

(d)

(a)

2.5

Xcap

0.2

DC: DC: DC: DC: UC: UC: UC: UC:

o

=20 =40o o =60 =80o =20o o =40 =60o o =80

*

deq (mm)

3.5

2

0.15

1.5

0.1 0.05

1 0

0.5

Time (ms)

1

0.5

1.5

0.2

(a)

0.4 Time (ms)

0.6

0.8

(b)

(c)

(d) Legends:

(e) : 20 ,

: 40 ,

: 60 ,

: 80

Fig. 5. Effect of contact angle on (a) bubble growth rate, (b) non-dimensional downstream (DC) and upstream (UC) cap locations, (c) bubble interface at points (a), (b), (c) and (d) marked in (a); and temporal variation of space averaged Nusselt Number for (d) bottom and (e) top wall. Inset figures show the axial variation of Nusselt number at 0.6 ms.

the effect of contact angles on velocity field. The velocity downstream to the bubble is greater in magnitude than the velocity upstream for each contact angle. This supports the argument that the liquid flow at the inlet supports the growth of the bubble downstream, whereas hinders the upstream bubble growth. The velocity is higher in magnitude near the interface as compared to the velocity elsewhere. The velocity near the interface also exhibits a decrease in magnitude with an increase in contact angle which suggests that, at a particular time, the velocity of interface decreases with increasing contact angle, i.e., the bubble growth rate decreases with increasing contact angle. This is again due to decrease in heat influx from the bottom wall leading to a decrease in evaporation flux which results in lower velocity near the interface. The effect of contact angle on thermal field can also be seen from the figure. The movement of interface disrupts the thermal boundary layer. Thus heat transfer is high near the interface. In case of 20° contact angle, the bubble approaches close to the top wall and disrupts the thermal boundary layer. This effect increases heat transfer from the top wall.

Fig. 5(d) shows the variation in space averaged Nusselt number at the bottom wall with time. Nusselt number at the bottom wall follows a decreasing trend as the vapour contact with the bottom wall increases with time and thus heat transfer decreases. With increase in contact angle, Nusselt number at the bottom wall decreases. This is because the length of the line occupied by the vapour increases with the increase in contact angle, which in turn, decreases heat transfer from the bottom wall. Fig. 5(d) shows the variation in space averaged Nusselt number at the top wall with time. Nusselt number decreases with time up to 0.6 ms, at which point the bubble starts to compress the thermal boundary layer at the top wall, and thus heat transfer increases. Also, the bubble growth rate is highest for 20°, as a result of which heat transfer at the top wall increases earliest for 20°. By comparing the above results with three simulations reported in Mukherjee et al. [17], the trends for temporal variation of Nu are identical. However, a slight deviation in the magnitude is observed due to the difference in temperature distribution around the bubble in the two cases. For the present two-dimensional as compared to the published three

726

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

Fig. 6. Bubble shapes, non-dimensional velocity and temperature contours for (a) h = 20° (b) h = 40° (c) h = 60° and (d) h = 80° at 0.6 ms.

[17] simulations, the current trend of variation for the temporal variation of here, and also for the temporal variation of equivalent bubble diameter in Fig. 3, indicate that the detailed parametric investigation in the present two-dimensional study probably captures the correct trend of variation for the more realistic three conditions. Three distinct stages of heat transfer from the top and bottom walls can be observed from the transient variation of Nusselt number (Fig. 5(e)), viz. (i) initial transient reduction, (ii) intermediate stabilization and (iii) final enhancement. During the initial phase, transient nature of flow in the microchannel causes reduction in the heat transfer at the top and bottom walls (Fig. 5(e)). The Nusselt number at the top wall in this phase also decreases rapidly as the thermal boundary layer is in developing stages. However, slight difference in Nu between the top and bottom walls exists due to the contact of the bubble with the bottom wall. For 20° contact angle, this initial transient phase can be observed up to 0.4 ms. In the next stage, i.e. intermediate stabilization, the overall heat transfer at top and bottom walls tends to stabilize, as the flow becomes fully developed. This phase exhibits slight drop in Nusselt number at the bottom and top walls as the gradual bubble growth alters the contact area and thermal boundary layers at the respective walls. At the top wall, the thermal boundary layer is compressed due to the growth of the bubble, however it is not disrupted during the stabilization stage. During the final stage, the bubble growth takes places very close to the top wall, and the thermal boundary layers are disrupted. This results in an increase in heat transfer along the top wall. In this phase, the contribution of heat transfer from the top wall is much higher as compared to that of the initial phase due to the proximity of the bubble to the top wall. During this phase, an increase in Nusselt number at the top wall is observed only for 20° after 0.6 ms as the bubble reaches close to the top wall (evident from Fig. 6). These three stages of heat transfer mechanism during bubble growth in microchannel are observed for all cases presented in the following sections. The inset in Fig. 5(d) shows the variation of Nusselt number with non-dimensional length along the channel ðX  Þ at the bottom

wall at 0.6 ms. The spikes in Nusselt number represent the position of the two interfaces at the bottom wall. For a given contact angle, between the two interfaces, Nusselt number decreases due to vapour contact at the bottom wall. Comparing the spike for Nusselt numbers at the downstream interface, the Nusselt number at the downstream is found to increase with the decrease in contact number due to the presence of a thin film of liquid which decreases with increasing contact angle. The variation of Nusselt number with non-dimensional length along the channel ðX  Þ at the top wall at 0.6 ms is shown in the inset of Fig. 5(e). The vapour bubble for contact angle of 20° reaches close to the top wall at around 0.6 ms and Nusselt number increases for a region where the thermal boundary layer is suppressed by the vapour bubble. The vapour bubble for 40° contact angle also reaches close to the top wall, but not as close as the vapour bubble as compared to contact angle of 20°. There is a very little difference between Nusselt numbers for 60° and 80°, because the vapour bubbles are too far away to affect the Nusselt number at the top wall. 5.3. Effect of surface tension Surface tension is one of the most important forces affecting the bubble dynamics in micro-domain. Since gravitational force is negligible in capillary scale, the interface dynamics is determined by the equilibrium between three forces, viz. evaporative flux force (the force due to the evaporative flux at the interface), the surface tension force and the inertia force due to the flow. Although the equilibrium between these forces was investigated in detail by Kandlikar [9], the presence of surface tension was not isolated. This is an important aspect as the magnitude of surface tension force along with the contact angle determines the equilibrium of interface. However, in the present study, the effect of surface tension on heat transfer and bubble dynamics is studied using three different values of surface tension which physically correspond to three different concentrations of surfactants: 0.03, 0.04 and 0.0589 N/m. Fig. 7 shows the effect of surface tension on bubble shapes, nondimensional velocity vectors and temperature contours at 1 ms.

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

727

remains unaffected by the change in surface tension up to 0.8 ms around which the top interface of the vapour bubble comes very close to the top wall. After 0.8 ms (final enhancement stage), Nusselt number at the top wall increases with an increase in surface tension, as in case of higher surface tension, the bubble growth will be more uniform resulting in compression of the thermal boundary layer at the top wall. The inset shows the spatial variation of Nusselt number at the top wall at 1 ms. The peak of Nusselt number is highest for 0.0589 N/m followed by 0.04 and 0.03 N/m, respectively. 5.4. Effect of wall superheat

Fig. 7. Bubble shapes, non-dimensional velocity and temperature contours at 1 ms, for (a) r = 0.03 N/m, (b) r = 0.04 N/m and (c) r = 0.0589 N/m.

The surface tension has negligible effect on bubble growth rates, but has an appreciable effect on the shape of the bubble. The longitudinal elongation of the bubble is maximum for a lower value of surface tension as seen in Fig. 7(a) and it decreases with an increase in the value of surface tension. This is expected as surface tension tries to decrease the surface area of the interface. There is a marginal difference between the downstream bubble cap locations for different values of surface tension. The vapour contact with the bottom wall remains almost same with different surface tension values. In accordance with prior observations (for different values of contact angle), for a particular value of surface tension, the velocity of downstream cap is higher as compared to the upstream cap. Moreover, vapour slug can be clearly observed in the microchannel. The figure also shows the temperature contours for different values of surface tension at 1 ms. It is clearly evident from the figure that the compression of thermal boundary layer at the top wall is highest in case of 0.0589 N/m. Fig. 8(a) shows the variation of space averaged Nusselt number at the bottom wall with time. Since there is no significant difference between the vapour contact lengths for different surface tension values, the heat transfer at bottom wall is almost unaffected by the change in surface tension. The inset shows the spatial variation of Nusselt number at the bottom wall at 1 ms. It is evident from the figure that there is a slight difference in Nusselt number at the bottom wall near the downstream interface. Comparing Nusselt number among different surface tensions, its value for 0.0589 N/m is found to be the highest followed by 0.04 and 0.03 N/m, respectively. Fig. 8(b) shows the variation of space averaged Nusselt number with time at the top wall. Nusselt number

The effect of wall superheat on Nusselt number, bubble shapes and velocity field is presented in this section. Wall superheat is an important parameter which affects the quantity of heat input from the wall which in turn affects the critical heat flux (CHF) and boiling instabilities. Fig. 9(a) shows the effect of wall superheat on bubble growth rates. It is clearly evident from the figure that the bubble growth rate increases with an increase in wall superheat. This is because, with the increase in wall superheat, more heat transfers through the bottom wall, which results in larger quantity of vapour causing an increase in equivalent diameter and growth rate. Fig. 9(b) depicts the effect of wall superheat on nondimensionalized bubble cap locations. With increasing wall superheat, the positions of both the upstream and downstream bubble caps move downstream at a given time. Also the variation in upstream bubble cap locations with time is insignificant as compared to that of the downstream bubble cap locations. Fig. 10 shows the bubble interface, non-dimensional velocity vectors and thermal field at 0.6 ms for different values of wall superheat. With an increase in wall superheat, the bubble size as well as the bubble elongation along the horizontal direction downstream increases. The vapour contact length also increases with an increase in wall superheat. The top interface of the vapour bubble in case of both 8 and 10 K wall superheat reaches near the top wall. For a particular wall superheat, the magnitude of velocity downstream is higher as compared with the magnitude of velocity upstream. The temperature inside the bulk of channel is higher for 5 K wall superheat as compared to 8 and 10 K. Fig. 9(c) presents the variation in space averaged Nusselt number with time at the bottom wall. It increases with an increase in wall superheat, however its value remains steady until 0.8 ms, which may be due to the increased vapour contact at the bottom wall. The inset shows the spatial variation in Nusselt number on the bottom wall at 0.6 ms. Comparing Nusselt number near the upstream interface, it can be seen from the figure that it is highest for 5 K followed by 8 and 10 K. There is no significant difference in Nu for 8 and 10 K wall superheat. This may be due to formation of thin film of liquid near upstream interface at 5 K, due to which heat transfer increases. Near the downstream interface, the maximum Nusselt number is obtained for 8 and 10 K together followed by 5 K, due to the presence of a thin film of liquid near the downstream interface. The variation in space averaged Nusselt number with time at the top wall is shown in Fig. 9(d). Nusselt number starts to increase near 0.7 ms (final enhancement stage) for both 8 and 10 K wall superheat. This is because the top interface of the bubble reaches very close to the wall which results in disruption of the thermal boundary layer at the top wall increasing heat transfer. For 5 K, there is not appreciable increase in Nusselt number. The inset shows the spatial variation in Nusselt number over the top wall at 0.6 ms. There is an increase in Nusselt number over the region where the thermal boundary layer is disrupted by the top interface of the bubble. Furthermore, the increase in heat transfer occurs over a larger area for higher values of wall superheat.

728

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

(a) Legends:

: 0.03 N/m,

(b) : 0.0589 N/m

: 0.04 N/m,

Fig. 8. Temporal variation in space averaged Nusselt number for (a) bottom and (b) top wall, at different values of surface tension coefficient. Inset figures show the axial variation of Nusselt number at 0.6 ms.

0.35

4

5K 8K 10K

0.3

3.5 3 2.5

Xcap

0.2

*

deq (mm)

0.25

DC:5K DC:8K DC:10K UC:5K UC:8K UC:10K

2

0.15

1.5 0.1 0.05

1 0

0.2

0.4 Time (ms)

0.6

0.8

0.5

(a)

(c) Legends:

0.2

0.4 Time (ms)

0.6

0.8

(b)

: 5 K,

: 8 K,

(d) : 10 K

Fig. 9. Effect of wall superheat on temporal variation of bubble (a) equivalent diameter, (b) non-dimensional downstream (DC) and upstream (UC) cap locations; and space averaged Nusselt Number for (c) bottom and (d) top wall. Inset figures show the axial variation of Nusselt number at 0.6 ms.

5.5. Effect of Reynolds number Fig. 11(a) shows the effect of Reynolds number on the bubble growth rate. With the increase in Reynolds number, the vapour bubble will reach the end of the channel faster and thus the duration of thermal contact with the wall decreases. There is a negligi-

ble difference between the equivalent bubble diameter for Re = 50 and 100, however, the bubble growth rate decreases with an increase in Reynolds number comparing the trends of Re = 50 and 200. Fig. 11(b) shows the effect of Reynolds number on the nondimensional bubble cap location with respect to time. It is evident

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

729

5.6. Effect of system pressure

Fig. 10. Bubble shapes, non-dimensional velocity and temperature contours at 0.6 ms, for wall superheat of (a) 5 K (b) 8 K and (c) 10 K.

from the figure that the velocity of downstream bubble cap is higher than that of upstream bubble cap for a given Reynolds number. With the increase in Reynolds number, both the upstream and downstream bubble caps move further downstream, which indicates the downstream motion of the bubble as a whole. Fig. 12 compares the bubble shapes, velocity and thermal fields for different Reynolds number at a particular instant of time (0.6 ms). The magnitude of velocity downstream is highest (0.63 m/s) in case of Re = 200 followed by Re = 100 (0.53 m/s) and 50 (0.38 m/s), respectively. Fig. 11(c) shows the space averaged Nusselt number as a function of time for different Reynolds number over the bottom wall, where marginal changes in Nu is observed. The inset shows the spatial variation of Nusselt number at 0.6 ms on the bottom wall for different Reynolds number. It can be noted that even though the spike in Nusselt number is maximum for Re = 50, the space averaged Nusselt number for the same Re is slightly lower in magnitude as shown in Fig. 11(c). Near the downstream and upstream locations at the bottom wall, Nusselt number is the highest for Re = 50 followed by 200 and 100, respectively. Fig. 11(d) shows the space averaged Nusselt number as a function of time for different Reynolds numbers on the top wall. For a given Reynolds number, Nusselt number at the top wall increases after the top interface of the vapour bubble reaches close to the top wall (final enhancement stage). The inset shows the spatial variation in Nusselt number at 0.6 ms over the top wall. The rise in Nusselt number indicates the disruption of thermal boundary layer by the top interface of the vapour bubble.

System pressure plays a very crucial role in determining the dynamics of fluid flow and heat transfer. Changing system pressure affects various thermophysical properties such as density, viscosity, thermal conductivity. Although the effect of these properties can be studied separately, it is important to know the way system responds if pressure is varied. Typically in microchannel applications having a range of fluids, the system pressures can be as high as 5 bars, therefore pressures of 101.3, 198.53, 361.3, 475.8 kPa are considered. These values of pressure correspond to the saturation temperatures of 100, 120, 140 and 160 °C respectively. The corresponding changes in boiling temperature of water are accounted with an equivalent increase in wall superheat and initial liquid temperatures. Fig. 13(a) shows the effect of system pressure on the bubble growth rate. It can be observed from the figure that the bubble growth rate decreases with an increase in system pressure. This is because pressure outside the vapour bubble opposes the growth of the vapour bubble due to evaporation. Fig. 13(b) depicts the non-dimensional bubble cap locations with respect to time. With the increase in system pressure, both the upstream and downstream bubble cap locations move upstream with respect to the microchannel. This is due to the decrease in size of the vapour bubble with an increase in system pressure. Fig. 14 compares the bubble shapes, velocity and thermal fields for different system pressures. The vapour contact length with the bottom wall decreases with increasing system pressure. As the system pressure increases, the magnitude of velocity downstream decreases. This can be attributed to the decrease in evaporative component of the velocity due to the rise in pressure. Thermal gradients near the downstream interface at the bottom wall are highest for 475.8 kPa due to smaller thermal boundary layer thickness. However, the thickness of thermal boundary layer does not vary significantly for lower pressures. Fig. 13(c) shows the space averaged Nusselt number as a function of time over bottom wall. With an increase in system pressure, the Nusselt number increases at the bottom wall due to decrease in vapour contact length with the bottom wall. The inset shows the spatial variation of Nusselt number at 0.6 ms over bottom wall. The peak of Nusselt number near the downstream interface at the bottom wall is highest for 475.8 kPa denoting the presence of high thermal gradients. Near the upstream interface at the bottom wall, the spike in Nusselt number is the lowest for 475.8 kPa. Fig. 13(d) shows the space averaged Nusselt number as a function of time over the top wall. The top interface of the vapour bubble reaches close to the top wall only in case of 101.33 kPa at 0.7 ms, after which heat transfer improves. Comparing Nusselt number among different system pressures, it can be concluded that Nusselt number at the top wall decreases with an increase in system pressure. The inset shows the spatial variation in Nusselt number at 0.6 ms over top wall. As mentioned earlier, only vapour bubble at 101.33 kPa is able to affect the Nusselt number at the top wall, which remains unaffected for other values of system pressures.

6. Conclusions Numerical study is performed with a bubble placed at the bottom of a two dimensional microchannel with heated top and bottom walls. Various parameters, namely static contact angle, surface tension, wall superheat, Reynolds number and system pressure are systematically varied and their effects on bubble size, bubble shape and Nusselt number on top and bottom walls are observed. It is found that for all parameters, velocity magnitude downstream is higher as compared to that of the upstream because the phase

730

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

4

0.35 0.3

3 2.5

Xcap

0.2

*

deq (mm)

0.25

DC:Re=50 DC:Re=100 DC:Re=200 UC:Re=50 UC:Re=100 UC:Re=200

3.5

Re=50 Re=100 Re=200

2

0.15

1.5 0.1

1

0.05 0

0.2

0.4 0.6 Time (ms)

0.8

1

0.5

0

(a)

(c) Legends:

0.2

0.4 0.6 Time (ms)

0.8

1

(b)

: Re = 50,

: Re = 100,

(d) : Re =200

Fig. 11. Effect of Reynolds number on temporal variation of bubble (a) equivalent diameter, (b) non-dimensional downstream (DC) and upstream (UC) cap locations; and space averaged Nusselt Number for (c) bottom and (d) top wall. Inset figures show the axial variation of Nusselt number at 0.6 ms.

Fig. 12. Bubble shapes, non-dimensional velocity and temperature contours at 0.6 ms, for (a) Re = 50 (b) Re = 100 and (c) Re = 200.

731

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

8

0.4 101.33 kPa 198.53 kPa 361.3 kPa 475.8 kPa

0.35 0.3

6 deq (mm)

deq (mm)

0.25 0.2

0.15

5 4 3

0.1

2

0.05 0

DC:101.33 kPa DC:198.53 kPa DC:361.3 kPa DC:475.8 kPa UC:101.33 kPa UC:198.53 kPa UC:361.3 kPa UC: 475.8 kPa

7

1 0

0.4

0.8

1.2 1.6 2 Time (ms)

2.4

2.8

3.2

(a)

Legends:

(c) : 101.33 kPa,

0

0.5

1 Time (ms)

1.5

2

(b)

: 198.53 kPa,

(d) : 361.3 kPa,

: 475.8 kPa

Fig. 13. Effect of system pressure on temporal variation of bubble (a) equivalent diameter and (b) non-dimensional downstream (DC) and upstream (UC) cap locations; and space averaged Nusselt Number for (c) bottom and (d) top wall. Inset figures show the axial variation of Nusselt number at 0.6 ms.

Fig. 14. Bubble shapes, non-dimensional velocity and temperature contours at 0.6 ms, for (a) 101.33 kPa (b) 198.53 kPa (c) 361.3 kPa and (d) 475.8 kPa.

732

G. Katiyar et al. / International Journal of Heat and Mass Transfer 101 (2016) 719–732

change component of the velocity is aided by the liquid flow near the downstream. By evaluating the temporal variation in Nusselt number in the microchannel, three distinct stages of heat transfer mechanism are identified in almost all parametric cases, namely, initial transient, intermediate stabilization and final enhancement. With the increase in contact angle, the bubble growth rate decreases. Nusselt number at the bottom wall is highest for contact angle of 20° and decreases with an increase in contact angle. It is found that there is a little effect of surface tension on the bubble growth rates, however the bubble shapes are influenced greatly by it. An increase in wall superheat causes the bubble growth rate to increase due to increase in the heat input from the walls. At the top wall, an early rise in Nusselt number can be observed for both 8 and 10 K as compared to 5 K due to very high bubble growth rates. More importantly, influence of system pressure on heat transfer characteristics and bubble growth is established. The bubble growth rates are found to decrease proportionately with the increasing system pressure. The numerical studies reveal important aspects of microchannel heat transfer and demonstrate their potential in performing further simulations and developing optimized microchannel designs. Acknowledgement The help and support received from Mr. Vinesh Gada, for the present in-house code, is gratefully acknowledged. References [1] E.M. Alawadhi, C.H. Amon, Thermal analysis of a PCM thermal control unit for portable electronics devices: experimental and numerical analysis, in: Eight Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems, San Diego, California, 2002. [2] http://www.itrs.net. [3] X.F. Peng, B.X. Wang, Forced convection and flow boiling heat transfer for liquid flowing through microchannels, Int. J. Heat Mass Transfer 36 (14) (1993) 3421–3427. [4] W. Qu, I. Mudawar, Flow boiling heat transfer in two-phase micro-channel heat sinks. I. Experimental investigation and assessment of correlation methods, Int. J. Heat Mass Transfer 46 (15) (2003) 2755–2771. [5] T. Hririchian, S.V. Garimella, Effects of channel dimensions, heat flux and mass flux on flow boiling regime in microchannels, Int. J. Multiph. Flow 35 (4) (2009) 349–362.

[6] S.S. Bertsch, E.A. Groll, S.V. Garimella, Effects of heat flux, mass flux, vapour quality and saturation temperature on flow boiling heat transfer in microchannels, Int. J. Multiph. Flow 35 (2) (2009) 142–154. [7] M.E. Steinke, S.G. Kandlikar, An experimental investigation of flow boiling characteristics of water in parallel microchannels, J. Heat Transfer ASME 126 (4) (2004) 518–526. [8] S.G. Kandlikar, W.K. Kuan, D.A. Willistein, J. Borrelli, Stabilization of flow boiling in microchannels using pressure drop elements and fabricated nucleation sites, J. Heat Transfer ASME 128 (4) (2006) 389–396. [9] S.G. Kandlikar, Heat transfer mechanisms during flow boiling in microchannels, J. Heat Transfer ASME 126 (2004) 8–16. [10] M. Magnini, B. Pulvirenti, J.R. Thome, Numerical investigation of hydrodynamics and heat transfer of elongated bubbles during flow boiling in a microchannel, Int. J. Heat Mass Transfer 59 (2013) 451–471. [11] Z.J. Edel, A. Mukherjee, Experimental investigation of vapour bubble growth during flow boiling in a microchannel, Int. J. Multiph. Flow 37 (10) (2011) 1257–1265. [12] G. Tomar, G. Biswas, A. Sharma, A. Agrawal, Numerical simulation of bubble growth in film boiling using coupled level-set and volume-of-fluid method, Phys. Fluids 17 (2005). 112103-1–112103-13. [13] V.H. Gada, A. Sharma, On a novel dual-grid level-set method for two-phase flow simulation, Numer. Heat Transfer Part B 59 (1) (2011) 26–57. [14] G. Son, V.K. Dhir, N. Ramanujapu, Dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface, J. Heat Transfer ASME 121 (3) (1999) 623–631. [15] V.K. Dhir, G.R. Warrier, E. Atkinol, Numerical simulation of pool boiling: a review, J. Heat Transfer ASME 135 (6) (2013). 061502-1–061502-17. [16] K.L. Pan, Z.J. Chen, Simulation of bubble dynamics in microchannel using a front-tracking method, J. Comput. Math. Appl. 67 (2) (2014) 290–306. [17] A. Mukherjee, S.G. Kandlikar, Z.J. Edel, Numerical study of bubble growth and wall heat transfer during flow boiling in a micro channel, Int. J. Heat Mass Transfer 54 (2011) 3702–3718. [18] Y.Y. Hsu, On the size range of active nucleation cavities on a heating surface, J. Heat Transfer 84 (1962) 207–216. [19] V.H. Gada, A. Sharma, On derivation and physical interpretation of level set method-based equations for two-phase flow simulations, Numer. Heat Transfer Part B 56 (4) (2009) 307–322. [20] G. Son, V.K. Dhir, A level set method for analysis of film boiling on an Immersed Solid Surface, Numer. Heat Transfer Part B 52 (2) (2007) 153–177. [21] 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. [22] Y.F. Yap, J.C. Chai, T.N. Wong, K.C. Toh, H.Y. Zhang, A global mass correction scheme for the level-set method, Numer. Heat Transfer Part B 50 (5) (2006) 455–472. [23] V.H. Gada, A novel level-set based CMFD methodology in 2D/3D Cartesian and cylindrical coordinates and its application for analysis of stratified flow and film boiling (Ph.D. thesis), Indian Institute of Technology, Bombay, 2012. [24] B.B. Mikic, W.M. Rohsenow, P. Griffith, On bubble growth rates, Int. J. Heat Mass Transfer 13 (1970) 657–666.