CFD Study on Micro Channel Heat Sink

Overview

A numerical study on heat transfer and fluid flow characteristics of a silicon micro channel heat sink with circular ribs on the channel walls. The ribs disrupt the thermal boundary layer to enhance convective heat transfer — a critical need as chip power densities push past 500 W in modern AI accelerators and GPUs.

The study simulates 9 Reynolds numbers (100–900) using ANSYS Fluent with conjugate heat transfer on a 2.3M-cell mesh. Results were validated against published literature for both Nusselt number and friction factor. The geometry was modelled in SolidWorks due to ANSYS SpaceClaim's inability to handle the micro-scale rib dimensions (25 µm radius).

The contour comparison plots were built by exporting section plane data from Fluent and stacking them in Matplotlib under a unified colour bar — something that cannot be done natively inside ANSYS.

Geometry & Problem Setup

Micro Channel with Circular Ribs

The computational domain represents a single periodic element of a heatsink with 50 parallel channels. The channel has a rectangular cross section (W = 0.1 mm, H = 0.2 mm) and a length of 20 mm. The solid domain (silicon) surrounds the channel on three sides while the top is modelled as an adiabatic wall.

Symmetry boundary conditions on both lateral walls of the solid domain replicate the full heatsink behaviour, reducing the computational cost to a single-channel slice. A constant heat flux of 1×106W/m21 \times 10^6 \, \text{W/m}^2 is applied to the bottom face.

Rib Arrangement

Circular ribs of 25 µm radius protrude from the vertical channel walls with a pitch of 0.4 mm. The sector angle (60° for this study) controls how far the rib extends into the flow. These features disrupt the developing boundary layer and introduce local velocity increases that enhance mixing and heat transfer.

Assumptions & Simplifications

Steady-State

No temporal variation — heatsink operates at thermal equilibrium with constant heat flux and flow rate.

Laminar Flow

Reynolds numbers 100–900, well below the critical Re = 2300 for internal flow. Laminar viscous model used.

Incompressible

Water as working fluid at low velocities — density treated as constant with no compressibility effects.

Constant Properties

Fluid and solid properties evaluated at inlet conditions. Justified by the small temperature rise in the fluid.

Boundary Layer Sizing

First Cell Height from y⁺ Requirement

The paper requires y⁺ < 1 at the solid-fluid interface. A conservative target of y+=0.8y^+ = 0.8 was used to derive the first boundary layer cell height analytically.

Starting from the y⁺ definition:

y+=yuνy^+ = \frac{y \, u^*}{\nu}

Substituting the friction velocity u=τw/ρu^* = \sqrt{\tau_w / \rho} and approximating wall shear stress with the Darcy friction factor for laminar flow (Cf=64/ReC_f = 64 / Re), the first cell height simplifies to:

y=0.8Dh412Rey = \frac{0.8 \, D_h}{4} \sqrt{\frac{1}{2 \, Re}}

where Dh=2WH/(W+H)=0.1333mmD_h = 2WH / (W + H) = 0.1333 \, \text{mm} is the hydraulic diameter. Using the most restrictive case (Re = 900) gives:

ymin=6.285×107m0.63μmy_\text{min} = 6.285 \times 10^{-7} \, \text{m} \approx 0.63 \, \mu\text{m}

This value was used as the first layer height for all simulations. A separate boundary layer on the solid side with a first height of 1 µm ensures adequate aspect ratio across the interface.

Computational Mesh

Cross Section

Cross-sectional view of the 2.3M-cell mesh showing the solid and fluid domains. The fluid region has significantly finer cells to resolve the velocity and thermal boundary layers, while the solid region uses a coarser mesh graded toward the interface.

Boundary Layers at Interface

Detail of the inflation layers on both sides of the solid-fluid interface. The fluid-side first cell height of 0.63 µm satisfies the y⁺ requirement, while the solid-side 1 µm first layer maintains reasonable aspect ratios across the boundary.

Near the Circular Ribs

Horizontal cut through the mesh at the rib locations. The 25 µm rib features required substantial local refinement to resolve the flow disruption they create. This was one of the main drivers of the overall cell count.

Surface Mesh at Rib

Surface mesh on the solid-fluid boundary near a circular rib. The extremely small boundary layer cells around the rib curvature led to challenging mesh quality metrics — minimum orthogonal quality of 0.25 and maximum skewness of 0.75, both near the acceptable limits.

Mesh Quality Metrics

MetricMeasuredRecommendedAcceptable
Total cells2,295,285
Min element quality0.127≥ 0.50.1–0.5
Max aspect ratio83.2< 10 (bulk)10–50
Min orthogonal quality0.250≥ 0.30.1–0.3
Max skewness0.750< 0.50.5–0.85

Simulation Settings

Numerical Schemes

Convective FluxSecond-order upwind — avoids numerical diffusion near ribs without QUICK oscillations
Diffusive FluxSecond-order central differencing — Fluent default for finite volume gradients
Viscous ModelLaminar — Re 100–900, below critical Re = 2300 for internal flow
Pressure-VelocitySIMPLE — segregated algorithm, lower memory than coupled solver
Pressure InterpolationSecond-order — consistent with other scheme choices
Convergence CriteriaContinuity: 10⁻³, Velocity: 10⁻⁶, Energy: 10⁻⁹

Grid Independence & Convergence

Grid Independence Study

Eight meshes from 908K to 2.3M cells were tested at Re = 100. Both the average Nusselt number and friction factor converge as the cell count increases, with successive variations dropping below 0.2%. The final 2.3M-cell mesh was selected for all subsequent simulations.

Residual Convergence

Residual plot for the Re = 100 case showing continuity, velocity components, and energy residuals dropping to their prescribed thresholds. Mass flow conservation was also monitored — the inlet-outlet mass flow difference converges to near zero, confirming global conservation.

Output Parameters

Nusselt Number & Friction Factor

The two key performance metrics were computed from the simulation data:

Average Nusselt Number

Nu=m˙Cp(TˉoTˉi)DhkAc(TˉcTˉf)\overline{Nu} = \frac{\dot{m} \, C_p \, (\bar{T}_o - \bar{T}_i) \, D_h}{k \, A_c \, (\bar{T}_c - \bar{T}_f)}

Computed from the energy balance of the fluid using mass average temperatures at inlet/outlet and the area-weighted average temperature of the solid-fluid interface.

Average Friction Factor

fˉ=ΔpDh2ρfLUm2\bar{f} = \frac{\Delta p \, D_h}{2 \, \rho_f \, L \, U_m^2}

Calculated from the pressure drop across the channel. Outlet pressure is the zero-gauge BC; inlet pressure is extracted as the area-weighted average.

Section Plane Results

Contour plots for all 9 Reynolds numbers (100–900) were generated by exporting section plane data from ANSYS Fluent and post-processing in Python with Matplotlib. Each subplot shares a unified colour bar so that the variation across Reynolds numbers can be compared directly — something that cannot be achieved natively in ANSYS, where each plot uses its own auto-scaled colour range.

Pressure — Mid Plane

Inlet pressure increases with Reynolds number as expected — while higher flow rates improve heat transfer, they demand proportionally higher pumping power. The outlet remains at zero gauge pressure across all cases.

Temperature — Mid Plane

The solid domain temperature increases along the flow direction as the fluid heats up, reducing the local temperature differential. At low Re, the fluid reaches higher temperatures but the solid also runs hotter; at high Re, efficient heat extraction keeps both domains cooler.

Temperature — Base Plane

Temperature contours at the bottom of the fluid channel (solid-fluid boundary). The distribution mirrors the mid plane trends but reveals the direct contact thermal profile where convective heat transfer occurs.

Velocity — Top Plane

The effect of the ribs is visible at low Reynolds numbers as periodic elliptical regions of increased velocity coinciding with rib locations. As Re increases, these regions elongate and eventually merge around Re ≈ 600. The velocity structure transitions from isolated disturbances to a developed core flow separated from the boundary layers.

Validation & Comparison

Nusselt Number vs. Reynolds Number

Good agreement at low Reynolds numbers, with increasing divergence at higher Re. The discrepancy is attributed to the boundary layer height being derived purely from flow dynamics (y⁺ calculation) without considering the thermal boundary layer thickness — leading to under-resolved thermal gradients at higher flow speeds.

Friction Factor vs. Reynolds Number

The friction factor shows excellent trend agreement with a constant offset from the literature values. This Re-independent offset suggests a systematic mesh-related difference rather than a physics modelling error — consistent with the grid independence study being performed only at Re = 100.

Key Findings

  • Circular ribs enhance heat transfer by disrupting the thermal boundary layer with minimal impact on pumping power.
  • Higher Reynolds numbers reduce both fluid bulk temperature rise and heatsink steady-state temperature, but at the cost of significantly increased pressure drop and pumping power.
  • The laminar viscous model may miss small-scale turbulence effects induced by the ribs, particularly the elliptical high-velocity regions observed in the velocity contours.
  • A more accurate approach for future work would consider the thermal boundary layer when sizing the first cell height, and use surface heat flux integration rather than outlet temperature for Nusselt number calculation.

Technology Stack

ANSYS FluentCFD solver and meshing
SolidWorksGeometry — SpaceClaim failed at 25 µm ribs
Python + MatplotlibContour plots with unified colour bars
ANSYS SpaceClaimNamed selections and volume extraction
Fluent Meshing2.3M-cell mesh with inflation layers

Nuwantha Kumara

Mechanical Engineering student passionate about software development, simulations, and creating impactful solutions.

Connect

© 2026 Nuwantha Kumara. All rights reserved.

Built withNext.js