Numerical modelling of heat transfer in rectangular microchannels

The paper presents both three and two-dimensional numerical analysis of convective heat transfer in microchannels. The three-dimensional geometry of the microchannel heat sink followed the details of the experimental facility used during a previous research step. The heat sink consisted of a very high aspect ratio rectangular microchannel. Two channel heights, namely 1mm and 0.3mm (0.1mm), were used for 3D (2D) numerical model respectively. Water was employed as the cooling liquid. The Reynolds number ranged from 200 to 3000. In the paper, the thermal entrance effect is analyzed in terms of heat transfer efficiency. Finally, the comparison between measured and computed heat flux and temperature fields is presented. Contrary to the experimental work, the numerical analysis did not reveal any significant scale effect in heat transfer in microchannel heat sink up to the smallest size considered (0.1 mm).


OBJECTIVES
Microchannel heat sinks are of special interest as high efficiency cooling devices due to undeniable advantages such as small coolant demand and dimensions, which directly result in high convective heat transfer coefficient.Some experimental studies [Peng et al. (1994), Qu et al. (2000), Gao et al. (2002)] indicate that the reduction of size may affect heat transfer in microchannels.These scale effects are difficult to interpret and the physical mechanisms are still not clarified (Guo and Li, 2003).Several numerical studies [Weisberg et al. (1992), Fedorov and Viskanta (2000), Qu and Mudawar (2002a, b), Toh et al. (2002), Tunc and Bayazitoglu (2002)] have been published in recent years on this subject.Among the possible sources of deviation from the laws governing the macroscale flows, effects of axial heat conduction in the walls and entrance effects in the microchannel test section are seriously suspected to contribute to the observed scale effect.The present paper devoted to the numerical modelling of the flow and heat transfer in microchannels is focused on the coupling of axial heat conduction in the walls with convection in the microchannels.Associated entrance effects are of special interest, because they are probably significant in most published data.
This study pertains to experimental investigations (Gao et al., 2002) conducted in large-span rectangular microchannels, which have revealed a significant reduction in the Nusselt number for microchannel height less than 500 µm.The geometry considered in the present numerical simulations is as close to the real situation as possible.It is therefore possible to compare numerical and experimental results and to test the relevance of conduction/convection coupling effects and associated entrance effects in the deviation of heat transfer results from the macroscale laws.

EXPERIMENTAL SET-UP
In the experimental set-up, the active channel walls were two hand-polished plane bronze blocks, which were separated by a stainless plate (of thickness e) with a hollowed out central part of width w (see dimensions in Fig. 1).Details of the facility may be found in Gao et al. (2002).The height e of the microchannel test section was determined by the plate thickness, which could be varied in the range mm 1 1 .0 ÷ by steps of 0.1 mm.Heating was provided by four electric cartridges, which were inserted inside the two blocks and surrounded by an insulating material.The inlet/outlet water temperatures (T in , T out ) were measured inside the upstream/ downstream containers, then the fluid bulk temperature T f distribution in the microchannel was approximated as a linear one.It is worth noting that the heater is slightly shorter than the microchannel length.Four thermocouples were mounted along the microchannel in the plane of symmetry with their tips located 1mm above the microchannel surface.Two pressure transducers were flushmounted in the wall of the inlet/outlet containers.Thus, the pressure at the microchannel entrance was corrected by the inertia term accounting for flow acceleration in the converging channel entrance.The data were interpreted by using the Reynolds number Re and the Nusselt number Nu, both based on the hydraulic diameter D h .
Since measurements of the local heat flux ϕ were not possible in the experimental study of Gao et al. (2002), the authors assumed that the power density dissipated into the test section was uniformly distributed along the active microchannel surfaces.Results from the experimental investigations showed that the longitudinal distribution of Nu as given by Shah and London (1978, referenced as SL hereafter) is well recovered for the channels of largest height in the laminar regime.However, for e < 500 µm, a significant decrease in Nu was observed for reduced channel size.

NUMERICAL MODEL
The geometry of microchannel heat sink shown in Fig. 1 followed the details of the experimental setup.
Taking advantage of the symmetry of the experimental facility, the numerical model was represented by half of the actual geometry.Heat transfer process occurring in such a complex geometry consists of coupled conduction in the solid and convection to the cooling fluid.The conventional Navier-Stokes and energy equations were simplified by the following assumptions: (1) fluid incompressibility (2) laminar flow character (3) negligible radiation heat transfer (4) negligible buoyancy (5) constant fluid and solid properties (6) negligible viscous heating Because of the assumptions of constant fluid properties and negligible buoyancy, the mass and momentum equations were not coupled with the energy equation.Therefore, the temperature field was calculated by solving the energy equation after the converged solution for the flow field was obtained.Such an approach allowed to significantly speed up the computational convergence.The momentum and continuity equations were solved for the following hydraulic boundary conditions: (1) uniform velocity profile was set at the inlet to the test section, (2) the flow was assumed to be fully developed at the test section exit plane (3D model) or at the channel exit plane (2D model).
The thermal boundary conditions are given as follows: (1) constant heat flux density at the walls of the heater cartridges (3D model, total heat flux Φ = 180 W), constant energy source inside the cartridge volume (2D model), (2) free convection heat transfer with surrounding air at the outside walls of the test section; the assumed values of the heat transfer coefficient and ambient temperature were equal to 10 W/m 2 K and 300 K, respectively.
The x-direction runs along the channel.The origin is at the channel entrance after the converging region.Numerical computations of the flow were carried out by using the commercial code Fluent 5.4.The equations were discretized by means of a second order upwinding finite volume method.Meshing was performed by using Gambit 1.3.Due to very complex geometry, the computation domain was split into several regions.
The domain corresponding to solid block was meshed using unstructured mesh.
Structured mesh was utilized for inlet/outlet containers and microchannel.The total number of cells was optimized from series of tests and equal to about 2x10 6 (2.7x10 6 ) for e = 1 (0.3) mm.The number of mesh nodes used in the microchannel was equal to 164 (300) x 20 x 30 for e = 1 (0.3) mm, respectively.Special care had to be taken with the cell aspect ratio especially in the entrance region of the test section, i.e. for the rapidly developing flow.In order to avoid a too large number of nodes, the grid applied in the microchannel entrance was non-uniformly distributed in x-direction.

Model validation
The numerical model was checked for laminar simultaneously thermally and hydraulically developing flow along a two-dimensional channel with symmetrical uniform heat flux surfaces.We used for comparison the results of Hwang and Fan (1964) reported by SL and fitted by Bejan and Sciubba (1992) using empirical formulas of the Churchill-Usagi (1972) type.For Pr = 0.7, the local Nusselt number is given by ( ) where the dimensionless abscissa is Pr Equation ( 1) assumes uniform inlet velocity and temperature profiles, i.e. it takes into consideration the socalled entrance effect.
Figure 2 shows an excellent agreement of the present simulation with Eq. ( 1).This provides confidence in the numerical technique and therefore in the results produced with the current numerical approach.

RESULTS AND DISCUSSION
Hydraulic properties Shah and London (1978) presented a comprehensive review of laminar forced convection in ducts of various cross-section, including rectangular channels.Their results will be used as the reference in this study.
The Fanning friction factor f is related to the pressure drop ∆p over the hydrodynamically developing length x via The Poiseuille number is defined as For laminar developing flow in two-dimensional channels of usual size, Shah and London (1978) proposed the following formula Figure 3 shows the variations of the Poiseuille number against x + for the two microchannel heights considered in this study.
For the complex geometry of the numerical simulation (including inlet container), Po values are significantly smaller than predicted by SL.This is obviously due to the finite thickness of the boundary layers at the channel entrance, as opposed to uniform velocity profile in the SL reference case.Furthermore, the distribution of Po depends Fig. 2. Variation of Nusselt number along the channel (Pr=0.7).Comparison with the formula of Bejan and Sciubba (1992).
on Re and channel height separately, contrary to SL reference case where the proper dimensionless group is x + .As a result, the distributions of Po(x + ) for both channel heights are different as shown in Fig. 3.It can also be observed that the pressure drop increases with channel height, although both distributions tend to the asymptotic value of laminar fully developed flow (Po = 24).
Heat transfer Temperature field.Figure 4 shows the temperature distribution in the plane of symmetry.The region of the highest temperature gradients corresponds to the presence of the microchannel insulation.It is worth emphasizing that insulation is not complete around the bronze blocks.There is a gap near the microchannel inlet allowing a part of the heat flux to be conducted to the surrounding blocks.That in turn leads to: (i) additional heat transfer to water in the inlet container and (ii) increased heat losses to ambient air.
Figure 4 shows that the wall temperature decreases in the last part of the microchannel (x/L > 0.73) due to the lack of heating in this region (see Fig. 1 to compare heater cartridge and microchannel lengths).
Nusselt number.The Nu number is defined as where the convective heat transfer coefficient h is given by Two-dimensionality checks showed that the wall heat flux ϕ num varied less than 5 % over two-third of the span.As a consequence, the wall heat flux, the wall and local fluid bulk temperatures (respectively T w , T f ) are taken in the plane of symmetry of the flow.
Figure 5 illustrates the variations of Nu against the dimensionless coordinate defined by Eq. ( 2).The numerical results are compared with the data obtained for uniform inlet velocity and temperature profiles and Pr = 6.This reference curve allows to identify and quantify the entrance effect.The Nusselt number was calculated along the channel for three values of the Reynolds number.Each curve in Fig. 5 corresponds to a particular value of Re.As can be seen in Fig. 5, the particular curves are clearly Re-dependent revealing relative shifts and intersections.The differences between these curves can obviously be explained by the Re-   The data may also be presented for a fixed value of x and varying Re, as in the experimental procedure.As an example, Figure 5 shows the variation of Nu for x/L = 0.44 (solid symbols).The comparison between the reference curve and the marked numerical points for this particular position in the microchannel allows for the following observations: (1) the entrance effect is weaker for the present results when e = 1 mm (2) the convergence of data points and the reference curve when e = 0.3 mm indicates that the entrance effect is so weak at this position in the channel that Nu is no more sensitive to the inlet conditions (3) the entrance effect is reduced for lower channel height Further explanation should be given for the first observation.
A non-uniform temperature profile was observed at channel inlet because a proportion of the total heating power was delivered to water in the converging part of the inlet container.
The close agreement of the data with the reference curve for x* > 0.005 suggests that the numerical results achieved for complex geometry can be physically interpreted and justified.In consequence, the numerical results were recognized as correct.

Experiment vs simulation
Figure 6 shows the distribution of the dimensionless temperature θ, defined as along the channel both for fluid and solid domains.The numerical solution for θ is considered at the location of the thermocouples.
Comparison between experimental results and numerical predictions for the water bulk temperature (squares in Fig. 6) shows that discrepancies are relatively small.The assumption of linear variation of T f is rather well verified.One can see that greater discrepancies occur in the solid region (circles in Fig. 6) resulting in higher values of the difference f w T T − for experimental results.For e = 1 mm, these discrepancies vary from about 5 % to 18 % with reference to experimental data, while for e = 0.3 mm they are larger and reach about 40 %.In this latter case, they account for the decrease in Nu observed in the experiments.
The numerical heat flux ϕ num was normalized by the value ϕ exp corresponding to uniform distribution over the heating surface, which was assumed during experiments.The three following regions may be observed on Fig. 7: (1) the microchannel entrance where strong values of ϕ num occur owing to very narrow velocity and temperature boundary layers in this region (2) the middle part of the microchannel where the ϕ num distribution can be assumed as approximately uniform (3) the exit part (around 30 % of the channel length) where ϕ num decreases due to the lack of heating and to conduction in the solid blocks The assumption used in the experiments neglected heat transfer in inlet and outlet containers.This explains why the numerical heat flux is smaller than the experimental one and, moreover, is non-uniformly distributed.
The location of the thermocouple tips could noticeably decrease the values of Nu, especially for small microchannel heights, owing to finite temperature gradient in the solid blocks.Nu was calculated in the experiments (Gao et al., 2002) with the temperature given by the thermocouples, which were about 1 mm away from the liquid/solid interface.In order to eliminate the influence of the thermocouples position in the determination of Nu, the experimental temperature was corrected using the numerical solution.The experimental wall temperature was estimated by using the wall temperature gradient in the solid.Furthermore, the actual wall heat flux is probably close to the numerical value.
Consequently, the experimental Nusselt number was recalculated with the numerical heat flux for want of a better estimation of ϕ.For e = 1 mm, which may be considered as a conventional size for a channel flow, above comparisons indicate that: (1) the wall temperature is slightly underestimated by the numerical model (Fig. 6a) (2) the heat flux is significantly overestimated (about 25 %) by the assumption used in the experiments due to neglected heat transfer in inlet/outlet containers (Fig. 7) (3) Nu is slightly overestimated (about 10 to 20 %) by the numerical model when Nu exp is calculated with ϕ num .
For e = 0.3 mm, comparisons show that: (1) the wall temperature is much higher than the numerical solution (Fig. 6b) (2) the heat flux ϕ num distribution is the same as in the previous case (3) The Nusselt number Nu exp is much lower (about 40 to 80 %) than the numerical value.
An important conclusion is therefore that the observed reduction in Nu exp cannot be accounted for by conduction effects, which would have distorted the temperature field and the corresponding heat flux distribution.
As a consequence, the reduction in Nu exp is not explained by the present numerical simulation.Further investigations are planned to take into account the variations of the fluid properties with tempearature in the numerical model.

2D NUMERICAL SIMULATIONS
The simulation performed for 2D geometry could bring obvious advantages as reduction in computer running time and simplification of modeled geometry.This section compares results obtained in 2D and 3D numerical simulations.The computations were achieved for the two extreme channel heights used in the experiments, i.e. 0.1 mm and 1 mm, with the corresponding Re number ranges: Re ∈ (200 ÷ 2166) for e = 1 mm Re ∈ (500 ÷ 3000) for e = 0.1 mm Figure 9 shows an excellent agreement between Nu distributions obtained with the 2D and 3D models.This suggests an application of 2D model for further quantitative investigations.

Hydraulic properties
3D simulations have shown that the Poiseuille number should be determined as a function of Re for given channel length L and height e, since the entrance effect was found to be Re-dependent.As a consequence, the computations were achieved for given values of L and e and varying Re, contrary to 3D simulations.It is, however, presented in Fig. 10 against L + , in consistency with the reduction of experimental data.Comparison between numerical and experimental data.
As can be seen in Fig. 10, the numerical results show a consistent distribution of Po(L + ) which is located below SL curve, as predicted by the 3D model and explained previously.
The entrance effects are found to be negligible when L + is higher than about 0.1 with corresponding values of the Re number: Re = 300 for e = 1 mm Re = 3000 for e = 0.1 mm The two sets of numerical results then overlap.

Heat transfer
Figure 11 shows the variations of Nu along the channel for the same two extreme channel heights (e = 1 mm and e = 0.1 mm).Computations were achieved by varying the Reynolds number as for hydrodynamics (ranges are given in previous subsection).Nu was estimated in the sections of the four thermocouples with the temperature at the fluid/solid interface.One can observe that Nu values are well fitted to the reference curve except for data corresponding to the section of the first thermocouple (x/L = 0.073, open circles in Fig. 11).For the conditions of the numerical simulation, Nu is therefore affected on a short distance by the finite thickness of the thermal boundary layer at the channel entrance.

CONCLUSIONS
The three-dimensional flow and associated heat transfer in a rectangular micro-channel heat sink were analyzed numerically.The geometry of the test section used in the experiments of Gao et al. (2002) was taken into account in 3D computations.The numerical simulation considered the coupling between convection in microchannels and conduction in the walls and in the complete solid material.The local heat transfer characteristics i.e. temperature, heat flux and Nusselt number, were obtained.The key findings of the present study are as follows: (1) solutions of numerical simulation using the continuum model (conventional Navier-Stokes equations and conventional energy equation) are in satisfactory agreement with published data on flow  and heat transfer in two-dimensional channels (2) thermal properties show that entrance effects are dependent on Reynolds number and channel height separately, contrary to Shah and London (1978) reference case (3) the numerical simulations turned out to be very helpful as a complement to the interpretation of experimental data, where complex measurements of heat flux and temperature field are obviously not possible (4) the numerical models used in this work assumed some simplifications e.g. did not take into account viscous heating, which could decrease the Nusselt number especially for channels of smaller heights than those investigated in the current work (5) the numerical simulations show that there is no size effect on heat transfer when the channel height is reduced from 1 mm down to 0.1 mm.As a result, the strong reduction in the Nusselt number observed in experiments cannot be explained by conduction effects due to the complex geometry, like axial conduction in the walls or lack of twodimensionality of the heat flux distribution.

Fig. 1 .
Fig. 1.Sketch of the test section

Fig. 3 .
Fig. 3. Variation of Poiseuille number along the channel.Influence of inlet velocity profile.

Fig. 5 .
Fig. 5. Variation of Nusselt number along the channel for two channel heights: a) e = 1 mm, b) e = 0.3 mm.

Fig. 6 .
Fig. 6.Variation of wall and fluid temperatures along the channel.Comparison between numerical and experimental results for two channel heights: a) e = 1 mm, b) e = 0.3 mm.
Figure 8 compares the numerical and experimental results of Nu.It should be noted that only the second and third thermocouples (x/L = 0.25 and 0.44) were considered in the experimental Nu values because they correspond to the region of nearly uniform heat flux.

Fig. 7 .
Fig. 7. Variation of dimensionless wall heat flux along the channel.

Fig. 11 .
Fig. 11.Variation of Nusselt number along the channel for two channel heights: a) e = 1 mm b) e = 0.1 mm.