Microchannel Heat Sinks Using Anisotropic Porous Media Approximation Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

Recent advances in semiconductor technology have led to significant increase in power densities in microelectronic equipment. According to the International Roadmap for Semiconductors (ITRS), the maximum power of single chip package will reach 173 W/cm2 in 2022 (ITRS, 2008). The growing thermal issue is causing distress due to shorter product lifetime, increased of build cost and increased of operating cost (Bailey, 2008). The increase of heat generated within the microelectronic equipment has led to the proposal of myriad cooling methods. One of the promising possible methods for removing high heat flux is by using forced convection microchannel heat sinks which are fabricated as an integral part of the silicon wafers (Li and Peterson, 2007).

The microchannel heat sinks were firstly proposed over 27 years ago by Tuckerman and Pease (1981, 1982). In their work, they constructed a silicon wafer and achieved maximum substrate temperature rise of 71°C with heat flux of 790W/cm2 using water as coolant. They found that the convective heat transfer coefficient, h at the substrate-coolant interface was the primary factor in thermal resistances reduction. However, their proposed optimization method was subjected to several constraints; the channel flow was assumed to be fully developed and laminar, while the Nusselt number, fin efficiency and pumping pressure were fixed.

Since this seminal paper, many numerical and experimental studies had been conducted. In 2007, Wei et al. (1984) presented experimental and numerical results of stacked microchannel heat sinks for liquid cooling. In their research, they studied the effects of flow direction and flow ratio in each microchannel and identified that within the range of the tested flow rate, counter flow arrangement provides better temperature uniformity, while parallel flow has the best performances in reducing the peak temperature. Recently, Ogedengbe et al. (2008) observed that axial conduction temperature has significant effect on the reduction of velocity gradient along the microchannel and the entropy production level is minimized as the flow approaching the fully-developed regime with counter flow-heated microfluidic wall. Dixit et al (2008) has proposed design of silicon micro/nanopillars based multilayer water-cooled heat sink that has enhanced the overall thermal performance of electronic. From their analytical results, it showed that the heat dissipation rate of heat sink with silicon pillars can be increased by 16% than that without silicon pillar. However, the pillars will easily collapse if the bending stress due to the flow exceeds the fracture strength of silicon (~7GPa).

Experimental works were conducted by Missaggia et al. (1989), Kleiner et al (1995), and Lasance et al. (2005). Theoretical and numerical works were performed by Samalam (1989) and Weisberg (1992) who investigated classical one-layer microchannel heat sink. Philips (1990) and Morini (2004) presented critical review paper on the numerical study and experimental work. However, most of the studies were limited to the conventional methods utilizing CFD codes. Even if it is feasible, mapping the whole microchannels can lead to very high computing cost which places a limit on the number of simulated channels and the dimension of apparatus. Such method request more extensive numerical effort to handle complex geometries which may involves other types of channel for fluid transportation such as right, obtuse or even sharp angle bends with occasionally rounded corners which are also common in microfluidic chips. (You et al., 2005). Shojaeian and Dibaji (2010) have simulated the slip flow regime for fully developed incompressible fluid flow and heat transfer through triangular microchannels. Their results reveal that the rarefaction decreases the Poiseuille number and the interaction between velocity slip and temperature jump will influence the Nusselt number.

In latter development, Koh and Colony (1986) modelled the microchannel heat sink as fluid-saturated porous medium, describing the flow using Darcy's law. This approach was later on developed by Tien and Kuo (1987) using the modified Darcy equation, which accounts for the boundary effect on the convection problems. In their approach, a two equation model which treats the fluid and the solid phase individually was utilized. Their work were extended by Kim and Kim (1999), whom introduced analytical solution for fluid flow and temperature distribution using analytical heat coefficient and permeability values. Their simulation results were in good agreement with the experimental undertaken, h, highlighting the successful use of the porous medium approach for thermal optimization and the prediction of the thermal resistance of the microchannel heat sinks.

In 2000, Kim and Kim (2000) presented a one equation model for heat transfer. However, the results obtained from the one equation model can be practically applied to microchannels only when the heat sink is highly porous and the thermal interaction between the fluid and solid phase is highly effective. Both papers by Kim and Kim (1999, 2000) have successfully predicted the fluid flow and temperature distribution based on analytical determination of the permeability and interstitial heat transfer coefficient. However, both papers were limited with several constraints: (1) the flow was assumed as fully developed, (2) the entrance effect of the flow which is important was neglected, (3) the heat coefficient was fixed, although temperature distribution was varying along different location along the microchannels. Do et al (2010) has proposed a modified similarity solution velocity and temperature distribution using fluid-saturated porous medium concept in the heat sink subject to a uniformly impinging jet. In their results, it has shown that the modified similarity solutions are suitable for predicting the pressure drops and the thermal resistance of the heat sink subject to uniformly impinging jet.

Recently, Amanifard et al (2007) and co workers modeled microchannels as fluid-saturated porous medium based on the Forchheimer-Brinkman-extended Darcy equation and the two equation model for heat transfer between solid and fluid phase. They found that the Nusselt number increased with the increasing aspect ratio. Their results also gave the required assurance of using the full Navier-stokes approach for the microchannels with hydraulic diameters about 100µm. However, their works were still limited to laminar, hydrodynamically and thermally fully developed flows.

In 2009, Ghazazvini and Shokouhmad (2009) have conducted a research study on using nanofluid as coolant of a microchannel heat sink. Fin model and the porous media approach were used to model the microchannel heat sink. To calculate the nanofluid heat transfer, a model based on Brownian-motion of nanoparticle was used. However, their comparison of their result was is imperfect because the fin model approach assumes that heat coefficient was constant, although the heat coefficient was varying along microchannel that will results in errors in predicting heat transfer of microchannel.

In this paper, the fluid flow and heat transfer for a single microchannel was simulated to investigate temperature contour and characteristic of the fluid flow. Simulation results for fluid flow for the single microchannel were then compared with prediction from the Brinkman equation model. To further the validation of the simulation, the results were examined and compared to the results done by Kim and Kim (1999-2000). The proposed method was later compared with the experimental results of Tuckerman (1984) using the same parameters and dimensions in their work.

Problem Description and Model Definition

In this study, water was used as the working fluid in the microchannel heat sinks operation, while silicon was used as the material of the microchannel. The microchannel has a width, Wch and depth, H. The walls that separate the channels are referred to as fin and has a width, Wfin. The top wall of the microchannels is assummed to be effectively insulated while the bottom wall receives a uniform heat flux, qw. The heat flux represents the heat needed to be dissipated from the electronic chip. Coolant water is forced to enter the open channel in the x-direction to achieve heat dissipation as shown in Fig. 1. The tb is the thickeness of silicon; whs is the total width of microchannels while Lhs is the total length of the microchannels.

Fig 1 Schematic diagram of the microchannel heat sink.

The water phase is assumed and treated as incompressible single phase-fluid. The steady-state governing equation can be written in vector form as follows:






, for fluid; (3)

, for conduction in substrate. (3)

The governing equations were solved using the commercial code COMSOL Multiphysics, which uses the finite element method in solving the governing equations.

In order to study the fluid flow and heat transfer of the single microchannel, the model employs the general steady convective heat transfer analysis with incompressible fluid as the medium. The domain is made of two sub-domains consisting of a cuboid and a L-shape as depicted in Fig. 2. Silicon was used to construct the microchannel structure. The dimensions for the single- microchannel is tabulated in Table 1 which corresponds to the aspect ratio H/Wc = 6 to allow comparison with the results obtained by Kim and co-workers (1999-2000). Due to symmetry of the problem, only half of the single microchannel was selected as computational domain as shown in Fig 3.

In order to ensure simulation accuracy, quadrilateral elements were used to generate the volume mesh, as shown in Fig.3 for the single microchannel configuration.

Fig. 2 Geometry for the Single Microchannel

Fig.3 Mesh for microchannel

Table 1 Dimensions of the microchannel.

Dimension (mm)











As proposed by Kim and co-workers (1999, 2000), the microchannel heat sink can be modelled as a porous media. The governing equation for the fluid flow and heat transfer is established by applying volume-averaged technique as shown below (Kim and Kim, 1999, 2000):

Fluid flow:


Energy in solid phase:


Energy in fluid phase:


Hence, the microchannel heat sink shown in Fig. 1 was modelled as a porous medium with porosity, ε, permeability, κ and the wetted area per volume, a can be expressed as (Bejan, 1984):


The parameters above were used in the simulation to model the porous media. Within the medium, the porous subdomains, the flow variables and fluid properties at the centroid of each mesh cell were defined as the value averaged over the volume of the mesh cell.

To investigate the fluid flow of microchannel as porous media, a steady-state simulation based on the Brinkman equation was performed. A configuration of microchannels was generated as shown in Fig. 4, with the dimension 1cm ´ 1cm ´ 300µm representing a single microchannel heat sink as fluid saturated porous media. The Brinkman equation was then employed throughout the whole porous media using the parameters proposed by Bejan (1984).

Fig. 4 Mapped quadrilateral mesh for microchannels

Physics of the Model

The material properties used to model the single microchannel are listed in Table 2 while the material properties of the porous media are listed in Table 3. The parameter for porosity, permeability and thermal conductivity were calculated as follow:




Table 2 Properties for single Microchannel


Thermal conductivity (W/m-K)

Specific heat capacity (J/kg-K)

Density (kg/m3)


(at 296K)








Table 3 Properties for porous media


Thermal conductivity (W/m-K)

Specific heat capacity (J/kg-K)

Density (kg/m3)

Porous media








These values were used for the porous medium simulation. In order to compare the single microchannel and porous media model, two cases have been conducted. For the first case, the same flow rate was used which is 6.2GPH and for the second case, simulations were undertaken under the limitation of pressure drop of 200kPa. The results were compared with the existing analytical solution (Kim and Kim, 1999, 2000).

The same dimension model of single microchannel and porous media was used to allow more effective comparison later on in section 4. For example, the same channel height and length were used in the single microchannel model and the porous media model. The aspect ratio of six was used in this paper.

Single Microchannel

Boundary Conditions

The bottom part of the L-shape was set to impart a constant convective heat flux of 0.42 MW/m2. The outer surface temperature of the boundary was set at the ambient temperature of 293 K. These settings were to mimic typical heat flux of a similar flip-chip. Temperature was specified at the inlet, while the outlet and the symmetry plane of the fluid domain were set as convective flux and thermal insulation, respectively.

For case 1, the inlet velocity was set at 2.17m/s and pressure outlet was set at 0 Pa respectively. The symmetry part was set to symmetry. All other parts were set to no slip condition.

For case 2, the pressure inlet and pressure outlet was set at 200kPa and 0 Pa respectively. The symmetry part was set to symmetry. All other part was set to no slip condition. From some preliminary calculations, the Reynolds numbers in both cases are lower than 1980 which justifies the consideration of laminar flow in both microchannel configurations.

Model Implementation and Computation

In this numerical study, mapped mesh was firstly generated at the inlet surface of the microchannel. Boundary layer mesh was also generated at the wall of the top and bottom part. Later on, the boundary layer mesh was generated with sufficient refinement to ensure the accuracy of the velocity gradient prediction near the wall. The remaining domain meshes was later swept with 10 number of element layers. Lagrange quadratic elements were used during all the meshing operation. Using a mapped mesh shown in Fig. 3, the microchannel model now consists of 3690 hexahedral elements. The simulation was executed in stationary mode.

Porous media

Boundary Conditions

Similar to the single microchannel, the bottom part was set at 0.42MW/m2. These boundary conditions were mention in section 3.1.1

For case 1, the inlet velocity was set at 4.35m/s. All other part was set to no slip condition.

For case 2, the pressure inlet was set at 200kpa and pressure outlet was set at 0Pa. All the other parts were set to no slip condition.

Model Implementation and Computation

For the porous media, a mapped quad mesh was generated as shown in Fig 4. A total of 264 elements were generated at the inlet surface and the surface was later swept with 8 number of element layers. Using this mapped mesh, the final meshed geometry consisted of 2,112 hexahedral elements.

Solving Method

Since the microchannel heat sink is a combination of many singular rectangular shaped microchannels, it is essential to solve the porous media in one direction. This implies that the porous media is anisotropic which allows the fluid to flow in only one direction. This can be done by solving the three dimensional flow by keeping the y-axis velocity and z-axis velocity to a small value which is almost zero. The proposed method will improve the accuracy of the fluid flow and heat transfer predictions. This directional treatment of the velocities was only used in the solution of the fluid flow, and was not applied in the solution of the heat transfer module. A comparison between the simulation using the ordinary Brinkman equation porous media and anisotropic Brinkman equation porous media is shown in Fig. 5. This simulation consists of one inlet at the left bottom and another outlet at the right top.

Results and Discussion



Fig. 5 Three Dimensional and Two Dimensional view of simulation of the (a) ordinary Brinkman equation porous media and (b) Anisotropic Brinkman equation porous media with x-direction.

From Fig. 5, the ordinary Brinkman was set as the benchmark for anisotropic porous media. Using the proposed method as mentioned earlier, it is obvious that we can simulate the flow rate in one direction which represents the true fluid flow behavior of microchannels heat sink. Setting the same mesh density, parameter and geometry, it took 427.266s computer time to solve the ordinary Brinkman equation while for the anisotropic brinkmann equation, it only took 49.45s to solve the fluid flow in microchannels. The proposed method can reduce the computing time by approximately one order of magnitude.

Since the inlet velocity or pressure setting described in section 3 was very high, it is impossible to visualize the entrance effect for fluid flow. Therefore, lower entrance value was used to visualize the ability of the proposed method to account for the entrance effect of fluid flow for both cases.

The results show that the entrance effect is accounted by the simulation. In Fig. 6 the characteristic of flow can be easily illustrated. The flow starts with uniform velocity at the inlet and becomes fully developed at the location in the range of 0.0015m - 0.002m.

For the fluid flow in porous media, the slice view was taken at the y-axis in the geometrical center. In Fig. 6, the flow starts to be fully developed at the range of 0 -0.001m. From Fig. 7, it can be clearly seen that the fluid flow in the porous media model takes shorter length to be fully developed than fluid flow in the single microchannel. This is due to the higher velocity set at the entrance of the porous media.

Fig. 6 Fluid flows at the Entrance of Single Microchannel (0.001m/s)

Fig. 7 Fluid flows at the Entrance of Porous media Microchannel (0.1m/s)

Fig. 8 shows the temperature distribution along the x- axial of single microchannel while fig. 9 shows the temperature distribution along the x- axial of porous media. Both of the models showed that the highest temperature occurs at the outlet of the microchannel.

Comparison between Single Microchannel and Porous Media:

Fig. 8 Temperature Distribution of Single Microchannel (200kPa)

Fig. 9 Temperature Distribution of Porous media Microchannel (200kPa)

Table 4. Comparison between single microchannel and porous media under the same pump capacity.


























Table 5. Comparison between single microchannel and porous media under the same pressure drop limitation




















Porous Media





Table 4 shows the comparison between the porous media and the single microchannel under the same pump capacity while Table 5 shows the comparison between the porous media and the single microchannel under the same pressure drop limitation. In treating the microchannels as fluid saturated porous media, the mean temperature and the maximum temperature can be accurately predicted with the discrepancy of 0.25-0.38 percent. In this section, prediction using the porous media approach reached good agreement with the single microchannel prediction. Overall thermal performance of the porous microchannel can be accurately predicted as shown in Table 4 and 5.

Analytical Solution

In this section, the result of case 2 will be validated with the analytical solution by Kim and Kim (1999, 2000) .The same aspect ratio of six and similar material were used for comparison. All the results reported were non-dimensionalized using the following dimensionless variable,


, where um is the mean velocity in fluid region, αs is the aspect ratio of the microchannel, H/Wc and Y is the dimensionless vertical coordinate. The value of Da was taken as 2.31 x 10-3.

For the single microchannel, um was calculated by dividing the volume (integration of the fluid domain) with the velocity field expression. The expression is shown in (12)


The value for um is 3.72m/s. Velocity data was extracted from the point (5mm,0µm,100µm) to point (5mm,0µm,400µm) which is at the mid of the microchannels length. To obtain the average in microchannel, the following relation was used:

(Poiseuille flow for parallel plate) (13)

For the porous media, um was also calculated by dividing the volume (integration of the fluid domain) with the velocity field expression following equation (12). The velocity data was extracted from the point, (5mm, 5mm, 100 µm) to the point (5mm, 5mm, 400µm). The value obtained for um was 1.83m/s.

The analytical result (Kim and Kim, 1999, 2000) for the dimensionless velocity and temperature is given by:






The mathematical expression for the coefficients appearing in the above solutions can be found in Kim and Kim (1999,2000).



Fig. 10. Comparison velocity with solutions from conventional methods (Da =2.31 x 10-3).

Fig. 10(a) depicts the results for the non-dimensionless velocities for the developed flow profile across the microchannels along the x-axis direction. Since the velocity profile is symmetry, only half of the fluid flow was plotted.

From Fig. 10, the Brinkman equation model deviates with an average of 2.84 percent from the results of Kim and Kim (1999, 2000). On the other hand, simulation from the single microchannel exhibited an average deviation of 1.08 percent from that of Kim and Kim (1999, 2000). The present results show good agreement with the numerical results of Kim and Kim numerical result which has been validated with numerical experimental result.

Using the correct permeability, the velocity of the fluid flow in statured porous medium can be obtained. Further comparison was made on the dimensionless temperature, as shown in Fig. 10(b). For the single microchannel, the temperature profile was been extracted from the point (5mm, 0, 100µm) to (5mm, 0, 400µm). For the porous media, temperature data was extracted from the point, (5mm, 5mm, 100 µm) to the point (5mm, 5mm, 400µm). From the results, it can be seen that temperature distribution for a single microchannel deviated significantly from that of Kim and Kim. This is because the result was obtained for only the fluid part without considering the solid part of the single microchannel. In the porous media, the temperature distribution is obtained by assuming the temperature of the fluid and solid phase are the same; the one equation and two equation model for heat transfer are use to validate the proposed method. In Fig. 10(b), the porous media approach exhibited similar results with that of Kim and Kim for Equation 16 and Equation 17 which represent the one equation for heat transfer and solid part for the two-equation model.

Fig. 11 Comparison temperature distribution porous media at different location with conventional methods

Fig. 11 shows the result of the dimensionless temperature distribution in the cross section of the porous media at different location at 0.002, 0.004, 0.006, 0.008 and 0.010m. The average of 5 of the location was plotted to validate with the result of Kim and Kim. From Fig. 11, the porous media approach shows good agreement with that of Kim and Kim with an average of only 3.83 percent deviation in the prediction.

The results showed that the temperature distribution varies along the microchannels. The analytical results by Kim and Kim showed a good average of temperature but it did not represent the exact temperature along the microchannel. The proposed method has the ability to predict the temperature distribution at different locations. This proposed method will be very important because it allows visualization of the temperature distribution which will be essential in dealing with hot spots in semiconductor packaging.

Besides that, the results from Kim and Kim were based on permeability and the interstitial heat transfer coefficients which were determined analytically. In the present study, we only need to determine the permeability, which represents the viscosity in the microchannels; and heat coefficient, which may vary was calculated during the simulation based on the porous media properties.

Experimental Results

The proposed method was validated with the results of Kim and Kim. However, the validation was only limited to velocity and temperature distribution along the z-axis. Therefore, it is necessary to compare the overall thermal performance with the available experimental data.

In the present work, overall thermal resistance for porous media has been compared with the experiments of Tuckerman (1984). The same dimension and properties of sample 81F9 were chosen for comparison purpose, shown in Table 6. The pressure drop and overall package dimension were maintained and the overall thermal resistance was calculated using the expression below:


Table 6. Comparison of thermal resistances between present study, and Tuckerman and Pease experimental results.




Present Study

Channel Dimensions

L (mm)



W (mm)
















Pressure Drop



Heat Flux (W/cm2)





Resistance (°C/W)



As shown in Tables 6, the values of total thermal resistance calculated from the present analysis are in good agreement with those obtained by Tuckerman (1984).


In the preceding analysis, a full three-dimensional simulation for fluid flow and heat transfer for single microchannel heat sink was developed. To validate the proposed porous media approach, comparison was made with past analytical and experimental research findings. The non-dimensionless form has been used in the comparison. The analysis indicated that the fluid flow predictions agreed well with that of past reports. In treating the microchannel as a fluid saturated porous media, the reasonable volume-averaged velocity and temperature distribution prediction were obtained. Thus, it may be concluded that the present methodology is able to predict the total thermal resistance of any microchannel heat sinks satisfactorily. In addition, the proposed method is not limited or constrained by assumptions on the laminar or fully developed nature (even thermally) of the flow. Thus, this approach can also be used developing flows. The porous media approach has the advantage in simulating the fluid behavior of microchannels with lower computation requirements when compared to the conventional single microchannel approach.


The authors gratefully acknowledge the financial support of the Ministry of Science, Technology and Innovation, Malaysia under the project 03-01-02-SF0318.