Modeling of pressure drop and heat transfer for flow boiling in a   mini/micro-channel of rectangular cross-section

Shashwat Jain; Prasanna Jayaramu; Sateesh Gedupudi

arXiv:1902.00264·physics.flu-dyn·February 4, 2019

Modeling of pressure drop and heat transfer for flow boiling in a mini/micro-channel of rectangular cross-section

Shashwat Jain, Prasanna Jayaramu, Sateesh Gedupudi

PDF

TL;DR

This study presents a 1D model for flow boiling in rectangular microchannels, accounting for pressure fluctuations, bubble dynamics, and dryout effects, validated against experimental data.

Contribution

It introduces a comprehensive model that includes pressure fluctuations and bubble growth effects on heat transfer in microchannels, enhancing existing predictive capabilities.

Findings

01

Model accurately predicts pressure drop and heat transfer coefficients.

02

Pressure fluctuations significantly influence heat transfer characteristics.

03

Flow reversal effects are also characterized.

Abstract

In the present study, a 1-dimensional model is proposed to estimate the pressure drop and heat transfer coefficient for flow boiling in a rectangular microchannel. The present work takes into account the pressure fluctuations caused due to the confined bubble growth and the effect of pressure fluctuations on the heat transfer characteristics. The heat transfer model considers five zones, namely, liquid slug, partially confined bubble, fully confined (elongated) bubble, partial dryout and full dry-out. The model incorporates the thinning of liquid film due to shear stress at liquid-vapour interface in addition to evaporation. The transient fluctuations in pressure and heat transfer coefficient, along with the time-averaged ones, are verified with the experimental data available in the literature. Heat transfer characteristics with flow reversal caused by inlet compressibility are also…

Figures40

Click any figure to enlarge with its caption.

Tables1

Table 1. Table 1: Various correlations for heat transfer coefficient

Reference	Correlation	Comments
Lazarek and Black [60]	$H T C_{t p} = \frac{k_{l}}{D_{h}} (30 R e_{l o}^{0.857} B o^{0.714})$	Mini channel correlations based on 728 points,
		$D_{h} = 3.1 m m$
		$G = 125 - 750 k g / m^{2} s, B o = 2.3 - 76 \times 10^{- 4}$
		$P = 1.3 - 4.1 b a r$
Kew and Cornwell [61]	$H T C_{t p} = \frac{k_{l}}{D_{h}} (30 R e_{l o}^{0.857} B o^{0.714}) {(1 - x)}^{- 0.143}$	Mini/micro channel correlation
		Based on circular data $D_{h} = 1.39 - 3.69 m m$ ,
		working fluid R141b
Sun and Mishima [55]	$H T C_{t p} = \frac{6 R e_{l o}^{1.05} B o^{0.54}}{W e_{l}^{0.191} {(ρ_{l} / ρ_{v})}^{0.142}} \frac{k_{l}}{D_{h}}$	Mini/micro channel correlation
		Based on 11 different fluid and 2505 data points
		$D_{h} = 0.21 - 6.05 m m$
		$G = 44 - 1500 k g / m^{2} s, q = 5 - 109 k W / m^{2}$
Li and Wu [56]	$H T C_{t p} = 334 B o^{0.3} {(B d R e_{l}^{0.36})}^{0.4} \frac{k_{l}}{D_{h}}$	Mini-micorchannel correlation
		Based on 3700 data points, 13 different fluids
		$D_{h} = 0.2 - 3 m m$
Mahmoud and Karayiannis [62]	For $D_{h} = 0.52 m m$ , $x \leq 0.3$	Mini/micro channel correlation
	$H T C_{t p} = 3320 \frac{B o^{0.63} W e_{l}^{0.2} R e_{l}^{0.11}}{C o^{0.6}} \frac{k_{l}}{D_{h}}$	Based on 8561 data points and R134a - working fluid
	For $D_{h} = 0.52 m m$ , $x > 0.3$	$D_{h} = 0.52 - 4.26 m m$ , $G = 100 - 500 k g / m^{2} s$
	$H T C_{t p} = 5324 {[\frac{B o^{0.3} W e_{l}^{0.25}}{N_{c o}^{0.25}}]}^{1.79} \frac{k_{l}}{D_{h}}$	$q = 2.4 - 175.4 k W / m^{2}$ , $P = 6 - 14 b a r$
Mahmoud and Karayiannis [57]	$H T C_{t p} = S_{n e w} H T C_{C o o p e r} + F_{n e w} H T C_{l}$	Mini/micro channel correlation
		Based on 5152 data points and R134a-working fluid
	$H T C_{C o o p e r} = 55 P_{r}^{0.12 - 0.434 l n (R_{p})} {(- l o g P_{r})}^{- 0.55} M_{w}^{- 0.5} q^{0.67}$	$D_{h} = 0.52 - 4.26 m m$ , $G = 100 - 700 k g / m^{2} s$
	$H T C_{l} = 4.36 k_{l} / D_{h}$ for $R e_{l} < 2000$	$q = 1.7 - 158 k W / m^{2}$ , $P = 6 - 14 b a r$
	$H T C_{l} = 0.023 R e_{l}^{0.8} P r_{l}^{0.4} \frac{k_{l}}{D_{h}}$ for $R e_{l} > 3000$
	$R e_{l} = \frac{(1 - x) G D_{h}}{μ_{l}}$
	$F_{n e w} = {(1 + \frac{2.812 C o^{- 0.408}}{X})}^{0.64}$
	$S = \frac{1}{1 + 2.56 \times 10^{- 6} {(R e_{l} F_{n e w}^{1.25})}^{1.17}}$
	$X = {(\frac{f_{l}}{f_{v}})}^{0.5} {(\frac{ρ_{v}}{ρ_{l}})}^{0.5} \frac{(1 - x)}{x}$
Lee and Mudawar [63]	For $0 < x \leq 0.05$ , $H T C_{t p} = 3.856 X^{0.267} H T C_{l}$	Microchannel correlation
	$H T C_{l} = N u_{3, \infty} k / D_{h}$	Based on water - 207 data points,
	$X_{v v} = {(\frac{μ_{l}}{μ_{v}})}^{0.5} {(\frac{1 - x}{x})}^{0.5} {(\frac{ρ_{v}}{ρ_{l}})}^{0.5}$	R134a - 111 data points
	$X_{v t} = {(\frac{f_{l} R e_{v}^{0.2}}{0.079})}^{0.5} (\frac{1 - x}{x}) {(\frac{ρ_{v}}{ρ_{l}})}^{0.5}$	$D_{h} = 349 μ m$
	$R e_{v} = G x D_{h} / μ_{v}$
	For $0.05 < x \leq 0.55$ , $H T C_{l} = 436.48 B o^{0.522} W e_{l}^{0.351} X^{0.665} H T C_{l}$
	For $0.55 < x \leq 1, H T C_{l} = M A X (108.6 X^{1.665} H T C_{v}, H T C_{v})$
	$H T C_{v} = N u_{v} k_{v} / D_{h}$
	$N u_{v} = 0.023 R e_{v}^{0.8} P r_{v}^{0.4}$

Equations122

T_{w a l l} - T_{s a t} = \frac{2 ( 1 + cos ( θ )) σ q}{ρ _{v} k _{l} i _{l v}}

T_{w a l l} - T_{s a t} = \frac{2 ( 1 + cos ( θ )) σ q}{ρ _{v} k _{l} i _{l v}}

R = \frac{12 k _{l} t _{g r o w}}{C 1 ρ _{l} C _{p, l} π} J a

R = \frac{12 k _{l} t _{g r o w}}{C 1 ρ _{l} C _{p, l} π} J a

J a = \frac{C _{p, l} Δ T _{s u p} ρ _{l}}{ρ _{v} i _{l v}}

J a = \frac{C _{p, l} Δ T _{s u p} ρ _{l}}{ρ _{v} i _{l v}}

q A_{h e a t} = \frac{d}{d t} (ρ_{v} A_{ba se} h i_{l v})

q A_{h e a t} = \frac{d}{d t} (ρ_{v} A_{ba se} h i_{l v})

τ = \frac{w h i _{l v} ρ _{v}}{p _{h} q}

τ = \frac{w h i _{l v} ρ _{v}}{p _{h} q}

\frac{d L _{b u b}}{d t} = \frac{L _{b u b}}{2 τ}

\frac{d L _{b u b}}{d t} = \frac{L _{b u b}}{2 τ}

(L_{b u b})_{t + Δ t} = (L_{b u b})_{t} + (\frac{L _{b u b}}{2 τ})_{t} Δ t

(L_{b u b})_{t + Δ t} = (L_{b u b})_{t} + (\frac{L _{b u b}}{2 τ})_{t} Δ t

U_{d} = U_{u} + (\frac{d A _{ba se}}{w d t})

U_{d} = U_{u} + (\frac{d A _{ba se}}{w d t})

U_{d} = U_{u} + (\frac{L _{b u b}^{2}}{w τ})

U_{d} = U_{u} + (\frac{L _{b u b}^{2}}{w τ})

z_{d} = z_{u} + L_{b u b}

z_{d} = z_{u} + L_{b u b}

(L_{b u b})_{t + Δ t} = (L_{b u b})_{t} + (\frac{L _{b u b}}{τ})_{t} Δ t

(L_{b u b})_{t + Δ t} = (L_{b u b})_{t} + (\frac{L _{b u b}}{τ})_{t} Δ t

U_{d} = U_{u} + (\frac{L _{b u b}}{τ})

U_{d} = U_{u} + (\frac{L _{b u b}}{τ})

z_{d} = z_{u} + L_{b u b}

z_{d} = z_{u} + L_{b u b}

\frac{d L _{b u b}}{d t} = - U_{u}

\frac{d L _{b u b}}{d t} = - U_{u}

\dot{V}_{in} + \frac{d V _{c}}{d t} = A_{cs} U_{i}

\dot{V}_{in} + \frac{d V _{c}}{d t} = A_{cs} U_{i}

P_{c} V_{c} = P_{0, ini} V_{0, ini} = co n s t an t

P_{c} V_{c} = P_{0, ini} V_{0, ini} = co n s t an t

\frac{d V _{c}}{d t} = - \frac{V _{c}}{P _{c}} \frac{d P _{c}}{d t}

\frac{d V _{c}}{d t} = - \frac{V _{c}}{P _{c}} \frac{d P _{c}}{d t}

P_{0} = P_{1} + (1.5 K) \frac{ρ _{l} U _{i}^{2}}{2}

P_{0} = P_{1} + (1.5 K) \frac{ρ _{l} U _{i}^{2}}{2}

\frac{d P _{0}}{d t} = \frac{d P _{1}}{d t} + (1.5 K) U_{i} \frac{d U _{i}}{d t}

\frac{d P _{0}}{d t} = \frac{d P _{1}}{d t} + (1.5 K) U_{i} \frac{d U _{i}}{d t}

\frac{d P _{1}}{d t} = - \frac{P _{c}}{V _{c}} \frac{d V _{c}}{d t} - (1.5 K) U_{i} \frac{d U _{i}}{d t}

\frac{d P _{1}}{d t} = - \frac{P _{c}}{V _{c}} \frac{d V _{c}}{d t} - (1.5 K) U_{i} \frac{d U _{i}}{d t}

\frac{d P _{1}}{d t} = - \frac{P _{0, ini}}{V _{0, ini}} (A_{cs} U_{i} - \dot{V}_{in}) - (1.5 K) U_{i} \frac{d U _{i}}{d t}

\frac{d P _{1}}{d t} = - \frac{P _{0, ini}}{V _{0, ini}} (A_{cs} U_{i} - \dot{V}_{in}) - (1.5 K) U_{i} \frac{d U _{i}}{d t}

Δ P_{v i s} = \frac{2 ρ _{l} f z _{n, u} U _{i}^{2}}{D _{h}} + i = m \sum n (\frac{2 ρ _{l} f L _{s l ug, i} U _{s l ug, i}^{2}}{D _{h}})_{i}

Δ P_{v i s} = \frac{2 ρ _{l} f z _{n, u} U _{i}^{2}}{D _{h}} + i = m \sum n (\frac{2 ρ _{l} f L _{s l ug, i} U _{s l ug, i}^{2}}{D _{h}})_{i}

(f R e) = 24 (1 - 1.3553 α + 1.9476 α^{2} - 1.7012 α^{3} + 0.9564 α^{4} - 0.2537 α^{5})

(f R e) = 24 (1 - 1.3553 α + 1.9476 α^{2} - 1.7012 α^{3} + 0.9564 α^{4} - 0.2537 α^{5})

Δ P_{a cc} = ρ_{l} z_{n, u} \frac{d U _{i}}{d t} + i = m \sum n (ρ_{l} L_{s l ug, i} \frac{d U _{s l ug, i}}{d t})_{i}

Δ P_{a cc} = ρ_{l} z_{n, u} \frac{d U _{i}}{d t} + i = m \sum n (ρ_{l} L_{s l ug, i} \frac{d U _{s l ug, i}}{d t})_{i}

N u_{4, \infty} = 8.235 (1 - 2.0421 α + 3.0853 α^{2} - 2.4765 α^{3} + 1.0578 α^{4} - 0.1861 α^{5})

N u_{4, \infty} = 8.235 (1 - 2.0421 α + 3.0853 α^{2} - 2.4765 α^{3} + 1.0578 α^{4} - 0.1861 α^{5})

N u_{3, \infty} = 8.235 (1 - 1.883 α + 3.767 α^{2} - 5.814 α^{3} + 5.361 α^{4} - 2 α^{5})

N u_{3, \infty} = 8.235 (1 - 1.883 α + 3.767 α^{2} - 5.814 α^{3} + 5.361 α^{4} - 2 α^{5})

N u_{3} = \frac{N u _{3, \infty}}{N u _{4, \infty}} N u_{4}

N u_{3} = \frac{N u _{3, \infty}}{N u _{4, \infty}} N u_{4}

N u = 0.00222 R e^{1.09} P r^{0.4}

N u = 0.00222 R e^{1.09} P r^{0.4}

H T C = \frac{k _{l}}{δ _{f i l m}}

H T C = \frac{k _{l}}{δ _{f i l m}}

\frac{δ _{f i l m}}{D _{h}} = C_{δ} (3 \frac{ν _{l}}{U D _{h}})^{0.84} (0.07 (\frac{ρ _{l} D _{h} U ^{2}}{σ})^{0.41})^{- 8} + 0. 1^{- 8}^{- 1/8}

\frac{δ _{f i l m}}{D _{h}} = C_{δ} (3 \frac{ν _{l}}{U D _{h}})^{0.84} (0.07 (\frac{ρ _{l} D _{h} U ^{2}}{σ})^{0.41})^{- 8} + 0. 1^{- 8}^{- 1/8}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Modeling of pressure drop and heat transfer for flow boiling in a mini/micro-channel of rectangular cross-section

Shashwat Jain, Prasanna Jayaramu and Sateesh Gedupudi111Corresponding author, Email: [email protected]

Heat Transfer and Thermal Power Laboratory, Department of Mechanical Engineering, IIT Madras, Chennai, India

Abstract

In the present study, a 1-dimensional model is proposed to estimate the pressure drop and heat transfer coefficient for flow boiling in a rectangular microchannel. The present work takes into account the pressure fluctuations caused due to the confined bubble growth and the effect of pressure fluctuations on the heat transfer characteristics. The heat transfer model considers five zones, namely, liquid slug, partially confined bubble, fully confined (elongated) bubble, partial dryout and full dry-out. The model incorporates the thinning of liquid film due to shear stress at liquid-vapour interface in addition to evaporation. The transient fluctuations in pressure and heat transfer coefficient, along with the time-averaged ones, are verified with the experimental data available in the literature. Heat transfer characteristics with flow reversal caused by inlet compressibility are also presented.

keywords:

Flow boiling; Microchannel; Heat transfer coefficient; Pressure drop; Flow reversal.

1 Introduction

Flow boiling in microchannel has been extensively researched in recent times due to its higher heat removal capacity, caused by its high heat transfer surface area to fluid flow volume ratio coupled with the utilization of latent heat of vaporization. Harirchian and Garimella [1] observed single phase, bubbly region, slug - plug, churn and annular flow patterns based on the operating conditions. The heat transfer capabilities are highly affected by the flow patterns. In general, the heat transfer coefficient associated with the slug-plug or the elongated bubble is higher due to the presence of thin liquid film surrounding the bubble. Mirmanto [2] experimentally investigated flow boiling in a single mini/micro-channel and observed a decrease in heat transfer coefficient with the increase in the local thermodynamic quality. Yin and Jia [3] observed an increasing and then a decreasing trend in heat transfer coefficient with local thermodynamic quality. Sobierska et al. [4] observed single phase, bubble and slug flow during boiling in a narrow channel. The maximum heat transfer coefficient was found in the slug region.

1.1 Transient variation

Jagirdar and Lee [5] investigated the transient local heat transfer coefficient at different locations in a single rectangular channel with dimensions as $0.42\times 2.54\times 25.4\,mm$ during subcooled flow boiling. The synchronized flow visualization showed that the variation in the heat transfer coefficient is due to the passing of liquid, vapour bubble (surrounded by liquid) and dry-out, cyclically over the given location. The increase in the local heat transfer coefficient was observed due to the passing of bubble with evaporating thin film. Similar observations were also reported by Bigham and Moghaddam [6] for flow boiling experiment conducted in a microchannel with nanostructures. A sudden drop was observed in wall temperature during the passing of the bubble. Jafari and Okutucu-Özyurt [7] performed Computational Fluid Dynamics (CFD) simulations using the phase-field method to study the hydrodynamics and heat transfer characteristics of an elongated bubble (present in superheated liquid) inside the microchannel. The results showed high heat transfer coefficient on the bottom and side walls in the confined bubble region. Sun et al. [8] demonstrated and formulated the effect of shear stress on the evaporating thin film inside a microtube. The results showed variation in the heat transfer coefficient due to the cyclic process of elongated bubble growth, dry out and liquid slug. The residence time of each region was also measured and it was found that with the increase in the heat flux, liquid slug time decreases and the elongated bubble or dry region time increases. Transient pressure fluctuations and heat transfer coefficient were observed by Barber et al. [9] for a single microchannel using n-Pentane as the working fluid. The confined bubble was observed due to the small aspect ratio and the bubble expanded both in the upstream and downstream directions. Pressure fluctuations were observed due to confinement of the bubble inside the channel. Rao and Peles [10] experimentally studied transient heat transfer coefficient in multiple flow regimes such as bubbly, slug and annular in a rectangular microchannel. The periodic variation in the heat transfer coefficient was observed with the peak due to explosive bubble growth and presence of evaporating thin film which is followed by a sharp decrease due to the emergence of local dry out.

Rapid confined bubble growth causes pressure fluctuations during flow boiling in mini/micro-channels. These fluctuations in the presence of upstream compressibility such as trapped non-condensables, pump characteristics, subcooled boiling in the upstream region/preheater Gedupudi et al. [11] and interaction between the neighbouring channels (in the case of parallel microchannels), as summarized by Prajapati and Bhandari [12], leads to flow reversal. Liu et al. [13] extensively studied the effect of upstream compressibility on flow boiling. Heat transfer coefficient was found to increase with increase in the compressible volume, for certain range of operating conditions. On the other hand, the inlet pressure was found to be oscillating with lower frequency. On further increasing the compressible volume, there was a decrease in pressure fluctuation frequency. Karayiannis and Mahmoud [14] pointed out the lack of data on the effect of flow reversal on heat transfer coefficient. Kandlikar et al. [15] observed pressure fluctuation and flow reversal due to the interaction between the neighbouring channels in a parallel microchannel heat sink. Kenning and Yan [16] and Yan and Lin [17] observed pressure fluctuations under both negligible and high upstream compressibility conditions. Barber et al. [18] studied bubble dynamics inside a rectangular microchannel with FC-72 as the working fluid and reported a sharp rise in pressure along with flow reversal. The experimental investigations carried by Gedupudi et al. [11] indicates that in the presence of upstream compressibility, the upstream end moved towards the inlet of the channel. On the other hand, no flow reversal was observed when the compressibility was removed. However, pressure fluctuations were observed in both the cases. Yin et al. [19] experimentally observed pressure fluctuations in the slug and bubbly flow patterns. Wang et al. [20] and Singh et al. [21] experimentally observed transient two-phase pressure drop, with the amplitude and frequency dependent on heat flux and mass flux.

1.2 Average heat transfer coefficient

Wang and Sefiane [22] experimentally studied the effect of heat flux, vapor quality and hydraulic diameter on flow boiling heat transfer in a rectangular microchannel using FC-72 as the working fluid. The local heat transfer coefficient initially increases due to thin film evaporation and then decreases due to the emergence of dry patches. Diaz and Schmidt [23] studied the variation of local heat transfer coefficient for different heat fluxes and working fluids. The heat transfer coefficient was found to increase with heat flux, indicating the dominance of the nucleation boiling regime. Temperature fluctuation of high amplitude was also observed for water in the low-quality region. Mirmanto [2] experimentally investigated the local pressure drop and heat transfer coefficient for flow boiling of water in a micro-channel. Non-linear variation in pressure drop was observed along the channel length. Local heat transfer coefficient was found to increase and then decrease along the channel. Effect of local pressure on the estimation of local heat transfer coefficient was also studied.

1.3 Prediction of heat transfer coefficient

Thome et al. [24] proposed a three-zone model consisting of the liquid zone, bubble with the thin liquid film and vapor zone for small tubes supplied with constant heat flux on the boundary and evaluated time-averaged heat transfer coefficient. The model utilized three adjustable parameters for initial film thickness, nucleation frequency and minimum film thickness, which were obtained from the experimental data set as explained by Dupont et al. [25]. This model does not consider the effect of local transient pressure during the bubble growth and the effect of shear stress on the thin film depletion. Wang et al. [26] extended the three-zone model (Thome et al. [24]) to a rectangular channel. The four zones they considered are liquid zone, elongated bubble zone, partial dry-out with liquid present in the corners and full dry-out region. The model was verified with the experiments conducted in a single rectangular channel with constant heat flux supplied through the channel walls. This model eliminated minimum thin film thickness as a parameter. This model also neglected the effect of local pressure conditions and effect of shear stress on the depletion of the thin film. The three-zone and four-zone models mentioned above were used to predict heat transfer coefficient over the entire range of thermodynamic quality. Harirchian and Garimella [27] estimated heat transfer coefficient based on the flow patterns observed during flow boiling. For the bubbly region, the correlation is used. For confined annular flow, the model has been developed based on the conservation of mass, momentum (along with the interfacial stresses) and energy equation. Similarly, three-zone model developed by Thome et al. [24] was utilized for slug-plug region but with a different initial film thickness relationship and the associated parameters. A similar flow pattern based approach was also utilized for the prediction of pressure drop across the channel. Antonsen [28] developed numerical models for the evaporation and condensation of refrigerants in a microtube with the combination of correlations, numerical methods and analytical model. Bubble growth model based on the heat flux is calculated for the cases with and without axial conduction. The model utilized the three-zone model for the estimation of heat transfer coefficient. Effects of different parameters such as wall material, nucleation frequency and mass flux have been studied. Jafari and Okutucu-Ozyurt [29] simulated flow boiling process in a microchannel considering an artificial cavity at the channel wall using Cahn - Hilliard method. The bubble is nucleated from the artificial cavity and it grows up till it reaches the channel height and then departs from the nucleation site. The velocity, temperature, and pressure distributions were observed inside the channel. The heat transfer performance increased by 8 times the single phase value due to the presence of the thin liquid film layer, which, later, upon the dry-out, led to a reduction in heat transfer capability. Periodic fluctuations in heat transfer coefficient were observed based on the bubble location. The heat removal capacity was observed to increase with the increase in the mass flux. Cerza et al. [30] performed an analytical investigation of the bubble growth and estimation of heat transfer. The model assumed a moving liquid layer along with the bubble with non-zero velocity and presence of the developing thermal profile. The model utilized the film thickness for the estimation of the heat transfer coefficient under various bubble velocity conditions. The model estimated the temperature profile of the liquid film below the bubble which caused the bubble growth. Magnini and Thome [31] improvised the three-zone model initially developed by Thome et al. [24] for the two-phase system with refrigerant as a fluid. The improved model considered the time dependent variation of the local film thickness along with the inclusion of thermal inertia of the liquid film. The recirculation of the liquid in between the trapped bubbles is also considered for the liquid slug heat transfer coefficient calculation. Li and Hrnjak [32] proposed a model for the bubble dynamics in a micro-tube under flow reversal condition and its effect on the heat transfer coefficient for R134a.

1.4 Objective of the current work

Most of the heat transfer and pressure drop models in the literature focus on microtubes (circular in cross section). But rectangular microchannel heat sinks, machined in copper or aluminum or etched in silicon, are important from the application point of view, considering the dissipation of high heat fluxes and the ease of manufacture. The four-zone model proposed by Wang et al. [26] for a rectangular channel does not consider partially confined bubble region and the influence of local pressure fluctuations, shear stress and flow reversals on heat transfer characteristics. The current work is aimed at predicting the transient heat transfer coefficient under the influence of local pressure fluctuation caused due to the confined growth of the bubble inside a rectangular micro-channel. Transient and time-averaged pressure drop and heat transfer coefficient are compared with the experimental data available in the literature for a single rectangular microchannel. Flow reversal caused due to explosive bubble growth and presence of upstream compressible volume and its effect on the heat transfer coefficient is studied. Effect of shear stress on transient heat transfer characteristics is also modeled.

2 Model description

The current model consists of two modules i.e., the bubble growth module and the heat transfer module. The first module estimates bubble length, bubble velocity, and pressure drop and the second module calculates the heat transfer coefficient based on the four zones, namely, (1) the liquid zone, (2) elongated bubble zone, (3) partial dry-out zone and (4) full dry-out zone as proposed by Wang et al. [26] and additionally considers partially confined bubble zone, as shown in Fig. 1. The liquid zone consists of the liquid trapped between the two bubbles. The second zone (partially confined bubble zone) consists of liquid bulk on either side the bubble as shown in Fig. 1. The third zone consists of thin liquid film surrounding bubble along with the bulk liquid present in the corner of the rectangular channel between the bubble and the wall. As the thin film depletes the third zone appears with the dry patch area along with the bulk liquid in the corner. Finally, complete evaporation of the corner liquid leads to full dry-out zone/vapor zone.

2.1 Assumptions

In order to simplify the model few assumptions are considered. These assumptions are in line with the assumptions made by Thome et al. [24], Wang et al. [26] and Gedupudi et al. [11].

Current model is 1-dimensional i.e. variations are only considered in the axial direction, i.e., along the flow direction. 2. 2.

The heat sink consist of a single channel with constant heat supplied through the three sides (two sides and the base) of the channel. 3. 3.

The bubble in the elongated region (liquid as thin film and liquid in the corner) is considered as a rectangle with rounded corners for simplification. 4. 4.

Fluid is considered to be saturated from the point of nucleation; this is possible for the cases with the flow boiling in the near saturated condition. 5. 5.

The bubbles are considered to be in the saturated condition, hence properties corresponding to local pressure are considered. 6. 6.

Initial bubble is assumed to be a cube with the length, breadth and height equal to the minimum dimension of the cross-section as assumed by Gedupudi et al. [11]. 7. 7.

Single nucleation site is considered in the present study. In the current case, the bubble grows in axial direction exponentially leading to increased velocity which reduces the thermal boundary layer thickness suppressing the bubble growth. Moreover passing of a long bubble in a shorter channel also decreases the possibility of multiple nucleation sites. The models developed by Thome et al. [24] and Wang et al. [26] also consider single nucleation site. 8. 8.

The length of liquid slug is equivalent to the liquid passing over the nucleation site during the interval between the successive nucleation, i.e., the liquid slug length is equal to the liquid velocity multiplied by the nucleation time period. 9. 9.

Minimum film thickness before the occurrence of the dry-out is of the order of surface roughness. This value varies with surface, so, if available, the specified value would be used. Else, it would be considered as a constant similar to the value proposed by Thome et al. [24]. 10. 10.

Change in property of fluid due to viscous dissipation as explained by Han and Lee [33] does not affect the properties of liquid significantly. 11. 11.

The simulation ends if there is no bubble in the channel, or the upstream end of the bubble reaches the channel entrance, possible in the case of inlet compressibility. 12. 12.

The contribution of surface tension to the pressure drop is negligible. The surface tension of water at [math]C is 0.058 N/m and the minor dimension of the channel considered in the study is 0.00024 m. Therefore the approximate pressure drop due to surface tension is $\sigma/R=0.058/0.00012=0.48\,kPa$ , which is negligible compared to the pressure drop due to friction and acceleration for the range of the parameters considered in the present study.

2.2 Properties determination

Variation of fluid properties with pressure is not considered in the models for the bubble growth presented by He et al. [34] and Gedupudi et al. [11]. Thome et al. [24] and Wang et al. [26] have also not considered the fluid property variation with local transient pressure in their heat transfer models. In the case of water, there is a significant difference in the amplitudes of pressure fluctuation if the variable properties are not considered (presented in the current study). This is due to the fact that as the pressure drop increases across the channel, the local pressure changes, which leads to a change in vapor properties that defines the bubble growth rate. These, in turn, influence the local velocity and pressure drop. The fluid properties are taken from Bell et al. [35]. The saturated properties corresponding to the local pressure are estimated at each time step. The pressure at the channel outlet is assumed to be 101 kPa. In the current study, the working fluid is water unless otherwise specified.

2.3 Location of nucleation site

The current transient model requires a nucleation site as an input parameter. This is one of the major challenges in the area of flow boiling as nucleation site is not deterministic and it depends on surface characteristics. Therefore, it is difficult to estimate the nucleation site accurately. Though, Liu et al. [36] summarizes relationships available in the literature and proposes the condition for nucleation based on wall superheat. In the present study, Eq. (1) is used to determine nucleation site.

[TABLE]

The location of nucleation site is determined iteratively. Initially for the specified inlet temperature and for the assumed channel inlet pressure equal to the exit pressure, the location where the wall temperature obtained from single-phase heat transfer coefficient correlation matches the wall temperature obtained from Eq. (1) is chosen as the initial nucleation site. Time-averaged two-phase pressure drop is then calculated. Saturation temperature corresponding to the obtained pressure at the assumed nucleation site is determined, and then the new wall temperature required for nucleation from Eq. (1) is obtained. Then the new location of the nucleation site is determined using single phase correlation based on the required wall temperature for nucleation. The estimation of the location of nucleation site and two-phase pressure drop is performed iteratively till the required wall temperature from Eq. (1) matches the obtained wall temperature from single phase correlation. Contact angle varies with the substrate and should be chosen accordingly. For example, it is around [math] for copper.

2.4 Nucleation frequency

The nucleation time period consists of two components i.e., waiting time and growth time. Bubble growth and its departure from the nucleation site disturbs the local wall temperature. It takes a while for the reformation of the thermal boundary layer, i.e., for the wall to attain the superheat required for nucleation. This duration is called as waiting period. This has been determined with the help of transient single phase simulation carried out using ANSYS FLUENT 16.1. Another component is bubble growth period and is determined using the correlation given by Mikic et al. [37]. The initial theoretical value obtained using this correlation is quite small and is of the order of $100\,\mu s$ . The experimental data obtained from Kadam and Kumar [38], Lee et al. [39] and Wang et al. [26] shows that the nucleation time period varies in the range 1-500 ms for the water based on both heat flux and mass flux. In the present case, with the increase in heat flux, the time period decreases as both the waiting and bubble growth time decreases. As the heat flux increases the wall reaches the condition of superheat at a faster rate. Similarly, with the increase in mass flux two forces are opposing one is the inertial force which detaches the bubble from the nucleation site this decreases the time period as also observed by Jafari and Okutucu-Ozyurt [29] and another is that the increases the time periods as due to large mass flux the thermal boundary layer thickness decreases.

[TABLE]

where,

[TABLE]

Here,

$C1$ is the adjustable factor to determine the growth time period. In the present study $C1=210$ is used.

2.5 Bubble growth

In the case of the flow boiling the bubble is initiated at the channel wall based on the wall superheat requirement. The initial stage consists of the free (unconfined) growth similar to the nucleate boiling. Once the bubble grows to the minor dimension of the channel the growth rate changes as observed by Yin and Jia [3]. Barber et al. [18] and Gedupudi et al. [11] also observed two different growth stages, namely partially confined growth where the bubble is confined by the minor dimension and grows widthwise and lengthwise and fully confined growth where the bubble is confined by the minor and major dimensions of the cross-section and the bubble grows only along the length. The stage where the vapour bubble passes through the channel outlet is called vapour venting. Energy balance for the bubble growth is given in Eq. (3).

[TABLE]

In the partially confined region, the area $A_{base}$ differs from the partially confined condition and fully confined condition. $A_{base}=L_{bub}\times L_{bub}$ for partially confined stage. $A_{base}=L_{bub}\times w$ for fully confined stage. $A_{heat}=L_{bub}\times L_{bub}$ for partially confined stage and $A_{heat}=p_{h}\times L_{bub}$ for fully confined stage. $p_{h}$ is the heated perimeter equal to $w$ for single-sided heating and equal to $2h\,+\,w$ for three-sided heating.

Time constant

Time constant is defined in Eq. (4).

[TABLE]

where $p_{h}$ is the heated perimeter. The value of $\tau$ is calculated for each time step and substituted in the equations to obtain the downstream end velocity, position and the bubble length.

2.5.1 Partially confined bubble growth

The growth rate is obtained by substituting the expression for $A_{base}$ in Eq. (3) and is given as

[TABLE]

Using Taylor series expansion, length of the bubble at t+ $\Delta$ t is estimated as in Eq. (6) where ‘t’ refers to the time from the start of the simulation.

[TABLE]

From continuity, the velocity of the downstream end of the bubble is obtained as shown in Eq. (7) and (8).

[TABLE]

Similarly, the location of the downstream end of the bubble is calculated by Eq. (9).

[TABLE]

2.5.2 Fully confined bubble growth

Similar to the partially confined bubble growth, the length, the velocity and the position of the bubble in fully confined growth are determined using Eq. (10), (11) and (12) respectively.

[TABLE]

2.5.3 Vapour venting

During vapour venting stage, as the bubble passes through the channel outlet, the bubble is driven by the velocity on the upstream end and is given by

[TABLE]

The exit vapour velocity is determined using either Eq. (8) or (11) as the case may be.

2.6 Upstream compressibility

Due to the presence of sources of upstream compressibility, there can be flow reversal in the channel. In the present study, trapped non-condensable in the channel inlet plenum is considered, with pressure $P_{c}$ and volume $V_{c}$ . Plenum is supplied with constant flow rate $\dot{V}_{in}$ . The channel inlet velocity can be calculated from the principle of continuity as.

[TABLE]

For the current case, non condensable gas is assumed to undergo isothermal compression.

[TABLE]

The rate of change of gas volume is determined using Eq. (16).

[TABLE]

Stagnation pressure $(P_{0})$ is given by

[TABLE]

Here $P_{c}$ is equal to $P_{0}$ , and $U_{i}$ is channel inlet velocity and $P_{1}$ is channel inlet pressure. $K=1$ for the forward flow and $K=0$ for the flow reversal.

Differentiating the pressure equation Eq. 17 gives

[TABLE]

Substituting the Eq. 16 in the Eq. 18.

[TABLE]

Substituting the value of $dV_{c}/dt$ from the Eq. 14.

[TABLE]

The initial compressible volume of the non-condensable gas is assumed to be a variable parameter and would be varied in order to study its effect.

2.7 Pressure drop

The channel pressure drop consists of the viscous pressure drop and the acceleration pressure drop across the liquid slugs. Pressure drops due to momentum change and friction across vapour bubble or vapour slug are assumed to be negligible Gedupudi et al. [11] as $\rho_{v}<<\rho_{l}$ . As explained in the Assumptions section, pressure drop across the vapor bubble and pressure drop due to surface tension are neglected.

2.7.1 Viscous pressure drop

The frictional pressure drop arises due to the motion of liquid slugs present in the channel. The total frictional pressure drop is the sum of the frictional pressure drops calculated for the individual liquid slugs. Apart from this, there will also be viscous pressure drop due to the liquid slug between the channel inlet and upstream end of the newest bubble.

[TABLE]

where,

where ‘m’ represents the liquid slug in between the exit and the last bubble, ‘n’ represents the liquid slug between the the newest bubble ‘n’ and ‘n-1’ bubble, as shown in Fig. 2.

The Fanning friction factor f used in Eq. (21) is for a rectangular cross section with the aspect ratio $\alpha$ and is obtained from Eq. (22). If the value f goes below 0.01, then f is taken as a constant equal to 0.01 to account for the turbulent flow, as considered in the previous study by Gedupudi et al. [11].

[TABLE]

2.7.2 Acceleration pressure drop

The acceleration pressure drop is calculated using the momentum conservation applied for all the liquid slugs.

[TABLE]

The first term in the above expression will be zero if there is no upstream compressibility.

2.8 Heat transfer coefficient

In the present work, the time-averaged heat transfer coefficient is calculated at each location. The heat transfer coefficient is calculated for each zone i.e. the liquid zone, partially confined bubble zone, elongated bubble zone, partial dry-out zone and full dry-out zone individually as they pass over a particular position cyclically. These are then time-averaged based on their residence time to obtain the local time-averaged heat transfer coefficient. The estimation of the heat transfer coefficients for different zones is described below.

2.8.1 Liquid slug

This is the liquid present between the two successive bubbles. The liquid slug moves with the velocity equal to the downstream end velocity of the bubble present on the upstream end or with the inlet velocity. In the literature various correlations have been used to estimate the heat transfer coefficient. Thome et al. [24] estimates the liquid slug heat transfer coefficient (HTC) with the help of the laminar and turbulent correlations for the circular tube case. In order to avoid the jump from the laminar to turbulent region, asymptotic approach is used. This method averages the laminar and turbulent region HTC. In the case of four zone model developed by Wang et al. [26], single phase laminar fully developed region correlation is used. He et al. [40] studied the heat transfer characteristics of liquid slug trapped between the two bubbles. Due to the movement of the liquid and no slip condition at the wall, internal vortex is seen. This vortex increases the convection and therefore increases the HTC. In the present study the fluid flow is considered laminar for $Re<\,1600$ and in the transition zone if $Re$ is in between 1600 and 3000 and turbulent for $Re\,>\,3000$ as used by the Kandlikar and Balasubramanian [41]. For the transition range, based on $Re$ , linear average of the heat transfer coefficient is used. Eq. (26) is used for the calculation of HTC for the rectangular channel in laminar region with the aspect ratio $(\alpha=h/w)$ .

For the laminar fully developed condition, the Nusselt number expressions for four-sided and three-sided heating configurations are as given in Eq. (24) and Eq. (25) respectively.

[TABLE]

Lee and Garimella [42] suggested heat transfer correlation for laminar thermally developing region for the rectangular channel with the aspect ratio $(\alpha)$ for the constant heat flux condition for four sided heat condition which was further modified by Phillips [43] for three side heating condition as given in Eq. (26).

[TABLE]

In the case of turbulent flow the relation, Eq. (27), proposed by Wu and Little [44] for a microchannel is used.

[TABLE]

2.8.2 Thin film region

As elucidated in the literature, evaporation of the thin film is the major contributor of the heat transfer coefficient. Hence it is very important to estimate the thin film thickness. Initial film thickness plays an important role in this process. Various researchers have proposed the relationship based on the experimental observations. Thome et al. [24] used the basic relationship provided by Moriyama and Inoue [45] and used an asymptotic method to estimate the film thickness. A factor was also used to correct the above relationship based on the experimental data available in the literature. Wang et al. [26] also used the similar equation as by Thome et al. [24] but the factor is based on the experiments conducted. Han and Shikazono [46] proposed the relationship for the initial film thickness for the bubble accelerating through the channel. The initial film thickness decreases with the increase in the acceleration due to the thinning of boundary layer. Film thickness at any location or at any instant of time is influenced by initial film thickness and rate of thinning. The latter depends on the rate of evaporation and shear stress. Thome et al. [24] and Wang et al. [26] considered only the evaporation of the thin liquid film, deducing the new film thickness with the help of the energy balance. On the contrary, Consolini [47] considered the effect of shear stress on the film thickness for a bubble growing inside a circular tube. Sun et al. [8] considered the thinning of the film due to evaporation and proposed an expression for the additional thinning of the thin film due to the shear stress. The relationship was derived analytically and a constant was used in order to match the experimental film thickness. Younes and Hassan [48] also employed the analytical approach by considering the shear stress at the liquid-vapour interface to calculate the heat transfer coefficient. The heat transfer coefficient is calculated by Eq. 28.

[TABLE]

The initial film thickness Eq. (29) for a vapour bubble in the two phase region is:

[TABLE]

Thome et al. [24] obtained the value of the constant $C_{\delta}$ equal to 0.29, by matching the model results with the experimental data for refrigerants boiling in a circular tube. The initial film thickness based on this constant is of the order of few microns. Wang et al. [26] used a value equal to 0.8117 by matching their model results with their experimental heat transfer data for water boiling in a rectangular channel ( $D_{h}$ = 0.137 mm). The initial film thickness based on this constant is in the range $1-2\,\mu m$ for the bubble velocity in the range $0.2-0.5\,m/s$ considered in their study. But in the present study with water as a working fluid, for the range of the parameters considered ( $D_{h}$ = 0.324 mm and 0.658 mm), the value of the constant $C_{\delta}$ equal to 3.48 results in a good match between the modeling results and the experimental heat transfer data from the literature. The initial film thickness based on this constant is in the range $10-30\,\mu m$ for the bubble velocity in the range $0.5-2\,m/s$ . The initial film thickness predictions made using this constant is in line with the experimental thin film thickness measured by Sun et al. [8] for water. Their experimental observations for flow boiling indicate initial film thickness of 25 and 31 $\mu m$ for the bubble velocities 3.09 m/s and 3.63 m/s and tube diameter 0.94 mm. The corresponding values obtained using the proposed constant (equal to 3.48) is 17 and 16 $\mu m$ respectively. Further, their experimental observations show that from the point of initial film thickness, there is a gradual decrease in local film thickness until it reaches dryout thickness which is around 2-3 $\mu m$ . The film thickness estimated by Gedupudi et al. [11] is in the range 20-30 $\mu m$ and that observed by Balasubramanian and Kandlikar [49] is near 50 $\mu m$ , all for water boiling in a mini/microchannel close to atmospheric pressure.

In the present study, both the heat transfer and shear stress effect are taken into account to evaluate the depletion of the thin film thickness. Thin film exists in partially confined bubble region, elongated bubble region and partial dryout region.

Partially confined bubble region

In this region, the bubble is partially confined by the minor channel dimension. The heat transfer coefficient is calculated considering the heat transfer through the thin film between the bottom wall and the bubble and the heat transfer through the liquid on either side of the bubble. The heat transfer coefficient is calculated using Eq. (30) to Eq. (32). The liquid heat transfer coefficient is estimated using the equations presented in Liquid slug section considering the effective liquid flow cross-section obtained after subtracting the bubble cross-sectional area from the channel cross-section. Hydraulic diameter is calculated using the effective cross-sectional area and the corresponding wetted perimeter.

[TABLE]

The $HTC_{film}$ is given by the Eq. 33.

[TABLE]

The $\delta_{film}$ is the film thickness over the flat region and its depletion is given by Eq. 34.

[TABLE]

The $\delta_{eq}$ is the equivalent thickness over the horizontal region which is given by Eq. 35.

[TABLE]

$A_{cs,bub}$ is calculated for hyperellipse with $n1=4$ , $a=0.5\times L_{bub}$ and $b=0.5\times h-\delta_{film}$ in Eq. 36 presented by Tamayol and Bahrami [50].

[TABLE]

Here, $\Gamma$ is the gamma function.

Elongated bubble region

This is the region where the cross- section has liquid film over the entire perimeter of the channel as shown in Fig. 4. Vapour bubble is assumed to be of rectangular shape with rounded corners, obtained with $n1=4$ in hyperellipse formula Tamayol and Bahrami [50]. The region between the curved portion at the corner and the channel walls form the bulk liquid region.

Here $\delta_{film}$ is the film thickness at the flat region of the channel.

[TABLE]

Here $p_{lv}$ is the perimeter of liquid and vapor interface at any cross-section used by Tamayol and Bahrami [50]. $p_{ch}$ is the channel perimeter and $p_{h}$ is the heated perimeter.

[TABLE]

The cross-sectional area for a rectangle with rounded corner is given by following equation and used by Tamayol and Bahrami [50]. Here $n1=4$ for rectangle with rounded corners. $a=0.5\times w-\delta_{film}$ and $b=0.5\times h-\delta_{film}$

[TABLE]

Here, $\Gamma$ is the gamma function.

The Eq. 40 is used to determine the heat transfer coefficient in the elongated bubble region.

[TABLE]

Partial dryout region

Once the thin film in the flat region of the elongated bubble is evaporated fully, i.e., $\delta_{film}=0$ , it leads to partial dryout as shown in Fig. 5. The corners of the channel still have the liquid which is not evaporated. Thus the equivalent heat transfer coefficient is the average of the heat transfer coefficient of the dryout region i.e. over the wall in contact with vapour and liquid heat transfer coefficient due to liquid in the corners of the channel. In the present case, the cross-sectional area occupied by vapour bubble is of rectangular in shape with rounded corners. In partial-dryout region, perimeter of bubble cross-section intersects with the channel wall. and the corner liquid area and lengths are calculated.

[TABLE]

Here $l_{l}$ refers to the length of the wall in contact with the liquid over the heated perimeter. $l_{l}\,=\,l_{a}\,+\,l_{b}\,+\,l_{c}=l_{a}\,+\,2l_{b}$ .

[TABLE]

$A_{l}$ is the area enclosed by the sides of the channel and the curve formed by the liquid-vapour interface at the corner. This curve is obtained by assuming a hyperellipse (Tamayol and Bahrami [50]) with $(0.5\times w+\delta_{corner})$ and $(0.5\times h+\delta_{corner})$ as major axis (a) and minor axis (b), as described in Tamayol and Bahrami [50]. The enclosed area is calculated by numerical integration.

[TABLE]

And the corner bulk liquid heat transfer coefficient is calculated by the Eq. (46).

[TABLE]

2.8.3 Full dry-out region

This is the region developed due to the complete evaporation of the thin film. The HTC can be calculated by means of the Eq. (24) to (27) using vapor properties.

2.8.4 Shear stress effect

Sun et al. [8] experimentally investigated the evolution of thin film with time, for adiabatic and diabatic cases. The depletion of the film thickness was observed to be faster in diabatic case than that in adiabatic case. They also derived the depletion rate of thin film thickness for circular tube. In the present study, effect of shear stress is extended to a rectangular channel.

From continuity,

[TABLE]

Differentiating the above equation with time and using $dz=U_{l,mean}dt$ results in Eq. (49). The rate of thin film depletion is defined by the Eq. (49).

[TABLE]

This equation requires the variation of the mean liquid velocity with the axial position. In order to estimate the value, the shear stress on the vapour side is equated to the shear stress on the liquid side at liquid-vapour interface, which is given by Eq. 50.

[TABLE]

The shear stress on the vapor side is given by the following Eq. (51).

[TABLE]

Similarly the shear stress due on the liquid is given by the following Eq. (52).

[TABLE]

Equating the shear stress relationship mentioned in Eq. (51) and (52) to obtain Eq. (53).

[TABLE]

In the Eq. 53, the spatial variation of the vapor velocity needs to be determined.

[TABLE]

In the Eq. (54), $u^{*}$ is a non-dimensional velocity calculated as described by Tamayol and Bahrami [50]. The equation below shows the relationship between the mean velocity and maximum velocity with the help of the velocity distribution along the cross-section.

[TABLE]

Differentiating Eq. (53) w.r.t dz results in the Eq. (56).

[TABLE]

Substituting the value of $U_{max}$ from Eq. (55) into Eq. (56) results in Eq. 57.

[TABLE]

From the energy balance the mean velocity variation vapor velocity is given by the Eq. (58).

[TABLE]

$dU_{v,mean}/dz$ from Eq. (58) is substituted into Eq. (57) and the resulting value of $dU_{v,mean}/dz$ into Eq. (49) to obtain the rate of depletion of the thin film due to shear stress.

In the case of partial dryout region , liquid is present in the corners and the effect of the shear stress on the depletion of the thin film is modeled as described below.

[TABLE]

$dU_{l}/dz$ is calculated using the Eq. (53) to (56). $dU^{*}/dr$ is calculated using hyperellipse with exponent as 4 and major and minor axes as mentioned in the section on partial dryout region. The limits of $\theta$ and $r$ in Eq. (55) and (57) correspond to intersection made by the hyperellipse with the rectangle.

In the present study, for pressure drop module, a single vapour bubble is assumed that comprises of all stages-partially confined, elongated, partial dryout and full dryout-considered in the heat transfer module. This will not affect pressure drop calculation as pressure drop across vapour bubbles is negligible, as mentioned in Pressure drop section. But these stages do influence the heat transfer coefficient and hence are considered in heat transfer module.

3 Solution procedure

In order to solve the equations presented in the bubble growth and heat transfer coefficient sections the channel is divided into several discretized elements with the length as dz. Forward time step marching is used with time step dt. Sequence of steps followed to obtain the solution is given below.

The first step involves the determination of the location of nucleation site as explained in section 2.3. 2. 2.

The nucleation time period is calculated by the procedure described in section 2.4. 3. 3.

Based on the time period and nucleation site, a partially confined bubble of dimension $h\times h\times h$ is placed at the nucleation site, as mention in Assumptions section. 4. 4.

Velocity at the downstream end of the bubble is calculated using Eq. (8) or (11) and bubble length based on Eq. (6) or Eq. (10) depending on the on the bubble growth stage. 5. 5.

Then frictional and acceleration pressure drops are calculated using Eq. (21) and Eq. (23) respectively. 6. 6.

Inlet stagnation pressure is calculated and then the inlet velocity using Eq. (17) if inlet compressibility is present. 7. 7.

Local heat transfer coefficients are calculated based on the zones as explained in section 2.8, considering the fluid property variation with the local transient pressure. 8. 8.

Time step is then incremented and the steps (4)-(7) are repeated until periodic behavior is achieved. 9. 9.

Time averaged and channel averaged pressure drop and heat transfer coefficient are calculated.

4 Results and discussion

4.1 Validation and verification

4.1.1 Time step and grid size independence study

Grid independence and time-step independence studies are carried out to make sure that the numerical errors do not affect the solution. Transient heat transfer coefficient and pressure drop variations are observed and it can be seen from Fig. 6 and 7 that the change is negligible for the grid size and time step smaller than $2\times 10^{-5}$ m and $2\times 10^{-5}$ s respectively.

4.1.2 Transient pressure and heat transfer coefficient fluctuations

Comparison of the modeled pressure fluctuation with the experimental data reported by Brutin and Tadrist [51] for n-Pentane is shown in Fig. 8a. The amplitude predicted by the model closely matches the experimental observation. The nucleation time period used in the model is the same as the observed fluctuation period. The nucleation time period may also depend on the surface characteristics which has not been taken into account in the present study. The location of the nucleation site is calculated as described in section 2.3. The modeled transient pressure has also been compared with the CFD prediction made by Zu et al. [52] for water and is shown in Fig. 8b. The trend matches with the 3-dimensional CFD result, but the amplitude obtained from model is lower than that predicted by CFD study due to the fact that the modeling considers fluid property variation with local pressure, but the CFD study assumes constant fluid property. The increase in pressure drop is due to the acceleration of the liquid slugs and the corresponding frictional pressure drop. As the liquid slug leaves the channel, the amount of accelerated liquid slug present in the channel decreases and hence the drop in pressure.

Fig. 9 shows the comparison of the modeled local transient heat transfer coefficient with the experimental data of Jagirdar and Lee [5] for water, under different operating conditions. The modeled heat transfer coefficient is in reasonably good agreement with the experimental data, both in trend and magnitude.

The local heat transfer coefficient variation from A to B is due to the acceleration of the liquid slug caused by the bubble growth on the upstream side. The jump from B to C is due to the change in zone from liquid slug to vapour bubble that causes a sudden increase in heat transfer coefficient. The gradual increase in heat transfer coefficient from C to D is attributed to the reduction in the thin film thickness caused by evaporation and shear stress. The drop from D to E is due to the change in zone from vapour bubble to liquid slug at that location. The heat transfer coefficient remains constant from E to F, as the liquid slug with constant velocity passes over the location. This continues until a bubble nucleates, which then causes an increase in the liquid velocity and hence the heat transfer coefficient (from A to B).

4.1.3 Time averaged pressure drop and heat transfer coefficient

In this section, time-averaged pressure drop and heat transfer coefficients are compared with the experimental data reported by Jayaramu et al. [53]. Fig. 10 shows the comparison for a channel of dimension $0.24\times 0.50\times 40$ mm, for the specified operating conditions. The modeled time-averaged pressure drop is within 25% of the experimental average pressure drop, as shown in Fig. 10a, and the time-averaged heat transfer coefficient within 24%, as in Fig. 10b. The maximum deviations for a channel of dimension $0.49\times 1.00\times 40$ mm are 20% and 5% for pressure drop and heat transfer coefficient respectively, as shown in Fig. 10. The deviations may perhaps be attributed to the change in flow patterns.

In the present model, the calculated quality is different from the thermodynamic quality. Liquid slugs are assumed to be superheated due to the heat transfer from the wall and the heat conducted between the liquid slug and the vapour bubble is neglected as the area of contact between the liquid slug and vapour bubble is assumed to be small compared to the total heat transfer area of liquid slug. The quality is estimated by the Eq. 60.

[TABLE]

The summation is over the time period at the channel exit. The variation in quality with heat flux is presented in Fig. 12. There is a reduction in the percentage of deviation between the thermodynamic quality and the actual quality as the heat flux increases. This is due to the decrease in the amount of liquid slug with the increase in frequency caused by the increase in heat flux. There is also a reduction in the percentage of deviation with the decrease in mass flux, due to the decrease in the amount of liquid slugs. Magnini and Thome [54] performed CFD simulations for confined bubble growth in a microtube under constant heat flux condition and reported 82% utilization of the heat flux for the bubble growth, which perhaps indicates the presence of superheated liquid and hence the deviation between the actual and thermodynamic quality. The observations from the present study are in line with their results.

Fig. 13 makes a comparison between the spatio-temporally averaged heat transfer from the model and that obtained from the existing steady flow correlations (Table 1) in the literature. For lower mass flux, the model prediction lies in between the ones obtained using Sun and Mishima [55] model and Li and Wu [56] model, and for higher it is in between Sun and Mishima [55] and Mahmoud and Karayiannis [57]. The possible reasons for deviation are the conditions used for the development of correlations and the influence of heating surface characteristics which has not been taken into account in the present model and also in the development of the correlations.

4.2 Different zones

The occurrence of five different zones during flow boiling depends on the channel dimension, inlet condition, mass flux and heat flux. Fig. 14a shows different zones passing cyclically over the location that is 1.25 mm away from the inlet. Initially the transients are affected by the initial condition during simulation and later the cycle repeats itself. Zone A-B is the liquid slug zone (un-accelerated), B-C is the liquid slug region (accelerated), the jump C-D represents transition from liquid to partial confinement, D-E the partially confined bubble region, E-F the jump from partial confinement to full confinement, F-G the fully confined or elongated bubble stage and G-H the transition from the bubble to the liquid slug zone. Fig. 14b shows different zones over a location that is 3.75 mm from inlet, for with and without shear and for different minimum liquid film thickness values for dryout. Zones A-B, B-C, C-D and D-E indicate the liquid slug zone (un-accelerated), the liquid slug zone (accelerated by partially confined bubble), the liquid slug zone (accelerated by fully confined bubble) and the transition from the liquid slug to the elongated bubble region, respectively. The zones E-F, E-F’ and E-F” represent elongated bubble region for the case with shear (3 $\mu m$ ), with shear (0.3 $\mu m$ ) and without shear (3 $\mu m$ ) respectively. The minimum liquid film thickness values considered in the present study are based on the experimental observations by Sun et al. [8]. Zones F-G, G-H, H-I, I-J and J-K represent the transition to partial dryout, the partial dryout region, the transition to full dryout, the full dryout region and the transition to the liquid slug region, respectively. Zones F’-K and F”-K indicate the transition from elongated bubble to liquid slug region. The differences between the two figures are the presence of partially confined bubble region and the absence of partial and full dryout in Fig. 14a as against Fig. 14b. The effect of shear stress is discussed in detail in the section 4.3.

4.3 Shear stress effect

In the process of bubble growth, the thin liquid film gets depleted due to the evaporation of film and also due to the shear stress acting on the liquid-vapor interface. The models presented by Thome et al. [24] and Magnini and Thome [31] for circular channel and Wang et al. [26] for rectangular channel, do not take into account the influence of shear. Fig. 15 shows the effect of shear stress on the variations of heat transfer coefficient and thin film thickness. Dashed lines indicate the heat transfer coefficient with shear and solid lines without shear. The spatially averaged and the time-averaged heat transfer coefficients with shear stress are higher than that without shear, as shown in Fig. 15a and Fig. 15b, respectively. This can be explained from Fig. 15c that shows the variation of local heat transfer coefficient with time for two different locations ( $0.2\,\times\,L_{ch}$ and $0.4\,\times\,L_{ch}$ from the inlet). With shear, the heat transfer coefficient during the passage of elongated bubble is higher compared to that without shear due to the higher depletion rate of thin film with shear as shown in Fig. 15d. The reduction in the thin film thickness is larger for the upstream location as the rate of depletion of thin film is directly proportional to the film thickness as indicated by the Eq. (49). This leads to a larger increase in the heat transfer coefficient for the upstream location.

4.4 Effect of inlet compressibility

This section presents the influence of inlet compressibility caused by trapped non-condensable gas on pressure fluctuations, flow reversals and heat transfer coefficients. Modeling the location of nucleation site and nucleation frequency with flow reversals is difficult due to the complex thermo-fluidics and hence the approaches presented in Sections 2.3 and 2.4 are not valid for the cases with flow reversals. Here, for simplicity, the nucleation site is assumed to be in the middle of the channel to account for the bubble growth in the upstream direction during flow reversal. Bubble is assumed to nucleate at regular intervals provided there is no bubble passing over the nucleation site at that instant of time. With inlet compressibility, there is a flow reversal, i.e., the bubble grows in both the upstream direction and the downstream direction as shown in Fig. 16a. The flow reversal can be explained using Eq. (14) and (20). The flow, after reversing for certain duration, changes its direction and moves forward. The extent of flow reversal is coupled with the change in the inlet (plenum or stagnation) pressure. Fig. 16b shows that there is a reduction in the amplitude of pressure fluctuation with the increase in compressible volume, which is in line with the experimental observations made by Liu et al. [13]. This can be attributed to the decrease in the net acceleration (and hence velocities) due to flow reversal. Fig. (16c) indicates higher amplitudes of heat transfer coefficient fluctuation and enhanced average heat transfer coefficient in the presence of instabilities. The spatial and temporal averaged values of heat transfer coefficient are 30.8, 37.0 and 39.3 $kW/m^{2}K$ for no compressible volume, 4 $V_{ch}$ and 8 $V_{ch}$ respectively. The maximum heat transfer coefficient is higher for the case with flow reversal due to larger residence period of the confined bubble leading to larger contribution from the thin film evaporation as shown in Fig. 17a to 17d. The minimum heat transfer coefficient with flow reversal is lower than that without flow reversal, due to lower velocities (caused by reduction in the net acceleration) of liquid slugs with flow reversals. Fig. 17a and 17b show the variation of local heat transfer coefficients with time on the upstream side of the nucleation site and Fig. 17c and 17d show the variation on downstream side of the nucleation site. It is clear that the thin film evaporation exists for a longer duration over a larger length for the cases with flow reversals, resulting in higher heat transfer coefficients. Higher heat transfer coefficients associated with the flow reversal cases are in line with the experimental observations made by Gedupudi et al. [11]. Kenning and Yan [16] reported no appreciable change in the heat transfer coefficient due to compressible volume for the range of parameters and dimensions considered in their experimental study. Liu et al. [13] reported higher heat transfer coefficients with compressible volume for saturated boiling of water for certain heat flux and mass flux and enhanced heat transfer coefficient only at thermodynamic quality greater than 0.05 for a different heat flux and mass flux. Li and Hrnjak [32] carried out the experimental investigation of the influence of periodic flow reversal on the performance of microchannel heat exchanger and reported higher heat transfer coefficients, mainly in the upstream part. The results obtained from the present model, shown in Fig. 17a and 17b, are very similar to their experimental observations.

Fig. 18a and Fig. 18b show the local transient variation of heat transfer coefficient and thin film thickness respectively. Larger residence time of the bubble for the case of flow reversal, as shown in Fig. 18b, causes complete dryout on the flat surfaces leaving the liquid film only at the corners, termed as partial dryout that results in lower local heat transfer coefficient as shown in Fig. 18a. In Fig. 18b, A-B indicates transition from liquid slug to elongated bubble, B-C depletion of thin film on flat face in elongated bubble region, C-D transition to partial dryout and D-E partial dryout. Tuo and Hrnjak [58] observed decrease in heat transfer coefficient due to the emergence of dryout during flow reversal. C’-D’ denotes transition from vapour bubble to liquid slug. The film thickness at C reaches minimum film thickness whereas at C’ it is higher than minimum film thickness. A-E and A’-D’ indicate the bubble residence period for flow reversal case and no flow reversal case respectively.

4.5 Effect of local pressure

In this section, the effect of local pressure on the pressure drop and the heat transfer coefficient is studied. In the models previously developed by Kew [59], Wang et al. [26], Thome et al. [24] and He et al. [34], fluid property variation with local transient pressure was not considered. Fig. 19a makes a comparison between the transient variation of channel pressure drop calculated considering the property variation with local pressure and assuming constant pressure equal to exit pressure. The peak pressure obtained considering the property (vapour density, latent heat, liquid density and liquid viscosity) variation with local pressure is much lower and the bubble residence time is higher than that with constant property case. This is due to the reduction in the bubble acceleration due to property variation with local pressure. Similarly, Fig. 19b the transient variation of channel heat transfer coefficient. The heat transfer coefficient obtained with variable property is lower than that with constant property, due to the changes in initial film thickness and the velocities of accelerated liquid slugs. Fig. 20a to 20d shows the transient variation of local heat transfer coefficient at different axial locations. The difference between the constant property case and variable property case, especially the peak values, increases along the flow direction, due to the increase in the difference between the accelerations (which influence thin film depletion rate and liquid slug velocity) obtained with constant property case and variable property case.

4.6 Effect of fluid properties

Fig. 21 shows the effect of fluid properties on the channel transient pressure drop and heat transfer coefficient for flow boiling of water ( $P_{e}$ = 101 and 19.8 kPa) and R134a ( $P_{e}$ = 76.6 kPa). For R134a, the expressions for nucleation frequency and initial film thickness are taken from Thome et al. [24]. The pressure drop and its amplitude for R134a is almost negligible as compared to that for water, due to the higher value of the product of vapour density and latent heat of vaporization for R134a. Hence, the assumption of constant property for R134a can be justified, but the same assumption cannot be used for water. The amplitude of pressure fluctuation for water with $P_{e}$ = 101 kPa is lower than that for water with $P_{e}$ = 19.8 kPa , as the product of vapour density and latent heat of vaporization is higher for water with $P_{e}$ = 101 kPa. However, the heat transfer coefficient for R134a is comparable, though lower than that for water. The heat transfer coefficient for water with $P_{e}$ = 76.6 kPa is higher than that for water with $P_{e}$ = 101 kPa due to the higher nucleation frequency and acceleration for the former.

5 Conclusions

In the present study, a 1-dimensional semi-mechanistic model combining pressure and heat transfer coefficient is proposed for flow boiling in a rectangular mini/micro-channel. The major contributions of the present study, based on the available literature, are the evaluation of transient local heat transfer coefficient in conjunction with the local transient pressure, inclusion of the effect of shear stress at liquid-vapour interface on heat transfer coefficient, addition of partial-confined bubble zone to heat transfer model and applying the developed model to evaluate the heat transfer characteristics under flow reversal condition, for flow boiling in a rectangular microchannel. The proposed model is validated with the transient and time-averaged data available in the literature. The major conclusions that can be drawn from the present work are as follows.

For high aspect ratio ( $w>>h$ ) and shorter channels, the duration over which the partial confinement occurs can be significant and therefore the present model incorporates the same into the heat transfer model. The local heat transfer coefficient is evaluated considering the heat transfer through the thin film and the surrounding bulk liquid. 2. 2.

The proposed model takes in account the evaporation as well as the shear stress at liquid-vapour interface for the evaluation of rate of thin film depletion during flow boiling in a rectangular mini/micro-channel. Results indicate that the inclusion of shear stress increases the heat transfer coefficient and also leads to early local dryout. 3. 3.

The present model has been applied to determine the heat transfer characteristics under flow reversal condition caused by inlet compressibility. Results show that the heat transfer coefficient slightly increases due to the larger contribution from the thin film evaporation resulting from the longer residence period of the elongated bubble. The model also demonstrates that the flow reversal can lead to partial / full dryout based on the operating conditions. 4. 4.

The proposed model incorporates the variation of fluid properties with local pressure. The estimated transient and time-averaged pressure drop and heat transfer coefficients are lower than those evaluated with constant pressure equal to the channel exit pressure. The effect of variable vapour properties is significant for water and is found to be negligible for refrigerants. 5. 5.

The present model for a rectangular channel neglects the transient conduction in the thin film. This can be a future study. The present model does not impose thermodynamic equilibrium, hence the estimated exit qualities are lower than the thermodynamic qualities, which seem to match with the CFD prediction in the literature. The presence of superheated liquid slugs, responsible for lowering the quality during flow boiling in mini/microchannels, needs further investigation.

Nomenclature

[TABLE]

Greek letters

[TABLE]

Subscripts

[TABLE]

Abbreviations

[TABLE]

Bibliography63

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Harirchian and Garimella [2010] T. Harirchian, S. V. Garimella, A comprehensive flow regime map for microchannel flow boiling with quantitative transition criteria, International Journal of Heat and Mass Transfer 53 (2010) 2694–2702. URL: http://dx.doi.org/10.1016/j.ijheatmasstransfer.2010.02.039 . doi: 10.1016/j.ijheatmasstransfer.2010.02.039 . · doi ↗
2Mirmanto [2016] M. Mirmanto, Local pressure measurements and heat transfer coefficients of flow boiling in a rectangular microchannel, Heat and Mass Transfer/Waerme- und Stoffuebertragung 52 (2016) 73–83. doi: 10.1007/s 00231-015-1671-7 . · doi ↗
3Yin and Jia [2016] L. Yin, L. Jia, Confined bubble growth and heat transfer characteristics during flow boiling in microchannel, International Journal of Heat and Mass Transfer 98 (2016) 114–123. URL: http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.02.063 . doi: 10.1016/j.ijheatmasstransfer.2016.02.063 . · doi ↗
4Sobierska et al. [2007] E. Sobierska, R. Kulenovic, R. Mertz, Heat transfer mechanism and flow pattern during flow boiling of water in a vertical narrow channel experimental results, International Journal of Thermal Sciences 46 (2007) 1172–1181.
5Jagirdar and Lee [2016] M. Jagirdar, P. S. Lee, Study of transient heat transfer and synchronized flow visualizations during sub-cooled flow boiling in a small aspect ratio microchannel, International Journal of Multiphase Flow 83 (2016) 254–266. doi: 10.1016/j.ijmultiphaseflow.2016.04.004 . · doi ↗
6Bigham and Moghaddam [2015] S. Bigham, S. Moghaddam, Microscale study of mechanisms of heat transfer during flow boiling in a microchannel, International Journal of Heat and Mass Transfer 88 (2015) 111–121. URL: http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.04.034 . doi: 10.1016/j.ijheatmasstransfer.2015.04.034 . · doi ↗
7Jafari and Okutucu-Özyurt [2015] R. Jafari, T. Okutucu-Özyurt, Phase-field modeling of vapor bubble growth in a microchannel, Journal of Computational Multiphase Flows 7 (2015) 143–158.
8Sun et al. [2018] Y. Sun, C. Guo, Y. Jiang, T. Wang, L. Zhang, Transient film thickness and microscale heat transfer during flow boiling in microchannels, International Journal of Heat and Mass Transfer 116 (2018) 458–470.