INTRODUCTION
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 cm^{2} in 2022 (http://www.itrs.net/).
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 790 W cm^{2} using water
as coolant. They found that the convective heat transfer coefficient, h at the
substratecoolant 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 study, many numerical and experimental studies had been
conducted. Wei et al. (2007) presented experimental
and numerical results of stacked microchannel heat sinks for liquid cooling.
In their research, they studied the effect 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 fullydeveloped regime with counter flowheated microfluidic wall. Dixit
et al. (2008) has proposed design of silicon micro/nanopillars based
multilayer watercooled 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 (~7 GPa).
Experimental works were conducted by Missaggia et al.
(1989), Kleiner et al. (1995) and Lasance
and Simons (2005). Theoretical and numerical works were performed by Samalam
(1989) and Weisberg (1992) who investigated classical
onelayer 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 fluidsaturated 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.
Kim et al. (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 study by Kim and Kim (1999) and Kim
et al. (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 fluidsaturated 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.
Amanifard et al. (2007) modeled microchannels
as fluidsaturated porous medium based on the ForchheimerBrinkmanextended
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
Navierstokes approach for the microchannels with hydraulic diameters about
100 μm. However, their works were still limited to laminar, hydrodynamically
and thermally fully developed flows. 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 Brownianmotion 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 study, 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) and Kim et
al. (2000).
The proposed method was later compared with the experimental results of Tuckerman
(1984) using the same parameters and dimensions in their study.
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, W_{ch} and depth, H. The walls that separate the channels are referred to as fin and has a width, W_{fin}. The top wall of the microchannels is assumed to be effectively insulated while the bottom wall receives a uniform heat flux, q_{w}. 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 xdirection to achieve heat dissipation as shown in Fig. 1. The t_{b} is the thickness of silicon; W_{hs} is the total width of microchannels while L_{hs} is the total length of the microchannels.
The water phase is assumed and treated as incompressible single phasefluid.
The steadystate governing equation can be written in vector form as follows:
Mass:
Momentum:
Energy:
For fluid:

Fig. 1: 
Schematic diagram of the microchannel heat sink 
For conduction in substrate:
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 subdomains consisting
of a cuboid and a Lshape as depicted in Fig. 2. Silicon was
used to construct the microchannel structure. The dimensions for the single
microchannel is shown in Table 1 which corresponds to the
aspect ratio H/Wc = 6 to allow comparison with the results obtained by Kim
and Kim (1999) and Kim et al. (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.
Table 1: 
Dimensions of the microchannel 


Fig. 2: 
Geometry for the single microchannel 

Fig. 3: 
Mesh for microchannel 
As proposed by Kim and Kim (1999) and Kim
et al. (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 volumeaveraged technique as shown (Kim and
Kim, 1999; Kim et al., 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 steadystate simulation based on the Brinkman equation was performed.
A configuration of microchannels was generated as shown in Fig.
4, with the dimension 1 cmx1 cmx300 μ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).
Physics of the model: The material properties used to model the single microchannel are shown in Table 2 while the material properties of the porous media are shown in Table 3.
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 7.83 cm^{3}
sec^{1} and for the second case, simulations were undertaken under
the limitation of pressure drop of 200 kPa.

Fig. 4: 
Mapped quadrilateral mesh for microchannels 
Table 2: 
Properties for single microchannel 

Table 3: 
Properties for porous media 

The results were compared with the existing analytical solution (Kim
and Kim, 1999; Kim et al., 2000). The same
dimension model of single microchannel and porous media was used to allow more
effective comparison. 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 study.
Single microchannel Boundary conditions: The bottom part of the Lshape was set to impart a constant convective heat flux of 0.42 MW m^{2}. 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 flipchip. 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.17 m sec^{1} 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 200 kPa 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 1,980 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 3,690 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.42 MW m^{2}. For case 1, the inlet velocity was set at 4.35 m sec^{1}. All other part was set to no slip condition. For case 2, the pressure inlet was set at 200 kPa and pressure outlet was set at 0 Pa. 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 yaxis velocity and zaxis 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.

Fig. 5: 
Three dimensional and two dimensional view of the simulation.
(a) Ordinary brinkman equation porous media and (b) anisotropic brinkman
equation porous media with xdirection 
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
Ordinary brinkman equation porous media: Anisotropic Brinkman equation porous media with xdirection.
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. Upon setting the same mesh density, parameter and geometry, it took 427.3 sev computer time to solve the ordinary Brinkman equation while for the anisotropic brinkmann equation, it only took 49.5 sec to solve the fluid flow in microchannels. The proposed method can reduce the computing time by approximately 90% reduction.

Fig. 6: 
Fluid flows at the entrance of single microchannel (0.001
m sec^{1}) 

Fig. 7: 
Fluid flows at the entrance of porous media microchannel (0.1
m sec^{1}) 

Fig. 8: 
Temperature distribution of single microchannel (200 kPa) 
Since, the inlet velocity or pressure setting described in previous section 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 visualized from the simulation as depicted in Fig. 6. The flow starts with uniform velocity at the inlet and becomes fully developed at the location in the range of 0.00150.002 m. By analyzing the entrance effect, it can be observed that the increase of normalized friction might cause by entrance effect. To eliminate the entrance effect, it is suggested that inlet pressure should be measured when the flow becomes fully developed.
For the fluid flow in porous media, the slice view was taken at the yaxis in the geometrical center. In Fig. 6, the flow starts to be fully developed at the range of 00.001 m. 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.
Figure 8 shows the temperature distribution along the xaxis
of single microchannel while Fig. 9 shows the temperature
distribution along the xaxis of porous media.

Fig. 9: 
Temperature distribution of porous media microchannel (200
kPa) 
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 

Both of the models showed that the highest temperature occurs at the outlet
of the microchannel.
Comparison between single microchannel and 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.250.38%. 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: Here, the result of case 2 will be validated with
the analytical solution by Kim and Kim (1999) and Kim
et al. (2000), which has been used as the benchmark for the remaining
analysis. The same aspect ratio of six and similar material were used for comparison.
All the results reported were nondimensionalized using the following dimensionless
variable,
where, u_{m} is the mean velocity in fluid region, α_{s} is the aspect ratio of the microchannel, H/W_{c} and Y is the dimensionless vertical coordinate. The value of Darcy number, Da was taken as 2.31x 10^{3}.
For the single microchannel, u_{m} was calculated by dividing the volume (integration of the fluid domain) with the velocity field expression. The expression is shown as below:
The value for u_{m} is 3.72 m sec^{1}. Velocity data was extracted from the point (5 mm, 0 μm, 100 μm) to point (5 mm, 0 μm, 400 μm) which is at the mid of the microchannels length. To obtain the average ⟨u⟩_{f} in microchannel, the following relation was used:
Poiseuille flow for parallel plate:
For the porous media, u_{m} was also calculated by dividing the volume (integration of the fluid domain) with the velocity field expression following Eq. 9. The velocity data was extracted from the point, (5 mm, 5 mm, 100 μm) to the point (5 mm, 5 mm, 400 μm). The value obtained for u_{m} was 1.83 m sec^{1}.
The analytical result (Kim and Kim, 1999; Kim
et al., 2000) for the dimensionless velocity and temperature is given
by:

Fig. 10: 
Comparison velocity with solutions from conventional methods
(Da = 2.31x10^{3}) 
The mathematical expression for the coefficients appearing in the above solutions
can be found by Kim and Kim (1999) and Kim
et al. (2000).
Figure 10a shows the results for the nondimensionless velocities
for the developed flow profile across the microchannels along the xaxis 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% from the results of the selected benchmark. On the other hand,
simulation from the single microchannel exhibited an average deviation of 1.08%
from the results of the selected benchmark. The present results show good agreement
with the numerical results of the selected benchmark 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. 10b. For the single microchannel,
the temperature profile was been extracted from the point (5 mm, 0, 100 μm)
to (5 mm, 0, 400 μm). For the porous media, temperature data was extracted
from the point, (5 mm, 5 mm, 100 μm) to the point (5 mm, 5 mm, 400 μm).

Fig. 11: 
Comparison temperature distribution porous media at different
location with conventional methods 
From the results, it can be seen that temperature distribution for a single
microchannel deviated significantly from results of selected benchmark. 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. 10b, the
porous media approach exhibited similar results with Kim and Kim for Eq.
14 and 15 which represent the one equation for heat transfer
and solid part for the twoequation model.
Figure 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.010 m. 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% 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 zaxis. Therefore, it is necessary to compare the overall thermal performance with the available experimental data.
In the present study, 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:
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).
Table 6: 
Comparison of thermal resistances between present study and
Tuckerman and Pease experimental results 

CONCLUSION
In the preceding analysis, a full threedimensional 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 nondimensionless 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 volumeaveraged 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.
ACKNOWLEDGMENT
The authors gratefully acknowledge the financial support of the Ministry of Science, Technology and Innovation, Malaysia under the project 030102SF0318.