Numerical Analysis of Two Phase Flow Boiling Heat Transfer through Micro-Channel

Numerical Analysis of Two Phase Flow Boiling Heat Transfer through Micro-Channel

Dalvir Singh

Dept. of Mechanical Engineering, Meerut Institute of Engineering Technology, Meerut India

Santosh Kumar Rai

Dept. of Mechanical Engineering, Meerut Institute of Engineering Technology, Meerut India

Sanjeev Kumar Shukla Dept. of Mechanical Engineering, Meerut Institute of Engineering

Technology, Meerut India

Abstract – Day by day rapid miniaturization of electronics devices in electronic industry have resulted in continuous raise of heat removal requirements for making the more efficient the electronic devices, thereby making a call for highly thermal efficient cooling technology. Two phase flow boiling heat transfer is the one of the potential solution for heat dissipation for electronic devices. For highly efficient performance and reliability of micro-electronic and mechanical systems (MEMS) devices, the working temperate should be keep below some specific temperature less than 90 °C. The purpose of the work is to develop a numerical model to predict the behaviour of two phase forced convection in micro-channels using de-ionize water (DIW) as the coolant with uniform heat flux in laminar and turbulent flow regime. A 1-D thermal-hydraulic – homogeneous equilibrium model (TH-HEM), consisting of mass, momentum and energy conservation equations and empirical correlation for laminar and turbulent flow, and properties of substance equation as dryness fraction relation with liquid and vapour density, are solved numerically using a finite-difference upwind scheme for single as well as two phase flow and algebraic equations. Finally, the TH-HEM model is used to evaluate the location of boiling front, its variation with uniform input power, Variation in dryness fraction in two phase region, pressure distribution of fluid along the axial length, wall and fluid temperature distribution along the axial length, for one fixed inlet temperature and fixed outlet pressure condition.

Keywords Two phase, Boiling, force convectio,; micro- channel, thermal hydraulic homogenoues equlibrium model

INTRODUCTION

With the development of micro fabrication technology, micro- fluidic systems have been increasingly used in different scientific disciplines such as biotechnology, physical and chemical sciences, electronic technologies, sensing technologies etc. Micro-channels are one of the essential geometry for micro-fluidic systems; therefore, the importances of convective transport phenomena in micro-channels and micro-channel structures have increased drastically. In recent years, a number of researchers have reported the heat transfer and pressure drop data for laminar and turbulent liquid or gas flow in micro-channels.

Flow boiling heat transfer in micro-channels is now the most popular topic in heat transfer, as the need has arisen for very high heat flux cooling for the new generation of computer chips. Understandings of macroscopic flow and heat transfer have reached a mature stage, but when it comes to micro-channel flow, flow becomes notably different and complex. The universal use of MEMS devices, micro-heat

exchangers, micro-fluidics, other biomedical applications like micro drug delivery, have opened a new field for research.

The quick development of excessive power electronics devices with the expanded scaling down of micro-electronic devices and thermal issues and high speed processing are influencing general system performance and electronic packaging. In the previous couple of decade, the requirements of uniform heat dissipation from electronics devices have increased drastically from 10 W/cm2 to 1000 W/cm2 which exceeds the system capability of present cooling technology "[1]". Tuckerman and Pease expressed that single phase convective cooling in micro-channels ought to be attainable for input power 1000 W/cm or more than 1000 W/cm, and exhibited a micro-channel heat sink that dissipate 790 W/cm with 71°C temperature increment at 600 mL/min mass flow rate "[2]".

The increase in heat flux, the temperature of MEMS system and electronic devices is also increase but the performance of the devices is decreased due to thermal issue. So improvement of performance of MEMS system and electronic devices are required to reduce thermal issue. So increase the performance of MEMS system and electronic devices, the maximum heat flux should remove from the device or reduce the thermal issue. That why, flow boiling heat transfer method is used to extraction of heat flux from the device, because in flow boiling due to extraction of latent heat devices can dissipate heat without rise in temperature. Also, convective heat transfer coefficient is 0.1 MW/m for laminar flow of water in a 250m square channel, whereas flow boiling heat transfer coefficient can exceed 1 MW/m2-°C, which makes more attractive for cooling of MEMS systems and electronic devices. In this present study, TH- HEM model has been adopted for flow boiling heat transfer through micro-channel. Numerical simulation is performed with a micro-channel dimension of

500 µm hydraulic diameter and having length of 3 cm. Numerical simulation is carried out to analyze of the two phase flow boiling through micro channels using DI water as the working fluid medium to predict the location of boiling front, pressure drop, variation of fluid coolant temperature and wall temperature along the length of micro channel.

With the increase in heat flux dissipation, device temperature also increases, limiting the heat dissipation capability [3-5]. But flow boiling due to extraction of latent heat component can dissipate heat without increase in temperature. Also, single

phase heat transfer coefficient for laminar flow of water in a 200 m square channel is around 0.1 MW/m2-°C, whereas flow boiling heat transfer coefficient can exceed 1 MW/m2-°C making two phase flow more attractive for cooling of electronic devices [6-10]. In the present study, TH-HEM model has been adopted for two phase flow through micro channel. Numerical analysis has been carried out on micro channels of having diameter 250 m and length 3 cm. Results of TH-HEM studies on two phase flow through micro channels using water as the fluid medium were analyzed to predict the location of boiling front, pressure drop, variation of fluid coolant temperature and wall temperature along the length of micro channel

Problem Definition

In two phases boiling heat transfer, Thermal-hydraulic Homogeneous Equilibrium Model (TH-HEM) and thermal- hydraulic Annular Flow Model (TH-AFM) are mostly adopted [5]. The mixture remains in equilibrium in TH-HEM while void fractions occur in TH-AFM. The TH-HEM works on the principle that the mixture (two phase) acts like pseudo single phase fluid whose mass/weight is similar to that of the mass/weight of vapour and liquid contents. At this stage only latent heat exchange within the phases takes place. Due to pressure change, variation of property occurs in the micro channel, resulting in complicated terms. Properties like density, viscosity and enthalpy behave dynamically and hence quality becomes important.

Numerical simulations were carried out with a micro-channel diameter of 250µm having length of 3 cm. DI Water was taken as coolant that enters the channel at ambient temperature of 25ºC. Constant heat flux was applied to the surface. Fig 1 shows schematic diagram of micro-channel with boundary condition on witch simulation was performed. DI Water extracts heat continuously, resulting in phase transition leading to a two phase predominant flow.

Figure 1: schematic iagram of micro channel

Mathematical Modeling of Th-Hem Assumption

Applying Drichlet boundary temperatur condition at inet and system pressure at oultlet of micro-channel shown in Fig 1, pressures and temperatur at discreet point along axial length of microchannel have been . All the other properties namely enthalpy of fluid, wall temperature, enthalpy, density, dryness frction were obtained by state of equation.

Detemination of Boiling Front

The equ. 1 is essentially an energy conservation equation. Here i1 is input heat flux supplied to the micro-channel, L is

length of micro-channel, while is mass flow rate of DI water. With pressure distribution known along the axial length

of the micro-channel, enthalpy of saturated liquid was evaluated using saturation properties. The location at which bulk fluid enthalpy obtained from equ. 1 and that of saturation enthalpy become equal and intersect to each other, give the location of boiling front known as Xc.

It is very important to obtain the location of boiling front. Beyond this point vapour quality is greater than zero hence two phase flow exist while before this point, flow is in single phase. Location of boiling front can be extracted by

if(x) = i1 +

W

Calculation of Dryness fraction( Vapour Quality)

(1)

characteristics of saturation enthalpy and bulk fluid enthalpy as shown in Fig 2. Effectiveness of cooling is obtained in terms of temperatures of wall and fluid. For the same two set of equation have been solved for two distinct zones.

Drynes fraction is the ratio of vapour mass to the total mass. It

is zero in sub cooled liquid, and varies in two phase from zero to 1 and more than 1 . Hence it was evaluated in the two phase region by equ. 2. Equation 2 establishes relationship between enthalpy of liquid, vapour contained, and mixture. Before boiling front, dryness fraction was zero but after boiling front, that is two phase zone the quality of liquid, Vapour quality (Drynes fraction) was determine by using equ. 2. Where if is bulk fluid enthalpy, is the saturation enthalpy at liquied zone and is the saturation enthalpy at vapour zone at given system pressure

() = (1 ())() + () () (2)

Determation of Pressure for Single Phase Zone

Pressure distribution in single phase zone was obtained by using the pressure equation 3. Where fsp is single phase friction

factor, G is mass flux in Kg/s-m2 and is DI water density in Kg/m3. fsp was be calculated using laminar flow correlation for flow inside a pipe fsp=64/Re for this. In calculations effect of pressure drop due to contraction has been ignored as cross section area was constant and smooth.

fsp G2 Xc

T = T + Q

w f hconv(wL + 2dL)

is fin efficiency and is heat transfer coefficient. if was computed using equ. 10.

dif

(9)

Psp =

2 Dh

(3)

W dx hconvPc(Tw Tf) = 0 (10)

,

Where G, Re and fsp are kg/s-m2 , kg/m3 Reynolds number, friction factor in single phase respectively.

Pressure Equation for Two Phase Region

The equ 4 is mainly pressure drop equation derived from Navier strokes momentum conservation equation. Pressure distribution was calculated by using equ. 4. Here in the equation the right hand side first term is known as signifies friction drop while second term signifies convective acceleration term. Here effect of gravity body force has been

Here W is mass flow rate, Pc is perimeter of cross section.

Heat Transfer Coefficient Calculation

Convective heat transfer coefficient ( hconv ) for single phase was computed using correlation for laminar flow through a tube for constant heat flux condition. In two phase region, heat transfer coefficient hl was computed using correlation of Kandlikar[15] for homogeneous flow through a given tube.

k

disregarded because micro-channel is horizontal. Also capillary effect i.e. surface tension effects have been not taking

hconv = 4.36 Dh

h

(11)

to simplify the problem. Hydraulic dryness Dh was used to make this case analogous to pipe flow. All the properties:

TP = 1Co2 + 3Bo4 (12)

hl

density, viscosity, enthalpy, friction factor of two phase (fTP) was dynamic so average properties were used.

1 x 0.8

C0 = ( x )

(p v 0.5

)

pl

is convention Number

dp fTP G2

dx = 2 D

d G2

+ dx ( ) (4)

1- 4 is constant dependent upon fluid properties, hl was found using Detous Boelter correlation [12-14].

h

Average density was computed by equ. 5.

= 0.023 0.8 0.4

(13)

1 = X +

m v

1 l

(5)

Numerical Simulation Algorithml

GDEs are used to performed simultaneous equation to be

Similarly average viscosity m and two phase friction factor fTP was computed using equ. 6 and 7 respectively. Here l and v are saturated liquid and saturated vapour density at specified pressure.

solved for pressure and temperature together for all node points. Numerical algorithm developed for solving these equations is given in the following sub section. Figure 2 is showing the schematic diagram of discrete physical domain with inlet and outlet boundary conditions.

m m(

= X v

v

+ (1 X) l l

) (6)

Algorithm

#.Find out if(j) at all grid points using given temperature at first node and equ. 1 by marching in forward direction.

f = 64 = 64 m

(7)

tp Re

G Dh

#.Assume pressure P(j) at all grid points maintaining negative

In equ (6) m and l are saturated vapour and liquid viscosity at given pressure.

Energy Equation use to calculat temperatur

Energy equation was used to calculate the bulk fluid temperature and wall temperatures. In single phase as well as two phase zone fluid temperature was computed using equ. 8

pressure gradient based on outflow Drichlet pressure BCs.

#.Evaluate saturation properties at all discretized points for assumed pressure P(j) like Tfsat(j), if(j), µl(j), µv(j), l (j), v(j), Prl(j), Rel(j).

#.Use if and if saturated to find out location of boiling front. At the location of boiling front, these values should be equal.

#.Find X(j) in two phase region using and properties of and using equ. 2.

T (x) = T

+ Q x

(8)

#.Find m(j) and µm(j) using equ. 5 and equ. 6 respectively for

f 1 W Cpl L

Here Cpl is saturated specific heat capacity at specified pressure. In two phase region temperature of fluid (Tf ) is equal to saturated temperature at specified pressure. Also, temperature of water (Tw ) was calculated using equ. 9.

two phase region.

#.Similarly find ftp and fsp using equ. 7.

#.Solve pressure equation for single phase and two phase region using discretized equ. 3 and equ. 4 respectively. Iterate the process till pressure P (i) gets converged.

#.Find Tf(j) for single phase region by using equ. 8. In two phase region it will be equal to Tf saturated.

#.Find hconv(j) using equa. 11 and equ. 12 for single phase and

two phase region respectively.

#.Find Tw(j) using equ. 9.

#.Find if using equ. 10. Iterate if till convergence is achieved.

Ti EXIT

j=1 j=2 Xc j=M

Figure 2: Discretization of physical domain.

RESULTS AND DISCUSSION

Numerical simulation was carried out for fixed inlet temperature and fixed outlet system pressure condition. The present work are presented in the form of location of boiling front, Its variation with heat fluid supplied Q, Variation in dryness fraction X in single and two phase region, pressure distribution of fluid along axial length, wall and fluid temperature distribution along the axial length.

Boiling Front Location

In the Fig 3 red curve shows bulk enthalpy charecteristics hile blue shows saturated enthalpy charecteristics at saturated pressure.The point of intersection gives boiling front location. The saturatin enthalpy first contant and then decreace minimal rate but bulk entholpy contounsly incerease hence saturated enthalpy and bulk enthalpy intersect to each other at a point. This point is known as boiling front location.

Figure 3: Boiling Front Location.

Variation of Boiling Front with Input power

With the increased input power supplied to the wall of micro- channel, input power supplied boiling will start early. Fig 4 shows intensely the shift of boiling front location towards inlet

section with the increase in input power supplied. With 0.5 W input power supplied it was at 8 mm while with 2 W input power supplies it was at 3 mm boiling front location. The rate of shift however would decrease and becomes almost constant; boiling front took almost fixed position beyond 2.5 W input powers supplied.

Figure 4: Variation of Boiling front with respect to input power supplied.

Vapour Quality along the Length of Channel

In single phase region dryness fraction is zero and only sub cooled liquid single face exists. In two phase region beyond Xc

, X increases almost linearly from 0 to 3 in downstream

locations as fluid will get more and more heat while moving in downstream direction. Fig 5 depicts almost liner variation of vapour quality X with channel length. It can be seen clearly that throughout two phase region faster heat transfer takes place as X has not crossed value of 1.0.

Figure 5: Variation in vapour quality in downstream direction points.

Pressure Variation along the Length

The variation of system pressure in single phase region is very low while that in two phase region is very high. The reason being in single phase region pressure drop is only due to friction while in later acceleration terms also contribute to pressure drop. Also friction loss is accompanied with a multiplier point whose value is greater. Pressure drop due to acceleration term dominates in two phase flow Fig 6.

Figure 6: Variation in Pressure distribution in downstream direction Points.

Pressure Drop Variation with Input power

Pressure drop (dP) increases with increase in input power because with increase in input power there is increase in convection current (more turbulence) which in turn increases pressure drop. Moreover as input power increases, two phase region would be larger and hence more pressure drop take place Fig 7.

Figure 7: Variation in pressure drop with respect to power input.

Bulk Temperature and Wall Temperature Variation Along the Length

Figure 8 shows in single phase region both bulk temperature (Tf) and wall temperature (Tw) increases because of increased fluid enthalpy as a result of heat input. In two phase region beyond Xc due to existence of saturation condition and to maintain favorable pressure gradient both Tf and Tw decreases.

Figure 8: Variation of Tw and Tf along micro-channel length

CONCLUSION

A numerical analysis was conceded out to estimate the existence of location of boiling front, pressure drop, and pressure drop variation with input power, variation of boiling front with input power, bulk temperature (Tf) and wall temperature (Tw) and variation of quality of fluid along the length of micro-channel. The present results of TH-HEM studies on two phase flow through micro- channels using DI water as the fluid medium were analyzed. GDEs analysis incorporates effects of fluid flow rate, heat transfer and power input. The GDEs were solved to predict the location of boiling front, pressure drop and variation of quality of fluid along the length of micro-channel.

• Increasing in power input location of boiling front decreases

• Increasing in power input pressure drop increases linearly.

• Variation of pressure along length first decrease with minimal rate in single phase and after reaching in two phase zone pressure is drastically decreases

• Bulk temperature (Tf) and wall temperature (Tw) increases up to single phase after reaching in two phase zone the bulk temperature and wall temperature decrease and both the temperature are almost equal

REFERENCES

1. Mudawar and Issam, 2001, Assessment of High-Heat Flux Thermal Management Schemes", Components and Packaging Technologies, IEEE Transactions no. 2, pp. 122-141.

2. Tuckerman, David, B. and Pease, R. F. W., 1981, "High Performance Heat Sinking for VLSI", Electron Device Letters, IEEE 2, no. 5, pp. 126-129.

3. Kandlikar, S. G., 2002, "Fundamental Issues Related to Flow Boiling in Minichannels and Microchannels", Experimental Thermal and Fluid Science 26, no. 2, pp. 389-407.

4. Thome and John, R., 2004, "Boiling in Microchannels: A Review of Experiment and Theory", International Journal of Heat and Fluid Flow 25, no. 2, pp. 128-139.

5. Kandlikar, S.G., 2004, Heat Transfer Mechanisms During Flow Boiling in Microchannels, Journal of Heat Transfer, ASME,

   Vol.126, pp. 8-16

6. Proceedings of the 4th International Conference on Nanochannels, Microchannels, and Minichannels, ASME, Ashman, Sean, and Kandlikar, S. G., 2006, "A Review of Manufacturing Processes for Microchannel Heat Exchanger Fabrication", pp. ASME, pp. 855-860.

7. Morgan, Chris, J, Vallance, R.R. and Marsh, E.R., 2006, "Micromachining and Microgrinding with Tools Fabricated By Micro Electro Discharge Machining", International Journal of Nano Manufacturing 1, No. 2, pp. 242-258.

8. Ghiaasiaan, S.M. and Abdel, S.I., 2001, Two Phase Flow in Microchannels, Advances in Heat Transfer 34, pp. 377-394.

9. Proceedings of the 2nd UK National Heat transfer Conference, 1988, Damianides, C.A. and Westwater, J.W., Two Phase Flow Patterns in a Compact Heat Exchanger And in Small Tubes, pp. 12571268.

10. Fukano, T. and Kariyasaki, A., 1993, Characteristics of Gas Liquid Two Phase Flow in Capillary Nuclear Engineering and Design 141, pp. 5968.

11. Vassighi, Arman, Sachdev, Manoj, 2006, Thermal and Power Management of Integrated Circuits. Integrated Circuits and Systems

12. Sarangi, R. K., Bhattacharya, A., and Prasher, R. S., 2009, "Numerical Modelling of Boiling Heat Transfer in Microchannels", Applied Thermal Engineering 29, no. 2 , pp. 300-309.

13. Todreas, N. E., and Kazimi, M. S., 1 st ed., 1990 Nuclear Systems: Thermal and Hydrulic Fundamentals, Taylors and Francis Publishers, Vol. 1., pp. 488

14. Incropera, F. P. and Dewitt, D. P. 2002, Fundamentals of Heat and Mass Transfer 6th ed. Wiley, New York, Chap.3, pp.149

15. Kandlikar, S.G.1990, "A General Correlation For Saturated Two Phase Flow Boiling Heat Transfer Inside Horizontal And Vertical Tubes", Journal of Heat Transfer 112, no. 1, pp. 219- 228.

Nomenclature

Tf bulk temperature

Tw wall temperature

Q heat fluid supplied

i1 heat supplied

L length of micro-channel

W mass flow rate

ii saturation enthalpies of liquid iv saturation enthalpies of vapour if bulk fluid enthalpy

fsp single phase friction factor G mass flux in Kg/s-m2

DIwater density in Kg/m3

Re Reynolds number

Dh Hydrolic dryness

Ftp friction factor of two phase

m Average density

m average viscosity

l saturated liquid density at specified pressure.

v saturated vapour density at specified pressure.

Cpl saturated specific heat capacity at specified pressure. fin efficiency

Pc perimeter of cross section