Performance of a flat grooved heat pipe with a localized heat load

Heat pipes are phase change heat transfer devices used in wide range of heat transport applications due to their high thermal transport capacities with low temperature differences. Heat pipes are especially preferred for electronic cooling applications and aerospace avionics to satisfy high heat transfer rate requirements. In this study, heat transfer and phase change mechanisms of working fluid are investigated and modeled using a 3-D thermal resistance network for multichannel flat grooved heat pipes. First, heat transfer and fluid flow are modeled in half of a single grooved structure due to symmetry, and is subjected to uniform heat flux. Radius of meniscus curvature and temperature distribution along the groove are calculated. Results are compared with experiments in the literature and show good agreement. The validated heat transfer and fluid flow models are extended to a multichannel model to observe performance of grooved heat pipes with localized heat sources, not covering the entire width, a vital feature for realistic simulation of operational devices. Predictions of the temperature distribution along the multichannel of the heat pipe are provided and the effect of the distribution of heat sources on the heat pipe is discussed.


Introduction
By the recent advances in technology, electronic components are getting smaller in size, and the corresponding power dissipated is getting higher and it becomes an important and critical issue for the proper utilization of these devices. Besides traditional heat transport devices and methods, heat pipes have been considered as a favorable alternative [1] not only for thermal management of electronic components [2] but also in space applications [3], heating, ventilating and air-conditioning systems [4] and nuclear applications [5].
As a result of increasing heat pipe utilization, a better understanding of heat pipe characterization is necessary to assess and improve thermal performance. However, it is a challenging task since heat pipes characterization involves simultaneous modelling of different phenomena including phase change free surface capillary flow and heat transfer modeling. Various studies have been done to predict the thermal performance of flat grooved heat pipes that suggest different phase change, flow and heat transfer models. Do et al. [6] developed a mathematical model for predicting the thermal performance of a flat micro heat pipe with a rectangular grooved wick structure. In their study, the effects of liquid-vapor interfacial shear stress, contact angle and the amount of liquid charge were investigated and it was concluded that some common assumptions used in previous studies may be misleading while predicting thermal performance of the heat pipe. * Corresponding author: refaz@metu.edu.tr Lefèvre et al. [7] developed a two-phase flow model and resistance network based thermal model to calculate the liquid and vapor pressures and velocities, the meniscus curvature radius and the temperature in the heat pipe container from the source to the sink region. Heat conduction in each cross section in liquid and solid regions was simulated using thermal resistances which were utilized to calculate the axial temperature distribution. By using two-phase model vapor velocities, the liquid and vapor pressures and the meniscus curvature radius are determined and used for calculate transversal resistances in heat transfer model. Odabaşı [8] developed a model using the 3-D heat transfer equations in both the solid and the liquid, coupled with a simplified 1-D momentum equation. A simplified form of the momentum equation along the heat pipe axis was formulated to calculate the variation of liquid vapor interface radius along the heat pipe which generates the capillary force necessary to drive the flow. Phase change heat transfer from micro region was calculated using the relation obtained from kinetic theory and phase change heat transfer both from micro region and macro region were included in the analysis. Zhang et al. [9] proposed a 1-D thermal resistance network to optimize the performance of a heat pipe with Ω-shaped grooves. In order to optimize the performance, the effect of structural parameters and total thermal resistance on heat transfer capability were investigated. Özçatalbaş et al. [10] investigated performance of flat grooved heat pipes under transient heat loads by using finite element method. They proposed a special purpose user subroutine to define effective thermal conductivity of heat pipe using experimental data provided in literature for different heat loads.
Beside analytical and numerical studies, numerous experimental studies concentrated on performance of heat pipes are conducted to investigate in literature. Alijani et al. [11] investigated the effect of filling ratio on thermal performance for four different heat fluxes experimentally by defining new three performance indicator. It is shown that groove geometry and filling ratio had an important effect on performance. Atay et al. [12] experimentally tested the performance of heat pipes with sintered and grooved wick structures for different heat loads, ambient temperatures and gravitational orientations. Dry-out limits of selected heat pipes were reported to guide thermal engineers for selection proper heat pipe on their application.
In this study, a 3-D thermal resistance network accounting for both the axial and transverse conduction heat transfer and a 1-D simplified flow model are developed and combined to obtained vapor temperature, contact angle and temperature distribution over the flat grooved heat pipe. Phase change mechanisms are introduced into model while calculating the related resistance value for evaporation and condensation. After the validation of the developed model for the uniformly distributed heat load, thermal resistance network is extended to allow the simulation of the performance of flat grooved heat pipes under non-uniform heating and cooling conditions. The effect of localized heat sources and heat sinks can be observed by applying the relevant boundary conditions to the corresponding nodes.

Heat transfer modeling
Heat transfer in a rectangular channel of flat grooved heat pipe is modeled by using thermal resistances for the steady-state operation. Channel geometry for a flat grooved heat pipe and liquid distribution both in the evaporation and condensation zones are shown in Figure  1.
By considering the geometry, loading condition and boundary conditions along the heat pipe, symmetry along the axial direction is used and the resistance network is constructed for a half groove and a half fin depicted with dashed lines in Figure 1. Flat grooved heat pipe (FGHP) divided into several control volumes and energy conversation equations are written for each node in every control volume. Nodes in each control volume are connected with corresponding thermal resistances. In addition, axial conduction through the solid is accounted by longitudinal resistances ܴ to model heat transfer between axially neighboring control volumes. Resistance network for control volumes with length, ∆y is shown in Figure 2.
Since the evaporation and condensation have different characteristics in terms of heat transfer, it should be checked whether evaporation or condensation takes place in the control volume in question. This is decided by comparing the wall temperature and vapor temperature, ܶ ௩ in each control volume.
Conservation equation are written for each node by using appropriate resistances. Evaporation and condensation phenomena are implemented to model by using the phase change resistances, namely ܴ ହ in evaporation and ܴ ଵ in condensation regions. These resistances are calculated by using phase change mass flow rates from the corresponding areas.
Evaporation mass flux in the thin film region can be found by equating mass flux obtained from kinetic theory to mass flux written in terms of heat flux [8]. To do so, evaporation mass flux can be found as Where ߜ is liquid film thickness, ܿ is accommodation coefficient, ‫ܯ‬ is the molecular weight, ℜ is the universal gas constant. In general, ܽ(ܶ ௪ − ܶ ௩ ) ≫ ܾ(ܲ − ܲ ௩ ) and assuming the liquid film thickness to vary linearly for heat conduction purposes, evaporating mass flux can be integrated over the liquid film thickness to obtain evaporating resistance, ܴ ହ .
Condensation mass flux in fin top region can be calculated by assuming a 4 th order polynomial profile for the fin top thickness, ߜ subjected to boundary conditions given below [8].  ݀δ/݀s = ݀ ଷ δ/d‫ݏ‬ ଷ = 0 at s = 0 ݀δ/݀s = −tan (π/2 − θ) at s = ‫ݓ‬ ݀ ଶ δ/d‫ݏ‬ ଶ = 0 at s = ‫ݓ‬ Similar to evaporation mass flux, using the fact that mass flux obtained from kinetic theory is equal to mass flux obtained from heat flux relation. An iterative procedure is employed to obtain the constant term of 4 th order profile. Once the fin top thickness profile is obtained, condensation mass flow rate per unit length, ݉ ᇱ can be calculated by taking following integral over fin width.

Flow modeling
A simplified form of the momentum equation along the FGHP is used to obtain the variation of liquid-vapor interface radius which generates the necessary capillary driving force for the flow. Force balance for liquid control volume can be written as follow where ‫ܣ‬ and ߩ is liquid cross sectional area and liquid density, respectively. ‫ܣ‬ ௩ and ߬ ௩ is cross sectional area and shear between liquid-vapor interface, respectively. ‫ܣ‬ ௪ and ߬ ௪ is cross sectional area and shear between the liquid and groove wall, respectively. Since the phase change occurs all along the heat pipe, changing in interfacial shear stress between liquid and vapor can be neglected [7]. Liquid pressure change can be obtained by taking the derivative of the Young-Laplace equation where ߪ is surface tension and assuming changing in vapor pressure along the heat pipe can be neglected. Using a force balance, the Young-Laplace equation and mass conservation for the control volume, the variation of the radius of curvature, which generates the capillary force, can be found.

Validation of the developed model
Validation of the current model is done by the comparing the two parameters, namely radius of curvature and wall temperature values with the numerical and experimental results provided in [7].
FGHP used in study is made of copper and working fluid is methanol. Dimensional and physical properties used in the model is provided in Table 1. All dimensions given in mm. The maximum difference for wall temperatures between current study and Lefèvre et al. [7] is about 1.3°C. Vapor temperature in [7] is given as 70°C and in the present study this is calculated as 69.8°C. Temperature distribution along the FGHP is given in Figure 4. The results for radius of curvature are given in Figure 5. Both wall temperatures and radii of curvature show good agreement with the numerical and experimental results given by Lefèvre et al. [7].

Extended multichannel model
In real life operations, heat pipe may not be exposed the uniformly distributed heat load all along its width as assumed in validated model. To simulate more realistic heating and cooling conditions encountered in applications, half-groove model extended to multichannel model that allowed to apply localized heat sources not covering the entire width of FGHP. A part of extended resistance network is shown Figure 6. FGHP with N number of grooves divided into M number of control volumes to construct resistance network analogy. Each half-groove model connected to its lateral neighbor with red colored resistance shown in Figure 6. As mentioned in the half groove model, each longitudinal neighbor volume also linked with each other with green colored longitudinal resistances to simulate axial heat conduction. In order to observe effect of localized heat source on heat pipe's performance, two different simulation is conducted. First simulation has a heat source and heat sink regions not covering the entire width of FGHP as seen in Figure-7. Heater with a dimension 40 x 20 mm located at the FGHP wall. Heat removing from the system is achieved using a heat sink dimensions of 50 x 28 mm. Physical parameters and dimensions are provided in Table-2. All dimensions given in mm. In the second simulation, nothing is changed but only it is assumed that heat source and heat sink are completely covering the width of FGHP with the same heat flux and lengths ‫ܮ‬ and ‫ܮ‬ , respectively. Aluminum and isopropyl alcohol are selected for FGHP and working fluid, respectively.

Results and discussion
In order to capture the effect of localized heat source, a multichannel model is simulated. Interested FGHP is made of aluminium with length 150 mm with heat source and sink which are not covering the entire width of FGHP as seen in Figure 7. Heat load q"=17000 W/m2, is applied as heat flux boundary condition to the nodes which stay inside the indicated red region in Figure 7. Similarly, convective heat transfer boundary condition with convective heat transfer coefficient, h=250 W/m2K and ambient temperature, T=10°C is imposed to nodes inside the blue region shown in Figure  7. Obtained temperature plot for FGHP is given in Figure 8. In longitudinal direction, max temperature difference is calculated as 14.3 °C between evaporator and condenser end of FGHP. In lateral direction, maximum temperature difference is equal to 6.4 °C.
In evaporation zone, the hot spot region can be seen around the heater area. Although the effect of heat source is getting decrease as moving away from it, it conducted through the base material, in this case it is aluminum, in the lateral direction.
In condensation zone, a non-uniform temperature distribution is also observed. However, temperature gradient is lower compare to evaporation zone due to wider heat sink region.
Temperature distribution through the FGHP for the second simulation with heat source and sink covering entire width of FGHP with the same length. ‫ܮ‬ and ‫ܮ‬ is provided in Figure 9. In Figure 9, a uniform temperature distribution in lateral direction is observed, as expected. In axial direction, temperature difference between evaporator and condenser end is calculated as 15.4 °C. Indeed, effect of localized heat load on thermal performance is monitored clearly from the temperature surface plots provided in Figure 8 & 9. Proposed extended model is able predict local hotspot regions and can be used to determine location and size of both heater and sink for reliable heat remove applications.