- Email: [email protected]
Evaporative heat transfer analysis of a micro loop heat pipe with rectangular grooves
International Journal of Thermal Sciences 79 (2014) 51e59
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Evaporative heat transfer analysis of a micro loop heat pipe with rectangular grooves M. Ghajar a, J. Darabi b, * a b
Techlegic LLC, USA Department of Mechanical Engineering, Southern Illinois University Edwardsville, Edwardsville, IL 62026, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 July 2013 Received in revised form 16 December 2013 Accepted 19 December 2013 Available online 1 February 2014
This study focuses on the thermal and capillary analysis of a micro loop heat pipe for the thermal management of electronic devices and systems. A model is developed using the principles of thin film evaporation to predict the evaporative heat transfer coefficient in grooved capillary structures. In addition, a micro-flow submodel is developed to compute the dry-out distance in rectangular capillary grooves. These submodels are incorporated into our previously developed system-level loop solver model to simulate the performance of a flat micro loop heat pipe. The integrated model predicts the thermal performance, evaporator surface temperature, and local and average heat transfer coefficients as a function of the applied heat load. The modeling results are verified by comparison with the experimental data for a similar device and a good agreement is obtained. Ó 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Micro loop heat pipe Evaporative heat transfer analysis Rectangular grooves Numerical model
1. Introduction Thermal systems miniaturization and the consequent rapid increase in power density of advanced microprocessors and electronic components have created a demand for achieving high heat dissipation rates. In most cases, such high heat fluxes cannot be easily dissipated using existing cooling techniques. Thus, research for new cooling technologies is critically important for space and defense applications. The quality of thermal management in electronic packaging directly affects the performance, cost, and reliability of electronic devices. A micro loop heat pipe (MLHP) can be a good candidate for such applications. Because no energy is needed to operate a loop heat pipe (LHP), it has been identified as an enabling technology for spacecraft thermal control. Conventional loop heat pipes have already been employed successfully in many applications, including spacecraft thermal control, electronics cooling, and ordinary commercial thermal devices. A loop heat pipe is a two-phase device, derived from heat pipe, which uses capillary forces developed inside its wicked evaporator to pump and circulate its working fluid in a closed loop, allowing heat transport over large distances at nearly constant temperatures. Major advantages of heat pipe include a very high thermal conductance, no pumping power requirements, no moving parts,
* Corresponding author. E-mail address: [email protected] (J. Darabi). 1290-0729/$ e see front matter Ó 2014 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ijthermalsci.2013.12.014
and relatively low pressure drops. More details on the principle of the heat pipes and loop heat pipes can be found elsewhere [1e7]. Loop Heat Pipes (LHPs) were invented in the former Soviet Union in the early 1980s [8,9]. LHPs have been used as baseline thermal control systems for a number of spacecrafts and commercial satellites due to their design simplicity and operation robustness [10e 15]. In addition to space applications, the high pumping head also makes LHP feasible in many terrestrial applications [12]. A comprehensive review of LHPs has been included in a previous publication of the authors [16] and other researchers [11e15]. A large number of analytical and experimental investigations have been conducted by various researchers to predict the thermal and fluid flow in conventional LHPs. However, a limited number of studies have focused on micro/miniature loop heat pipes [17e21]. Advances in microfabrication techniques have enabled the development of miniaturized, low-cost, and integrated MLHPs. A MLHP can be directly attached to the surface of an electronic component for which localized cooling is required. In an earlier work, the authors developed a FORTRAN-based system-level numerical model to simulate a planar MLHP and predict its thermal behavior as a function of the applied heat load [16]. Later, an improved model was developed to predict the evaporator surface temperature by assuming a typical evaporative heat transfer coefficient [22]. However, the evaporative heat transfer coefficient depends on several operational and geometrical parameters such as heat flux, evaporator temperature, fluid properties, and the geometry of the evaporator and the grooves. The
52
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
Nomenclature
l
A AR Dg h k K _ m Pvlv Pvlvt Q Qappl q_ q00 R RT T U V
n s
W x
dispersion coefficient, N m aspect ratio of the grooves distance between two grooves, m depth of the groves, m thermal conductivity, W/m K curvature, m1 mass flow rate, kg/s vapor pressure of the thin film, Pa saturation pressure at the film temperature, Pa heat transferred per unit length of the groove, W/m heat load, W heat flow rate, W/s heat flux, W/m2 ideal gas constant, N m/K thermal resistance, m2 K/W temperature, K fluid velocity, m/s phase velocity or velocity of fluid in phase change normal to the interface, m/s width of the groove, m coordinate along the groove, m
Greek symbols d film thickness, m q meniscus contact angle, rad
present model has been built upon prior work of the authors [23,24] to simulate the MLHP under various operating conditions and geometrical parameters. In this work, a submodel is developed to calculate the evaporative heat transfer coefficient using the principles of thin-film evaporation and incorporated into our previously developed systemelevel loop solver to further improve the simulation tool. Parameters of interest such as the local heat transfer coefficient, micro region characteristics, and dry-out distance were calculated iteratively using a procedure developed for the main loop solver [22]. 2. Model formulation Fig. 1 depicts a schematic illustration of a planar MLHP. The loop consists of an evaporator, a condenser, a compensation chamber, and liquid and vapor lines. A sectional view of the compensation chamber, liquid line, and evaporator, which are connected by an
coordinate perpendicular to the groove used in the micro-region, m kinematic viscosity, m2/s surface tension, N/m
Subscripts 0 initial appl applied ave average c capillary cc compensation chamber d disjoining e evaporator g grooves i interfacial iv interface (vapor side) l liquid ll liquid line lv liquidevapor n normal s surface sat saturation tot total T total v vapor w wall, wick, width
array of parallel microgrooves is shown in Fig. 1b. It is widely accepted that the temperature difference between the two ends of a wicking structure in a cylindrical evaporator can be estimated using the following equation [25,26]:
Te Tll ¼
dT dP
sat
DPtot DPg
(1)
where Te is the temperature of the capillary grooves adjacent to the evaporator, and Tll is the temperature of the capillary grooves adjacent to the liquid line, DPtot is the total pressure drop of the loop, and DPg represents the pressure drop of the wicking structure. Similarly, Eq. (1) can be used to estimate the temperature difference between the two ends of the wicking structure in a MLHP. However, for a MHLP, the presence of a conductive substrate can result in significantly higher heat conduction to the liquid line under a more complex transport pattern. The authors addressed this issue by
Fig. 1. (a) Schematic diagram of a planar MLHP; (b) sectional view (AeA)- the compensation chamber, liquid line, and evaporator are connected by microgrooves (not to scale).
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
estimating heat conduction to the compensation chamber using a finite difference method [16,22]. In this study, the principals of thin-film evaporation are incorporated into a micro flow solver to determine the wetted length of the groove during evaporation. Fig. 2 shows a schematic illustration of the liquid profile along a rectangular groove in the evaporator. As the liquid flows downstream towards the vapor line, the crosssection of the liquid film in the grooves deceases due to evaporation. At any given cross section, the evaporative heat transfer coefficient at the interlines is the main contributor to the local heat transfer coefficient. Thus, the local heat transfer coefficient increases as the number of interline regions increase. The interline regions of the thin-film are indicated by circles in Fig. 2b. An average normal evaporative heat transfer coefficient is defined as [27]
hn ¼
q00 ave Ts Tsat
(2)
where q00 ave is the average applied heat flux and Ts is the temperature at the inner surface of the evaporator. The vapor pressure in the groove is assumed to be constant and the heat transfer from the dry-out regions is assumed to be negligible. To calculate the substrate surface temperature, the local heat transfer coefficient must be known. The thin film evaporation theory can be used to calculate this parameter. The relationship between the vapor pressure of the thin film and the saturation pressure is given by [28]
ln
Pvlv Vl ¼ P Pvlvt RTSat d
(3)
where Pvlv and Pvlvt are the vapor pressure of the thin film and saturation pressure at the film temperature, and Pd is the disjoining pressure. For polar liquids, the disjoining pressure is given by [29]
Pd ¼
d ln 3
A
d
d0
(4)
where A is the dispersion coefficient and d is the liquid film thickness. The heat transfer through the liquid film is primarily dominated by conduction since the film thickness is very small. The interfacial resistance to heat transfer has been previously studied by several researchers [30e33] and is not repeated here for brevity. The following fourth order differential equation has been developed for the film thickness [27]. This equation requires four initial conditions for d. For instance, the value of d along with its three derivatives should be known at x ¼ 0.
hfg d Tw Tiv dP d3 c ¼ RT 3nl dl dl
53
(5)
Obtaining the required initial conditions for d to solve equation (5) is a challenging task. A recent approach by Sartre et al. [34] attempts to solve an alternative form of equation (5) in which the following system of differential equations are solved simultaneously:
dd 0 ¼ d dl dd ¼ dl
(6)
1þd
02 1:5
s
A Pc 3
d
(7)
dPc 3ml Q ¼ 3 dl hfg rl d
(8)
dQ ¼ ðTw Tiv ÞR1 T dl
(9)
The advantage of this approach is that the initial conditions at x ¼ 0 can be more intuitively understood as follows:
d ¼ d0 ;
K ¼ 0;
and
d0 ¼ 0
(10)
The radius of film curvature is constant in the intrinsic meniscus region and outside the microregion. This fact provides the remaining required boundary condition. The wall superheat temperature and the evaporative heat transfer coefficient are locally interdependent. Therefore, an iterative procedure was used to address this issue in which the values of the evaporative heat transfer coefficient and the wall temperature are updated alternatively. The evaporative heat transfer coefficient is calculated by the micro-region model while the wall temperature is determined by solving a steady-state heat conduction equation.
V2 T ¼ 0
(11)
The temperature variation in the solid domain is expected to be small for this microscale system. The alternate direction implicit (ADI) scheme was used in the finite difference scheme. All variables in the computer program were chosen as double precision. The mesh was refined to ensure convergence for a range of applied heats. Furthermore, analytical calculations and commercial software were used to verify the accuracy of the model for several cases. The system of equations (6)e(9) was solved numerically using a fourth order Runge-Kutta method. The meniscus profile was
Fig. 2. Thin film profile along the grooves in the evaporator.
54
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
constructed using the calculated meniscus radius and contact angle obtained from the microregion code. The local heat transfer coefficient of the thin film outside the microregion was calculated based on a 1-D conduction through the liquid meniscus. The localized temperature and evaporative heat flux information was obtained from the converged solution for the micro/macro regions. The local evaporative heat transfer coefficient was then calculated from:
hðlÞ ¼
q00 ðlÞ ; Tw ðlÞ Tiv ðlÞ
(12)
where q00 ðlÞ is the local heat flux and Tw(l) and Tiv(l) are the local wall temperature and local interface temperature in the evaporator. A cross-sectionally averaged heat transfer coefficient, h(x), was then obtained by integrating local heat transfer coefficient over l and was used in the loop solver. In this study, the triangular region in deep grooves was not included because of its short length compared to the rectangular region. For shallower grooves, it is necessary to consider both the rectangular and the triangular sections by combining the triangular (not the subject of this paper) and rectangular models. Fig. 3 shows a comparison between the triangular and rectangular geometries. Traditionally, for triangular grooves, the cross sectional areas for the liquid and vapor passage are expressed as functions of rm. This is because in a triangular geometry, the height of the liquid meniscus, Lw, can be geometrically related to the radius of curvature and contact angle of the meniscus. However, such a geometrical relation does not exist for a rectangular geometry. In a rectangular geometry, a similar analysis requires additional information about the liquid meniscus in the groove because rm is not related to the wetted height of the groove. Nilson et al. [35] proposed that the curvature of the meniscus in a rectangular micro groove only depends on the width of the groove. While this can be a reasonable assumption for wider rectangular channels, it does not necessarily result in an interesting solution for a micro groove for the following reason. It is acceptable to consider both as a fact and a precondition for the development of the above continuum-based formulations that the liquid film profile is a C1 continuous function. Therefore, it must be a smooth extension of the micro-region film thickness profile at the intrinsic meniscus junctions with the micro-region film. If the meniscus curvature is only a function of the groove width as in Refs. [35], the curvature at a point where the intrinsic meniscus starts to form, is constant for any wall temperature as long as the channel width remains constant. However, the
Fig. 4. Mass and momentum balance in a given cross-sectional element of the groove. The momentum terms due to shear stress at the interface are not shown.
contact angle which determines the slope of the microregion film profile is a function of the wall temperature. As a result, the C1 continuity condition at the junction cannot be satisfied. Thus, rm should be related to both the groove width and the contact angle. If the contact angle is available from the solution of the micro region model, the rm is known and the wetted height, Lw, can be calculated from the micro-flow model [36]. In the present work, a uniform meniscus radius is used by assuming a uniform contact angle obtained from the micro-region model in a single cross-section of the groove. This is a reasonable assumption for calculating the capillary pressure in microgrooves [3e5], allowing to determine the capillary driving force from the Young-Laplace equation. Ideally, the micro-region solver should be run for several cross-sections along the groove to calculate the contact angle and the meniscus radius. In a separate publication, the authors demonstrated that a constant contact angle assumption along the groove was valid [36]. The meniscus half-angle of curvature and the height of the meniscus can be geometrically determined using Fig. 3b:
Fig. 3. Comparison of liquid menisci formed in a triangular microgroove (a) and a rectangular microgroove (b).
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
q ¼ 2 sin1 Lw ¼
w 2rm
(13)
1 1 2 Al þ rm ðq sin qÞ w 2
(14)
The continuity and the momentum equations are derived for the one-dimensional two phase evaporating flow as shown in Fig. 4. The governing equations for this case are as follows:
d ½U A ¼ Vil Lci dx l l
(15)
d ½Uv Av ¼ Viv Lci dx
(16)
d P 1 Ul2 Al þ Al l ¼ ½slw ðw þ 2Lw Þ þ si Lci þ gAl sinð4Þ rl rl dx
d Pv U 2 Av þ Av rv dx v
¼
1
rv
½svw ðw þ 2Lw Þ si Lci
AT ¼ Al þ Av ¼ wh þ De w þ dg Pv ¼ Pl þ
s rm
rm
Lci ¼ 2rm q
where Lci is the length of the liquid surface curve at a given crosssection and 4 is the inclination angle of the grooved plate. The effect of the inclination is negligible in a micro loop heat pipe. Unlike a heat pipe, where the liquid and vapor flow in opposite directions, vapor flows concurrently with the liquid in a loop heat pipe. Therefore, the shear stress at the interface facilitates the flow of liquid due to a higher vapor velocity. The boundary conditions for the flow are similar to previous heat pipe studies [34,37,38]. The system of differential equations was solved using a Runge-Kutta scheme. The shear stresses at the wall and interface were calculated using the friction factor concept and known exposed areas at the wall and the liquidevapor interface. Equations (15)e(20) were solved to determine the liquid cross-section. The height of the meniscus at the wall was then calculated using Eq. (14). The friction factors for the grooves were calculated from Shah’s correlation [39]. The boundary conditions for liquid and vapor phase velocities are dependent of the liquid and vapor mass flow rates in the grooves [34,37,38].
rm ðrm q w cos qÞ 2
(17)
Al ðx ¼ 0Þ ¼ wh
(18)
Av ðx ¼ 0Þ ¼ AT eAl ðx ¼ 0Þ
(19)
Ul ðx ¼ 0Þ ¼
_ m whrl
(27)
Uv ðx ¼ 0Þ ¼
rl Ul Al rv Av
(28)
(20)
The unknowns are Ul, Uv, Al, Av, Pl, and Pv. The rm, Lci, Vil, and Viv are given in the following relations:
w ¼ 2 cos ðaÞ
55
Pv ¼ Pl þ s=rm
(21)
where Pl is an arbitrary base pressure.
(22)
3. Solution algorithm
Vil ¼
qf w rl hfg Lci
(23)
Viv ¼
rl V rv il
(24)
Fig. 5. Comparison of the predicted heat removal capability with the experimental data of Hsu et al. [40].
(25) (26)
(29)
The model consists of three parts. The first part which is based on the authors’ previous model [23] has remained mainly unchanged and is referred to as the “loop solver” in this paper. Briefly, the loop solver calculates the pressure drops in all segments of the loop and the temperature difference at the two ends of the wicking structure between the liquid line and the evaporator. Once, the liquid line temperature is known, a function is defined as the difference between the calculated temperatures from the available temperature field obtained from the finite difference model. By varying the saturation temperature, this function is minimized to determine the operating temperature. The second part is referred to as “the micro-region solver (MRS)” and provides an in-situ calculation of the evaporative heat transfer coefficient in the evaporator grooves. The MRS consists of two submodels. One submodel calculates the evaporative heat transfer coefficient and the other updates the wall temperature. The third part, which is referred to as the “micro-flow solver (MFS)” estimates the dry-out length by solving one-dimensional two-phase flow equations (continuity, momentum, and energy for both phases). An iterative procedure was used in a system-level loop solver model [23]. Using an initial wall temperature, the MRS determines the thin film profile and updates the heat transfer coefficient and the interfacial temperature. For a given applied heat flux, the wall temperature inside the evaporator is then calculated using a finite difference routine. This process is repeated until a converged solution is obtained for the wall temperature. After this procedure ends, the local heat transfer coefficient and wetted length are calculated from the micro-flow solver. The average heat transfer
56
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
coefficient in the grooved structures is determined by integrating Eq. (12) and used in the next iteration of the main loop solver.
Table 1 Geometrical dimensions used in the model e the base geometry. Parameter
Abbreviation
Dimension
Unit
Width of the groove Thickness of the groove wall Length of the groove Depth of the groove Length of the vapor line Length of the liquid line Width of the vapor line Width of the evaporator Length of the evaporator Depth of the liquid and vapor line, condenser, CC, and evaporator
W Dg Lg h Lv Lll Wv We Le hd
5 2 500 50 20 30 2 2 4 200
mm mm mm mm mm mm mm mm mm mm
4. Results and discussion The accuracy of the system-level model was verified by comparison with the experimental data of Hus et al. [40] using methanol as a working fluid. As shown in Fig. 5, the model predicts the experimentally measured evaporator temperature with good accuracy. Unless otherwise noted, the rest of the results in this section are presented for a micro loop heat pipe of the geometrical dimensions listed in Table 1 using methanol as the working fluid. It was demonstrated in a previous publication of the authors [22] that the device performance reached a saturation point as the heat transfer coefficient in the evaporator exceeded certain values. At very high heat transfer coefficients, the frictional forces present in the liquid and vapor lines and other components increase and have
Fig. 6. Profiles of various parameters in the micro-region using methanol as the working fluid.
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
57
Fig. 9. Comparison of the modeling results with the similar experimental data from Ref. [42]. Fig. 7. Variation of contact angle vs. DT ¼ Tw Tsat.
a negative impact on the overall device performance. Thus, the overall heat removal capability for a given operating temperature may not necessarily improve at high heat transfer coefficients. However, the evaporator wall temperature, which is the most critical parameter, is a strong function of the heat transfer coefficient inside the evaporator and can be greatly reduced with increasing the heat transfer coefficient. The results section consists of three separate parts pertaining to three parts of the model. The results of the micro-region model are presented first. In a micro groove, the groove length is several orders of magnitude larger than the length scale of its crosssection. Thus, the calculations are essentially based on the assumption of a one-dimensional heat transfer from the inner wall of the groove to the liquidevapor interface. In the second part, the micro-flow model predicts the groove behavior during evaporation using a one-dimensional hydrodynamic model. The purpose of this investigation is to ensure the grooves are capable of removing the applied heat load without forming a complete dryout region. Therefore, this sub-model is solved at each iteration of the loop solver and serves as a continuous examiner of the main model. This study investigates the behavior of the device at steady-state condition. A situation where a complete dry-out occurs in a section of the groove should be treated as a timedependent problem and is beyond the scope of this study. Finally, in the last part, the results of the integrated model which consists of the micro-region, micro-flow, and loop solver submodels are demonstrated.
Fig. 10. Variation of the wetted height along the groove at various groove depths. Please refer to Fig. 3b for definition of Lw.
4.1. Micro-region modeling results The results of the micro-region model using methanol as the working fluid are presented in Fig. 6. The general behavior of the results is similar to previously published data for ammonia [29,41]. As shown in Fig. 6a, the film thickness remains at the initial value to
Fig. 8. Variation of the heat transfer coefficient and wall temperature in a full cross-section of a 10 mm by 10 mm groove.
58
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
heat transfer to the evaporator is calculated in this part of the model and used to determine the mass flow rate. 4.2. Micro-flow modeling results
Fig. 11. Maximum tolerable heat flux as a function aspect ratio at various groove widths and a wall thickness of 5 mm.
mark the non-evaporating adsorbed layer before increasing sharply to specify the location of the interline. As shown in Fig. 6b, the distance with high evaporation heat flux, is limited to a few tenth of a micron which includes the interline. The interline region also corresponds to the region with a high evaporative heat transfer coefficient as shown in Fig. 6c. The significant part of the evaporation phenomena occurs in this region. Fig. 6d shows that the curvature reaches a maximum at the interline followed by formation of a constant curvature meniscus. The radius of this meniscus depends on the contact angle and the groove dimensions. As shown in Fig. 6e and f, the capillary pressure and the interface temperature behave quite similar since they are related linearly. The interface temperature approaches the saturation temperature in regions outside the interline. The variation of the contact angle as a function of the difference between the wall and the saturation temperature is illustrated in Fig. 7. As mentioned before, the contact angle data are used to determine the meniscus radius. The film thickness and subsequently the heat transfer coefficient outside the micro-region are determined by adopting a 1-D conduction approach and a known meniscus radius. The heat transfer coefficient and wall temperature profiles for a 10 mm by 10 mm groove are presented in Fig. 8. As expected, the minimum wall temperature occurs at locations with highest heat transfer coefficient values. The
In order to verify our model prediction and methodology, the simulation results were compared to available experimental data for an array of grooves with open triangular channels machined into a flat plate [42]. Fig. 9 shows a comparison between the modeling results and the experimental data for the dry-out distance as a function of the axial inclination angle. The results indicate that the dry-out distance is in good agreement with the experimental data and follows the same general trend. The liquid mass flow rate at a given cross-section along the groove decreases as it evaporates and enters the vapor phase. In addition, the liquid velocity decreases further downstream of the groove since the viscous forces overcome the inertia forces. Eventually, if the heat flux is high enough, a stationary condition is formed for the liquid phase at the dry-out point. Extremely high vapor velocities close to or exceeding the sonic limit are detrimental to the function of these devices and need be monitored in the design. Such a limit occurs at a heat load of approximately 6 W for the baseline dimensions listed in Table 1. The effect of the groove depth on the dry-out length is depicted in (Fig. 10). It is observed that as the groove depth increases, the dry-out occurs further away from the groove inlet. To better understand this phenomenon, the maximum heat flux that a groove can bear before a dry-out occurs as a function of the groove aspect ratio (depth/width) is shown in Fig. 11. The results indicate that for the range of the aspect ratios investigated, a groove with a width of 5 mm can deliver the highest heat flux before it experiences a dryout. However, for large aspect ratios, a narrower groove is expected to perform better since the capillary force is higher for a narrow groove due to smaller width, W (i.e. Pc w 1/W). For a constant width, a higher aspect ratio corresponds to a larger cross-sectional area. Thus, the friction forces decrease as the aspect ratio increases. This is the primary reason that the graphs approach constant values at large aspect ratios. The graphs intersect because far from the capillary limit of the device, a larger width (i.e. a larger crosssection) constitutes a lower contribution of the friction forces. 4.3. MLHP Modeling results Fig. 12 presents the cross-sectionally averaged heat transfer coefficient and maximum evaporator surface temperature as a function of the heat load. The results indicate that the crosssectionally averaged heat transfer coefficients decrease as the heat load increases. This is attributed to the change in the contact angle of the meniscus and subsequently a change in the meniscus radius and the thickness of the extended film. The evaporator surface temperature at a heat load of approximately 7 W exceeds the allowable temperature for electronics cooling applications (w100 C) for the geometrical dimensions used in this study. The heat removal capability is expected to improve by shortening or widening the vapor line in the loop. Also, our previous studies have shown that there exist optimum groove dimensions in the wicking structure. A parametric study similar to previous publications of the authors [16,22] can be performed to determine the optimal dimensions as well as the most suitable working fluids. 5. Conclusions
Fig. 12. Variation of the cross-sectionally averaged evaporative heat transfer coefficient and the maximum surface temperature of the evaporator as a function of the total heat load.
A new model was developed to simulate a micro loop heat pipe with rectangular capillary grooves. This model utilizes an in-situ calculation of the evaporative heat transfer using the principles of
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
the thin film evaporation. In addition, a micro flow solver submodel was incorporated into the model to estimate the length of the wetted evaporating region. This model is not intended to simulate situations where dry-out occurs in the evaporator. Such a case requires a transient solution rather than a steady state solution presented here. This integrated model which is built upon previously developed models by the authors is capable of predicting the operating characteristics and evaporative heat transfer coefficient of the micro loop heat pipe for any geometrical dimensions and operational parameters. As in the previous publications of the authors [16,22], the performance of the device is highly dependent on the geometrical dimensions and choice of working fluid. Thus, the optimization of the device geometry, particularly the groove dimensions are necessary to design an improved MLHP device. References [1] S.W. Chi, Heat Pipe Theory and Practice, Hemisphere Publishing Co, Washington, D.C, 1976. [2] P.D. Dunn, D.A. Reay, Heat Pipes, second ed., Pergamon Press, 1978. [3] G.P. Peterson, An Introduction to Heat Pipes: Modeling, Testing, and Applications, John Wiley & Sons, Inc., New York, 1994. [4] J.T. Dickey, G.P. Peterson, Experimental and analytical investigation of a capillary pumped loop, J. Thermophys. Heat Transfer 8 (3) (1994) 602e607. [5] A. Faghri, Heat Pipe Science and Technology, Taylor & Francis, Inc, Washington DC, 1995. [6] L.L. Vasiliev, Micro and miniature heat pipes e electronic component coolers, Appl. Therm. Eng. 28 (4) (2008) 266e273. [7] S. Launay, V. Sartre, J. Bonjour, Parametric analysis of loop heat pipe operation: a literature review, Int. J. Therm. Sci. 46 (7) (2007) 621e636. [8] Y.F. Maidanik, S.V. Vershinin, V.F. Kholodov, Y.F. Dolgirev, US Patent No. 4515209, May 7, 1985. [9] V.M. Kiseev, A.G. Belonogov, N.P. Pogorelov, Development of two-phase loops with capillary pumps, SAE 972482, July 14e17, 1997. [10] Y.F. Maidanik, Y.G. Ferchtater, K.A. Goncharov, Capillary Pump Loop for the Systems of Thermal Regulation of Spacecraft, 1991. ESA SP-324. [11] Y.F. Maidanik, Y.G. Ferchtater, Theoretical basis and classification of loop heat pipes and capillary pumped loops, in: 10th International Heat Pipe Conference, Stuttgart, Germany, 1997. [12] J. Ku, Operating characteristics of loop heat pipes, SAE Paper (1999). No. 199901-2007. [13] D.A. Wolf, D.M. Donald, A.L. Phillips, Loop heat pipes e their performance and potential, in: SAE Paper No. 941575, 24th ICES, Germany, 1994. [14] A.L. Phillips, J.E. Fale, N.J. Gernert, D.B. Sarraf, W.J. Bienert, Loop Heat Pipe Qualification for High Vibration and High-G Environments, 1998. AIAA-98e 0885. [15] A. Hoelke, H.T. Henderson, F.M. Gerner, M. Kazmierczak, Analysis of the Heat Transfer Capacity of a Micromachined Loop Heat Pipe, International Mechanical Engineering Congress and Exposition, IMECE, HTD-Vol. 364-3, ASME Publication, Nashville, Tennessee, November 15e20, 1999, 53e60. [16] M. Ghajar, J. Darabi, N. Crews Jr., A hybrid CFD-mathematical model for simulation of a MEMS loop heat pipe for electronics cooling applications, J. Micromech. Microeng. 15 (2005) 313e321. [17] M.A. Chernysheva, S.V. Vershinin, Y.F. Maydanik, Operating temperature and distribution of a working fluid in LHP, Int. J. Heat Mass Transfer 50 (2007) 2704e2713. [18] R. Singh, A. Akbarzadeh, M. Mochizuki, Operational characteristics of a miniature loop heat pipe with flat evaporator, Int. J. Therm. Sci. 47 (11) (2008) 1413e1562.
59
[19] S.F. Wang, W.B. Zhang, X.F. Zhang, J. Chen, Study on start-up characteristics of loop heat pipe under low-power, Int. J. Heat Mass Transfer 54 (4) (2010) 1002e1007. [20] J. Li, D. Wang, G.P. Peterson, Experimental studies on a high performance compact loop heat pipe with a square flat evaporator, Appl. Therm. Eng. 30 (6e7) (2010) 741e752. [21] X. Zhang, X. Li, S. Wang, Three-dimensional simulation on heat transfer in the flat evaporator of miniature loop heat pipe, Int. J. Heat Mass Transfer 54 (2012) 188e198. [22] M. Ghajar, J. Darabi, Numerical modeling of evaporator surface temperature of a micro loop heat pipe at steady-state condition, J. Micromech. Microeng. 15 (2005) 1963e1971. [23] M. Ghajar, J. Darabi, Integration of a thin-film evaporation submodel with a micro loop heat pipe solver, in: Proceedings of the IMECE2006, 2006 ASME International Mechanical Engineering Congress and Exposition, November 5e 10, 2006. Chicago, Illinois, USA. [24] M. Ghajar, J. Darabi, Simulation of steady state surface temperature response of a novel micro loop heat pipe to various heat loads, in: Proceedings of the ASME Summer Heat Transfer Conference, HT, 2005, pp. 633e640. [25] T. Kaya, T. Hoang, Mathematical Modeling of Loop Heat Pipes, 1999. AIAA Paper 99-0477. [26] T. Hoang, T. Kaya, Mathematical Modeling of Loop Heat Pipes with Two-phase Pressure Drop, 1999. AIAA Paper 99e3448. [27] P.C. Stephan, C.A. Busse, Analysis of heat transfer coefficient of grooved heat pipe evaporator walls, Int. J. Heat Mass Transfer 35 (2) (1992) 383e391. [28] B.V. Derjaguin, Z.M. Zorin, Optical study of the adsorption and surface condensation of vapors in the vicinity of saturation on a smooth surface, in: 2nd Int’l Congr. Surf. Activity (London), 2, 1957. [29] P.G. deGennes, Wetting: statics and dynamics, Rev. Mod. Phys. 57 (3) (1985). Part 1. [30] P.C. Wayner Jr., Y.K. Kao, L.V. LaCroix, The interline heat transfer coefficient of an evaporative thin film, Int. J. Heat Mass Transfer 19 (1976) 487e492. [31] B.V. Derjaguin, S.V. Nerpin, N.V. Churayev, Effect of film transfer upon evaporation of liquids from capillaries, Bull. Rilem 29 (1965). [32] J.G. Truong, P.C. Wayner Jr., Effects of capillary and van der Waals dispersion forces on equilibrium profile of a wetting liquid: theory and experiment, J. Chem. Phys. 87 (7) (1987) 4180e4188. [33] V.P. Carey, A.P. Wemhff, Disjoining pressure effects in ultra-thin liquid films in micropassages- comparison of thermodynamic theory with predictions of molecular dynamics simulations, J. Heat Transfer 128 (2006) 1276e1284. [34] V. Sartre, M.C. Zaghdouni, M. Lallemand, Effect of interfacial phenomena on evaporative heat transfer in micro heat pipes, Int. J. Therm. Sci. 39 (2000) 498e504. [35] R.H. Nilson, S.W. Tchikanda, S.K. Griffiths, M.J. Martinez, Steady evaporating flow in rectangular microchannels, Int. J. Heat Mass Transfer 49 (2006) 1603e 1618. [36] M. Ghajar, J. Darabi, Modeling of heat transfer characteristics of rectangular microgrooves in evaporators, in: AIAA Conference Paper 2007-4154, Miami, FL, 2007. [37] J.P. Longtin, B. Badran, F.M. Gerner, A one-dimensional model of a micro heat pipe during steady-state operation, J. Heat Transfer 116 (1994) 709e715. [38] M. Karthikeyan, J. Huang, J. Plawsky, P.C. Wayner Jr., Experimental study and modeling of intermediate section of the non-isothermal vapor bubble, Trans. ASME 120 (1998) 166e173. [39] R.K. Shah, M.K. Bhatti, in: S. Kacac, R.K. Shah, W. Aung (Eds.), Laminar Convective Heat Transfer in Ducts, Handbook of Single Phase Convective Heat Transfer, Wiley, New York, 1987. [40] Chin-Chun Hsu, Shung-Wen Kang, Tung-Fu Hou, Performance testing of micro loop heat pipes, Tamkang J. Sci. Eng. 8 (2) (2005) 123e132. [41] I.E. Dzyaloshinskii, E.M. Lifshitz, L.V. Pitaevskii, The general theory of van der Waals forces, Adv. Phys. 10 (1961) 165. [42] G.P. Peterson, J.M. Ha, Capillary performance of evaporating flow in micro grooves: an approximate analytical approach and experimental investigation, J. Heat Transfer 120 (1998) 743e751.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Evaporative heat transfer analysis of a micro loop heat pipe with rectangular grooves M. Ghajar a, J. Darabi b, * a b
Techlegic LLC, USA Department of Mechanical Engineering, Southern Illinois University Edwardsville, Edwardsville, IL 62026, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 July 2013 Received in revised form 16 December 2013 Accepted 19 December 2013 Available online 1 February 2014
This study focuses on the thermal and capillary analysis of a micro loop heat pipe for the thermal management of electronic devices and systems. A model is developed using the principles of thin film evaporation to predict the evaporative heat transfer coefficient in grooved capillary structures. In addition, a micro-flow submodel is developed to compute the dry-out distance in rectangular capillary grooves. These submodels are incorporated into our previously developed system-level loop solver model to simulate the performance of a flat micro loop heat pipe. The integrated model predicts the thermal performance, evaporator surface temperature, and local and average heat transfer coefficients as a function of the applied heat load. The modeling results are verified by comparison with the experimental data for a similar device and a good agreement is obtained. Ó 2014 Elsevier Masson SAS. All rights reserved.
Keywords: Micro loop heat pipe Evaporative heat transfer analysis Rectangular grooves Numerical model
1. Introduction Thermal systems miniaturization and the consequent rapid increase in power density of advanced microprocessors and electronic components have created a demand for achieving high heat dissipation rates. In most cases, such high heat fluxes cannot be easily dissipated using existing cooling techniques. Thus, research for new cooling technologies is critically important for space and defense applications. The quality of thermal management in electronic packaging directly affects the performance, cost, and reliability of electronic devices. A micro loop heat pipe (MLHP) can be a good candidate for such applications. Because no energy is needed to operate a loop heat pipe (LHP), it has been identified as an enabling technology for spacecraft thermal control. Conventional loop heat pipes have already been employed successfully in many applications, including spacecraft thermal control, electronics cooling, and ordinary commercial thermal devices. A loop heat pipe is a two-phase device, derived from heat pipe, which uses capillary forces developed inside its wicked evaporator to pump and circulate its working fluid in a closed loop, allowing heat transport over large distances at nearly constant temperatures. Major advantages of heat pipe include a very high thermal conductance, no pumping power requirements, no moving parts,
* Corresponding author. E-mail address: [email protected] (J. Darabi). 1290-0729/$ e see front matter Ó 2014 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ijthermalsci.2013.12.014
and relatively low pressure drops. More details on the principle of the heat pipes and loop heat pipes can be found elsewhere [1e7]. Loop Heat Pipes (LHPs) were invented in the former Soviet Union in the early 1980s [8,9]. LHPs have been used as baseline thermal control systems for a number of spacecrafts and commercial satellites due to their design simplicity and operation robustness [10e 15]. In addition to space applications, the high pumping head also makes LHP feasible in many terrestrial applications [12]. A comprehensive review of LHPs has been included in a previous publication of the authors [16] and other researchers [11e15]. A large number of analytical and experimental investigations have been conducted by various researchers to predict the thermal and fluid flow in conventional LHPs. However, a limited number of studies have focused on micro/miniature loop heat pipes [17e21]. Advances in microfabrication techniques have enabled the development of miniaturized, low-cost, and integrated MLHPs. A MLHP can be directly attached to the surface of an electronic component for which localized cooling is required. In an earlier work, the authors developed a FORTRAN-based system-level numerical model to simulate a planar MLHP and predict its thermal behavior as a function of the applied heat load [16]. Later, an improved model was developed to predict the evaporator surface temperature by assuming a typical evaporative heat transfer coefficient [22]. However, the evaporative heat transfer coefficient depends on several operational and geometrical parameters such as heat flux, evaporator temperature, fluid properties, and the geometry of the evaporator and the grooves. The
52
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
Nomenclature
l
A AR Dg h k K _ m Pvlv Pvlvt Q Qappl q_ q00 R RT T U V
n s
W x
dispersion coefficient, N m aspect ratio of the grooves distance between two grooves, m depth of the groves, m thermal conductivity, W/m K curvature, m1 mass flow rate, kg/s vapor pressure of the thin film, Pa saturation pressure at the film temperature, Pa heat transferred per unit length of the groove, W/m heat load, W heat flow rate, W/s heat flux, W/m2 ideal gas constant, N m/K thermal resistance, m2 K/W temperature, K fluid velocity, m/s phase velocity or velocity of fluid in phase change normal to the interface, m/s width of the groove, m coordinate along the groove, m
Greek symbols d film thickness, m q meniscus contact angle, rad
present model has been built upon prior work of the authors [23,24] to simulate the MLHP under various operating conditions and geometrical parameters. In this work, a submodel is developed to calculate the evaporative heat transfer coefficient using the principles of thin-film evaporation and incorporated into our previously developed systemelevel loop solver to further improve the simulation tool. Parameters of interest such as the local heat transfer coefficient, micro region characteristics, and dry-out distance were calculated iteratively using a procedure developed for the main loop solver [22]. 2. Model formulation Fig. 1 depicts a schematic illustration of a planar MLHP. The loop consists of an evaporator, a condenser, a compensation chamber, and liquid and vapor lines. A sectional view of the compensation chamber, liquid line, and evaporator, which are connected by an
coordinate perpendicular to the groove used in the micro-region, m kinematic viscosity, m2/s surface tension, N/m
Subscripts 0 initial appl applied ave average c capillary cc compensation chamber d disjoining e evaporator g grooves i interfacial iv interface (vapor side) l liquid ll liquid line lv liquidevapor n normal s surface sat saturation tot total T total v vapor w wall, wick, width
array of parallel microgrooves is shown in Fig. 1b. It is widely accepted that the temperature difference between the two ends of a wicking structure in a cylindrical evaporator can be estimated using the following equation [25,26]:
Te Tll ¼
dT dP
sat
DPtot DPg
(1)
where Te is the temperature of the capillary grooves adjacent to the evaporator, and Tll is the temperature of the capillary grooves adjacent to the liquid line, DPtot is the total pressure drop of the loop, and DPg represents the pressure drop of the wicking structure. Similarly, Eq. (1) can be used to estimate the temperature difference between the two ends of the wicking structure in a MLHP. However, for a MHLP, the presence of a conductive substrate can result in significantly higher heat conduction to the liquid line under a more complex transport pattern. The authors addressed this issue by
Fig. 1. (a) Schematic diagram of a planar MLHP; (b) sectional view (AeA)- the compensation chamber, liquid line, and evaporator are connected by microgrooves (not to scale).
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
estimating heat conduction to the compensation chamber using a finite difference method [16,22]. In this study, the principals of thin-film evaporation are incorporated into a micro flow solver to determine the wetted length of the groove during evaporation. Fig. 2 shows a schematic illustration of the liquid profile along a rectangular groove in the evaporator. As the liquid flows downstream towards the vapor line, the crosssection of the liquid film in the grooves deceases due to evaporation. At any given cross section, the evaporative heat transfer coefficient at the interlines is the main contributor to the local heat transfer coefficient. Thus, the local heat transfer coefficient increases as the number of interline regions increase. The interline regions of the thin-film are indicated by circles in Fig. 2b. An average normal evaporative heat transfer coefficient is defined as [27]
hn ¼
q00 ave Ts Tsat
(2)
where q00 ave is the average applied heat flux and Ts is the temperature at the inner surface of the evaporator. The vapor pressure in the groove is assumed to be constant and the heat transfer from the dry-out regions is assumed to be negligible. To calculate the substrate surface temperature, the local heat transfer coefficient must be known. The thin film evaporation theory can be used to calculate this parameter. The relationship between the vapor pressure of the thin film and the saturation pressure is given by [28]
ln
Pvlv Vl ¼ P Pvlvt RTSat d
(3)
where Pvlv and Pvlvt are the vapor pressure of the thin film and saturation pressure at the film temperature, and Pd is the disjoining pressure. For polar liquids, the disjoining pressure is given by [29]
Pd ¼
d ln 3
A
d
d0
(4)
where A is the dispersion coefficient and d is the liquid film thickness. The heat transfer through the liquid film is primarily dominated by conduction since the film thickness is very small. The interfacial resistance to heat transfer has been previously studied by several researchers [30e33] and is not repeated here for brevity. The following fourth order differential equation has been developed for the film thickness [27]. This equation requires four initial conditions for d. For instance, the value of d along with its three derivatives should be known at x ¼ 0.
hfg d Tw Tiv dP d3 c ¼ RT 3nl dl dl
53
(5)
Obtaining the required initial conditions for d to solve equation (5) is a challenging task. A recent approach by Sartre et al. [34] attempts to solve an alternative form of equation (5) in which the following system of differential equations are solved simultaneously:
dd 0 ¼ d dl dd ¼ dl
(6)
1þd
02 1:5
s
A Pc 3
d
(7)
dPc 3ml Q ¼ 3 dl hfg rl d
(8)
dQ ¼ ðTw Tiv ÞR1 T dl
(9)
The advantage of this approach is that the initial conditions at x ¼ 0 can be more intuitively understood as follows:
d ¼ d0 ;
K ¼ 0;
and
d0 ¼ 0
(10)
The radius of film curvature is constant in the intrinsic meniscus region and outside the microregion. This fact provides the remaining required boundary condition. The wall superheat temperature and the evaporative heat transfer coefficient are locally interdependent. Therefore, an iterative procedure was used to address this issue in which the values of the evaporative heat transfer coefficient and the wall temperature are updated alternatively. The evaporative heat transfer coefficient is calculated by the micro-region model while the wall temperature is determined by solving a steady-state heat conduction equation.
V2 T ¼ 0
(11)
The temperature variation in the solid domain is expected to be small for this microscale system. The alternate direction implicit (ADI) scheme was used in the finite difference scheme. All variables in the computer program were chosen as double precision. The mesh was refined to ensure convergence for a range of applied heats. Furthermore, analytical calculations and commercial software were used to verify the accuracy of the model for several cases. The system of equations (6)e(9) was solved numerically using a fourth order Runge-Kutta method. The meniscus profile was
Fig. 2. Thin film profile along the grooves in the evaporator.
54
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
constructed using the calculated meniscus radius and contact angle obtained from the microregion code. The local heat transfer coefficient of the thin film outside the microregion was calculated based on a 1-D conduction through the liquid meniscus. The localized temperature and evaporative heat flux information was obtained from the converged solution for the micro/macro regions. The local evaporative heat transfer coefficient was then calculated from:
hðlÞ ¼
q00 ðlÞ ; Tw ðlÞ Tiv ðlÞ
(12)
where q00 ðlÞ is the local heat flux and Tw(l) and Tiv(l) are the local wall temperature and local interface temperature in the evaporator. A cross-sectionally averaged heat transfer coefficient, h(x), was then obtained by integrating local heat transfer coefficient over l and was used in the loop solver. In this study, the triangular region in deep grooves was not included because of its short length compared to the rectangular region. For shallower grooves, it is necessary to consider both the rectangular and the triangular sections by combining the triangular (not the subject of this paper) and rectangular models. Fig. 3 shows a comparison between the triangular and rectangular geometries. Traditionally, for triangular grooves, the cross sectional areas for the liquid and vapor passage are expressed as functions of rm. This is because in a triangular geometry, the height of the liquid meniscus, Lw, can be geometrically related to the radius of curvature and contact angle of the meniscus. However, such a geometrical relation does not exist for a rectangular geometry. In a rectangular geometry, a similar analysis requires additional information about the liquid meniscus in the groove because rm is not related to the wetted height of the groove. Nilson et al. [35] proposed that the curvature of the meniscus in a rectangular micro groove only depends on the width of the groove. While this can be a reasonable assumption for wider rectangular channels, it does not necessarily result in an interesting solution for a micro groove for the following reason. It is acceptable to consider both as a fact and a precondition for the development of the above continuum-based formulations that the liquid film profile is a C1 continuous function. Therefore, it must be a smooth extension of the micro-region film thickness profile at the intrinsic meniscus junctions with the micro-region film. If the meniscus curvature is only a function of the groove width as in Refs. [35], the curvature at a point where the intrinsic meniscus starts to form, is constant for any wall temperature as long as the channel width remains constant. However, the
Fig. 4. Mass and momentum balance in a given cross-sectional element of the groove. The momentum terms due to shear stress at the interface are not shown.
contact angle which determines the slope of the microregion film profile is a function of the wall temperature. As a result, the C1 continuity condition at the junction cannot be satisfied. Thus, rm should be related to both the groove width and the contact angle. If the contact angle is available from the solution of the micro region model, the rm is known and the wetted height, Lw, can be calculated from the micro-flow model [36]. In the present work, a uniform meniscus radius is used by assuming a uniform contact angle obtained from the micro-region model in a single cross-section of the groove. This is a reasonable assumption for calculating the capillary pressure in microgrooves [3e5], allowing to determine the capillary driving force from the Young-Laplace equation. Ideally, the micro-region solver should be run for several cross-sections along the groove to calculate the contact angle and the meniscus radius. In a separate publication, the authors demonstrated that a constant contact angle assumption along the groove was valid [36]. The meniscus half-angle of curvature and the height of the meniscus can be geometrically determined using Fig. 3b:
Fig. 3. Comparison of liquid menisci formed in a triangular microgroove (a) and a rectangular microgroove (b).
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
q ¼ 2 sin1 Lw ¼
w 2rm
(13)
1 1 2 Al þ rm ðq sin qÞ w 2
(14)
The continuity and the momentum equations are derived for the one-dimensional two phase evaporating flow as shown in Fig. 4. The governing equations for this case are as follows:
d ½U A ¼ Vil Lci dx l l
(15)
d ½Uv Av ¼ Viv Lci dx
(16)
d P 1 Ul2 Al þ Al l ¼ ½slw ðw þ 2Lw Þ þ si Lci þ gAl sinð4Þ rl rl dx
d Pv U 2 Av þ Av rv dx v
¼
1
rv
½svw ðw þ 2Lw Þ si Lci
AT ¼ Al þ Av ¼ wh þ De w þ dg Pv ¼ Pl þ
s rm
rm
Lci ¼ 2rm q
where Lci is the length of the liquid surface curve at a given crosssection and 4 is the inclination angle of the grooved plate. The effect of the inclination is negligible in a micro loop heat pipe. Unlike a heat pipe, where the liquid and vapor flow in opposite directions, vapor flows concurrently with the liquid in a loop heat pipe. Therefore, the shear stress at the interface facilitates the flow of liquid due to a higher vapor velocity. The boundary conditions for the flow are similar to previous heat pipe studies [34,37,38]. The system of differential equations was solved using a Runge-Kutta scheme. The shear stresses at the wall and interface were calculated using the friction factor concept and known exposed areas at the wall and the liquidevapor interface. Equations (15)e(20) were solved to determine the liquid cross-section. The height of the meniscus at the wall was then calculated using Eq. (14). The friction factors for the grooves were calculated from Shah’s correlation [39]. The boundary conditions for liquid and vapor phase velocities are dependent of the liquid and vapor mass flow rates in the grooves [34,37,38].
rm ðrm q w cos qÞ 2
(17)
Al ðx ¼ 0Þ ¼ wh
(18)
Av ðx ¼ 0Þ ¼ AT eAl ðx ¼ 0Þ
(19)
Ul ðx ¼ 0Þ ¼
_ m whrl
(27)
Uv ðx ¼ 0Þ ¼
rl Ul Al rv Av
(28)
(20)
The unknowns are Ul, Uv, Al, Av, Pl, and Pv. The rm, Lci, Vil, and Viv are given in the following relations:
w ¼ 2 cos ðaÞ
55
Pv ¼ Pl þ s=rm
(21)
where Pl is an arbitrary base pressure.
(22)
3. Solution algorithm
Vil ¼
qf w rl hfg Lci
(23)
Viv ¼
rl V rv il
(24)
Fig. 5. Comparison of the predicted heat removal capability with the experimental data of Hsu et al. [40].
(25) (26)
(29)
The model consists of three parts. The first part which is based on the authors’ previous model [23] has remained mainly unchanged and is referred to as the “loop solver” in this paper. Briefly, the loop solver calculates the pressure drops in all segments of the loop and the temperature difference at the two ends of the wicking structure between the liquid line and the evaporator. Once, the liquid line temperature is known, a function is defined as the difference between the calculated temperatures from the available temperature field obtained from the finite difference model. By varying the saturation temperature, this function is minimized to determine the operating temperature. The second part is referred to as “the micro-region solver (MRS)” and provides an in-situ calculation of the evaporative heat transfer coefficient in the evaporator grooves. The MRS consists of two submodels. One submodel calculates the evaporative heat transfer coefficient and the other updates the wall temperature. The third part, which is referred to as the “micro-flow solver (MFS)” estimates the dry-out length by solving one-dimensional two-phase flow equations (continuity, momentum, and energy for both phases). An iterative procedure was used in a system-level loop solver model [23]. Using an initial wall temperature, the MRS determines the thin film profile and updates the heat transfer coefficient and the interfacial temperature. For a given applied heat flux, the wall temperature inside the evaporator is then calculated using a finite difference routine. This process is repeated until a converged solution is obtained for the wall temperature. After this procedure ends, the local heat transfer coefficient and wetted length are calculated from the micro-flow solver. The average heat transfer
56
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
coefficient in the grooved structures is determined by integrating Eq. (12) and used in the next iteration of the main loop solver.
Table 1 Geometrical dimensions used in the model e the base geometry. Parameter
Abbreviation
Dimension
Unit
Width of the groove Thickness of the groove wall Length of the groove Depth of the groove Length of the vapor line Length of the liquid line Width of the vapor line Width of the evaporator Length of the evaporator Depth of the liquid and vapor line, condenser, CC, and evaporator
W Dg Lg h Lv Lll Wv We Le hd
5 2 500 50 20 30 2 2 4 200
mm mm mm mm mm mm mm mm mm mm
4. Results and discussion The accuracy of the system-level model was verified by comparison with the experimental data of Hus et al. [40] using methanol as a working fluid. As shown in Fig. 5, the model predicts the experimentally measured evaporator temperature with good accuracy. Unless otherwise noted, the rest of the results in this section are presented for a micro loop heat pipe of the geometrical dimensions listed in Table 1 using methanol as the working fluid. It was demonstrated in a previous publication of the authors [22] that the device performance reached a saturation point as the heat transfer coefficient in the evaporator exceeded certain values. At very high heat transfer coefficients, the frictional forces present in the liquid and vapor lines and other components increase and have
Fig. 6. Profiles of various parameters in the micro-region using methanol as the working fluid.
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
57
Fig. 9. Comparison of the modeling results with the similar experimental data from Ref. [42]. Fig. 7. Variation of contact angle vs. DT ¼ Tw Tsat.
a negative impact on the overall device performance. Thus, the overall heat removal capability for a given operating temperature may not necessarily improve at high heat transfer coefficients. However, the evaporator wall temperature, which is the most critical parameter, is a strong function of the heat transfer coefficient inside the evaporator and can be greatly reduced with increasing the heat transfer coefficient. The results section consists of three separate parts pertaining to three parts of the model. The results of the micro-region model are presented first. In a micro groove, the groove length is several orders of magnitude larger than the length scale of its crosssection. Thus, the calculations are essentially based on the assumption of a one-dimensional heat transfer from the inner wall of the groove to the liquidevapor interface. In the second part, the micro-flow model predicts the groove behavior during evaporation using a one-dimensional hydrodynamic model. The purpose of this investigation is to ensure the grooves are capable of removing the applied heat load without forming a complete dryout region. Therefore, this sub-model is solved at each iteration of the loop solver and serves as a continuous examiner of the main model. This study investigates the behavior of the device at steady-state condition. A situation where a complete dry-out occurs in a section of the groove should be treated as a timedependent problem and is beyond the scope of this study. Finally, in the last part, the results of the integrated model which consists of the micro-region, micro-flow, and loop solver submodels are demonstrated.
Fig. 10. Variation of the wetted height along the groove at various groove depths. Please refer to Fig. 3b for definition of Lw.
4.1. Micro-region modeling results The results of the micro-region model using methanol as the working fluid are presented in Fig. 6. The general behavior of the results is similar to previously published data for ammonia [29,41]. As shown in Fig. 6a, the film thickness remains at the initial value to
Fig. 8. Variation of the heat transfer coefficient and wall temperature in a full cross-section of a 10 mm by 10 mm groove.
58
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
heat transfer to the evaporator is calculated in this part of the model and used to determine the mass flow rate. 4.2. Micro-flow modeling results
Fig. 11. Maximum tolerable heat flux as a function aspect ratio at various groove widths and a wall thickness of 5 mm.
mark the non-evaporating adsorbed layer before increasing sharply to specify the location of the interline. As shown in Fig. 6b, the distance with high evaporation heat flux, is limited to a few tenth of a micron which includes the interline. The interline region also corresponds to the region with a high evaporative heat transfer coefficient as shown in Fig. 6c. The significant part of the evaporation phenomena occurs in this region. Fig. 6d shows that the curvature reaches a maximum at the interline followed by formation of a constant curvature meniscus. The radius of this meniscus depends on the contact angle and the groove dimensions. As shown in Fig. 6e and f, the capillary pressure and the interface temperature behave quite similar since they are related linearly. The interface temperature approaches the saturation temperature in regions outside the interline. The variation of the contact angle as a function of the difference between the wall and the saturation temperature is illustrated in Fig. 7. As mentioned before, the contact angle data are used to determine the meniscus radius. The film thickness and subsequently the heat transfer coefficient outside the micro-region are determined by adopting a 1-D conduction approach and a known meniscus radius. The heat transfer coefficient and wall temperature profiles for a 10 mm by 10 mm groove are presented in Fig. 8. As expected, the minimum wall temperature occurs at locations with highest heat transfer coefficient values. The
In order to verify our model prediction and methodology, the simulation results were compared to available experimental data for an array of grooves with open triangular channels machined into a flat plate [42]. Fig. 9 shows a comparison between the modeling results and the experimental data for the dry-out distance as a function of the axial inclination angle. The results indicate that the dry-out distance is in good agreement with the experimental data and follows the same general trend. The liquid mass flow rate at a given cross-section along the groove decreases as it evaporates and enters the vapor phase. In addition, the liquid velocity decreases further downstream of the groove since the viscous forces overcome the inertia forces. Eventually, if the heat flux is high enough, a stationary condition is formed for the liquid phase at the dry-out point. Extremely high vapor velocities close to or exceeding the sonic limit are detrimental to the function of these devices and need be monitored in the design. Such a limit occurs at a heat load of approximately 6 W for the baseline dimensions listed in Table 1. The effect of the groove depth on the dry-out length is depicted in (Fig. 10). It is observed that as the groove depth increases, the dry-out occurs further away from the groove inlet. To better understand this phenomenon, the maximum heat flux that a groove can bear before a dry-out occurs as a function of the groove aspect ratio (depth/width) is shown in Fig. 11. The results indicate that for the range of the aspect ratios investigated, a groove with a width of 5 mm can deliver the highest heat flux before it experiences a dryout. However, for large aspect ratios, a narrower groove is expected to perform better since the capillary force is higher for a narrow groove due to smaller width, W (i.e. Pc w 1/W). For a constant width, a higher aspect ratio corresponds to a larger cross-sectional area. Thus, the friction forces decrease as the aspect ratio increases. This is the primary reason that the graphs approach constant values at large aspect ratios. The graphs intersect because far from the capillary limit of the device, a larger width (i.e. a larger crosssection) constitutes a lower contribution of the friction forces. 4.3. MLHP Modeling results Fig. 12 presents the cross-sectionally averaged heat transfer coefficient and maximum evaporator surface temperature as a function of the heat load. The results indicate that the crosssectionally averaged heat transfer coefficients decrease as the heat load increases. This is attributed to the change in the contact angle of the meniscus and subsequently a change in the meniscus radius and the thickness of the extended film. The evaporator surface temperature at a heat load of approximately 7 W exceeds the allowable temperature for electronics cooling applications (w100 C) for the geometrical dimensions used in this study. The heat removal capability is expected to improve by shortening or widening the vapor line in the loop. Also, our previous studies have shown that there exist optimum groove dimensions in the wicking structure. A parametric study similar to previous publications of the authors [16,22] can be performed to determine the optimal dimensions as well as the most suitable working fluids. 5. Conclusions
Fig. 12. Variation of the cross-sectionally averaged evaporative heat transfer coefficient and the maximum surface temperature of the evaporator as a function of the total heat load.
A new model was developed to simulate a micro loop heat pipe with rectangular capillary grooves. This model utilizes an in-situ calculation of the evaporative heat transfer using the principles of
M. Ghajar, J. Darabi / International Journal of Thermal Sciences 79 (2014) 51e59
the thin film evaporation. In addition, a micro flow solver submodel was incorporated into the model to estimate the length of the wetted evaporating region. This model is not intended to simulate situations where dry-out occurs in the evaporator. Such a case requires a transient solution rather than a steady state solution presented here. This integrated model which is built upon previously developed models by the authors is capable of predicting the operating characteristics and evaporative heat transfer coefficient of the micro loop heat pipe for any geometrical dimensions and operational parameters. As in the previous publications of the authors [16,22], the performance of the device is highly dependent on the geometrical dimensions and choice of working fluid. Thus, the optimization of the device geometry, particularly the groove dimensions are necessary to design an improved MLHP device. References [1] S.W. Chi, Heat Pipe Theory and Practice, Hemisphere Publishing Co, Washington, D.C, 1976. [2] P.D. Dunn, D.A. Reay, Heat Pipes, second ed., Pergamon Press, 1978. [3] G.P. Peterson, An Introduction to Heat Pipes: Modeling, Testing, and Applications, John Wiley & Sons, Inc., New York, 1994. [4] J.T. Dickey, G.P. Peterson, Experimental and analytical investigation of a capillary pumped loop, J. Thermophys. Heat Transfer 8 (3) (1994) 602e607. [5] A. Faghri, Heat Pipe Science and Technology, Taylor & Francis, Inc, Washington DC, 1995. [6] L.L. Vasiliev, Micro and miniature heat pipes e electronic component coolers, Appl. Therm. Eng. 28 (4) (2008) 266e273. [7] S. Launay, V. Sartre, J. Bonjour, Parametric analysis of loop heat pipe operation: a literature review, Int. J. Therm. Sci. 46 (7) (2007) 621e636. [8] Y.F. Maidanik, S.V. Vershinin, V.F. Kholodov, Y.F. Dolgirev, US Patent No. 4515209, May 7, 1985. [9] V.M. Kiseev, A.G. Belonogov, N.P. Pogorelov, Development of two-phase loops with capillary pumps, SAE 972482, July 14e17, 1997. [10] Y.F. Maidanik, Y.G. Ferchtater, K.A. Goncharov, Capillary Pump Loop for the Systems of Thermal Regulation of Spacecraft, 1991. ESA SP-324. [11] Y.F. Maidanik, Y.G. Ferchtater, Theoretical basis and classification of loop heat pipes and capillary pumped loops, in: 10th International Heat Pipe Conference, Stuttgart, Germany, 1997. [12] J. Ku, Operating characteristics of loop heat pipes, SAE Paper (1999). No. 199901-2007. [13] D.A. Wolf, D.M. Donald, A.L. Phillips, Loop heat pipes e their performance and potential, in: SAE Paper No. 941575, 24th ICES, Germany, 1994. [14] A.L. Phillips, J.E. Fale, N.J. Gernert, D.B. Sarraf, W.J. Bienert, Loop Heat Pipe Qualification for High Vibration and High-G Environments, 1998. AIAA-98e 0885. [15] A. Hoelke, H.T. Henderson, F.M. Gerner, M. Kazmierczak, Analysis of the Heat Transfer Capacity of a Micromachined Loop Heat Pipe, International Mechanical Engineering Congress and Exposition, IMECE, HTD-Vol. 364-3, ASME Publication, Nashville, Tennessee, November 15e20, 1999, 53e60. [16] M. Ghajar, J. Darabi, N. Crews Jr., A hybrid CFD-mathematical model for simulation of a MEMS loop heat pipe for electronics cooling applications, J. Micromech. Microeng. 15 (2005) 313e321. [17] M.A. Chernysheva, S.V. Vershinin, Y.F. Maydanik, Operating temperature and distribution of a working fluid in LHP, Int. J. Heat Mass Transfer 50 (2007) 2704e2713. [18] R. Singh, A. Akbarzadeh, M. Mochizuki, Operational characteristics of a miniature loop heat pipe with flat evaporator, Int. J. Therm. Sci. 47 (11) (2008) 1413e1562.
59
[19] S.F. Wang, W.B. Zhang, X.F. Zhang, J. Chen, Study on start-up characteristics of loop heat pipe under low-power, Int. J. Heat Mass Transfer 54 (4) (2010) 1002e1007. [20] J. Li, D. Wang, G.P. Peterson, Experimental studies on a high performance compact loop heat pipe with a square flat evaporator, Appl. Therm. Eng. 30 (6e7) (2010) 741e752. [21] X. Zhang, X. Li, S. Wang, Three-dimensional simulation on heat transfer in the flat evaporator of miniature loop heat pipe, Int. J. Heat Mass Transfer 54 (2012) 188e198. [22] M. Ghajar, J. Darabi, Numerical modeling of evaporator surface temperature of a micro loop heat pipe at steady-state condition, J. Micromech. Microeng. 15 (2005) 1963e1971. [23] M. Ghajar, J. Darabi, Integration of a thin-film evaporation submodel with a micro loop heat pipe solver, in: Proceedings of the IMECE2006, 2006 ASME International Mechanical Engineering Congress and Exposition, November 5e 10, 2006. Chicago, Illinois, USA. [24] M. Ghajar, J. Darabi, Simulation of steady state surface temperature response of a novel micro loop heat pipe to various heat loads, in: Proceedings of the ASME Summer Heat Transfer Conference, HT, 2005, pp. 633e640. [25] T. Kaya, T. Hoang, Mathematical Modeling of Loop Heat Pipes, 1999. AIAA Paper 99-0477. [26] T. Hoang, T. Kaya, Mathematical Modeling of Loop Heat Pipes with Two-phase Pressure Drop, 1999. AIAA Paper 99e3448. [27] P.C. Stephan, C.A. Busse, Analysis of heat transfer coefficient of grooved heat pipe evaporator walls, Int. J. Heat Mass Transfer 35 (2) (1992) 383e391. [28] B.V. Derjaguin, Z.M. Zorin, Optical study of the adsorption and surface condensation of vapors in the vicinity of saturation on a smooth surface, in: 2nd Int’l Congr. Surf. Activity (London), 2, 1957. [29] P.G. deGennes, Wetting: statics and dynamics, Rev. Mod. Phys. 57 (3) (1985). Part 1. [30] P.C. Wayner Jr., Y.K. Kao, L.V. LaCroix, The interline heat transfer coefficient of an evaporative thin film, Int. J. Heat Mass Transfer 19 (1976) 487e492. [31] B.V. Derjaguin, S.V. Nerpin, N.V. Churayev, Effect of film transfer upon evaporation of liquids from capillaries, Bull. Rilem 29 (1965). [32] J.G. Truong, P.C. Wayner Jr., Effects of capillary and van der Waals dispersion forces on equilibrium profile of a wetting liquid: theory and experiment, J. Chem. Phys. 87 (7) (1987) 4180e4188. [33] V.P. Carey, A.P. Wemhff, Disjoining pressure effects in ultra-thin liquid films in micropassages- comparison of thermodynamic theory with predictions of molecular dynamics simulations, J. Heat Transfer 128 (2006) 1276e1284. [34] V. Sartre, M.C. Zaghdouni, M. Lallemand, Effect of interfacial phenomena on evaporative heat transfer in micro heat pipes, Int. J. Therm. Sci. 39 (2000) 498e504. [35] R.H. Nilson, S.W. Tchikanda, S.K. Griffiths, M.J. Martinez, Steady evaporating flow in rectangular microchannels, Int. J. Heat Mass Transfer 49 (2006) 1603e 1618. [36] M. Ghajar, J. Darabi, Modeling of heat transfer characteristics of rectangular microgrooves in evaporators, in: AIAA Conference Paper 2007-4154, Miami, FL, 2007. [37] J.P. Longtin, B. Badran, F.M. Gerner, A one-dimensional model of a micro heat pipe during steady-state operation, J. Heat Transfer 116 (1994) 709e715. [38] M. Karthikeyan, J. Huang, J. Plawsky, P.C. Wayner Jr., Experimental study and modeling of intermediate section of the non-isothermal vapor bubble, Trans. ASME 120 (1998) 166e173. [39] R.K. Shah, M.K. Bhatti, in: S. Kacac, R.K. Shah, W. Aung (Eds.), Laminar Convective Heat Transfer in Ducts, Handbook of Single Phase Convective Heat Transfer, Wiley, New York, 1987. [40] Chin-Chun Hsu, Shung-Wen Kang, Tung-Fu Hou, Performance testing of micro loop heat pipes, Tamkang J. Sci. Eng. 8 (2) (2005) 123e132. [41] I.E. Dzyaloshinskii, E.M. Lifshitz, L.V. Pitaevskii, The general theory of van der Waals forces, Adv. Phys. 10 (1961) 165. [42] G.P. Peterson, J.M. Ha, Capillary performance of evaporating flow in micro grooves: an approximate analytical approach and experimental investigation, J. Heat Transfer 120 (1998) 743e751.