Abstract

Silicon-based embedded microchannel with three-dimensional (3D) manifold (MF) μ-cooler offers lower pressure drop and increased heat removal capability (>1 kW/cm2) for microprocessors and power electronics cooling using single-phase water. In this paper, we present a thermal–fluidic numerical analysis of silicon-embedded microchannel cooling. We develop a full-scale computational fluid dynamics (CFD) model of a large footprint (24 × 24 mm2) device having embedded microchannels and a 3D manifold. It is found that the pressure/velocity distributions at three different critical regions inside the inlet manifold have a significant impact on the temperature distribution. A previous study reported a shift of the chip temperature hot-spot at high flow rates; this study delves deep into the flow and pressure variations within the MF and cold plate (CP) that leads to this shift. This study also investigates the degree of flow maldistribution, first between the manifold channels caused by the plenum and then between the cold plate channels caused by individual MF channels. Finally, this study concludes with a comparison between two different 3D manifold inlet channel heights. The comparison reveals that the manifold with 1.5 mm thickness can reduce the pressure drop by a factor of 4 while maintaining the same thermal resistance of 0.04 K cm2/W, thus indicating an increase in the coefficient of performance (COP) by a factor of 4, compared with a manifold thickness of 0.7 mm.

Introduction

With the increasing power densities of electronic devices, active liquid cooling is regarded as an efficient cooling technique that can deal with the high heat flux beyond 1 kW/cm2 [13]. Microchannel liquid cooling pioneered by Tuckerman and Pease [3] shows that it can remove heat flux of 790 W/cm2 over a 10 × 10 mm2 heated area. However, the major drawback of conventional microchannel heat sinks in Fig. 1(a) is the large pressure drop along the straight microchannels. Higher heat transfer coefficient h requires a small hydraulic diameter Dh, but at a cost of higher pressure drop ΔP (h ∼ 1/Dh and ΔP ∼ L/Dh4). As a solution to this issue, embedded microchannels with three-dimensional (3D) manifold cooler (EMMC) show great potential to reduce the pressure drop and improve the total thermal performance [4]. The EMMC schematic in Fig. 1(b) shows that an added manifold (MF) with alternating inlets and outlets can route the fluid vertically by delivering the cold liquid from the top and collecting the hot fluid through the adjacent outlet channel. The fluid elements only travel a short distance L1 in the cold plate (CP) channel with small hydraulic diameter as compared to the much longer distance L in traditional coolers.

Fig. 1
Microchannel cooler concept (a) conventional two-dimensional microchannel cooler with fluid travel length L and (b) embedded microchannel cooler with 3D manifold with reduced fluid travel length L1
Fig. 1
Microchannel cooler concept (a) conventional two-dimensional microchannel cooler with fluid travel length L and (b) embedded microchannel cooler with 3D manifold with reduced fluid travel length L1
Close modal

In order to predict thermofluidic behavior of an EMMC geometry, empirical correlations and computational fluid dynamics (CFD) are commonly used. Empirical correlations based on the Nusselt and friction factor correlation can provide the fast prediction of the thermal/hydraulic performance of the microchannel cooler. However, simple correlations developed for straight channels and simplified U and L-bends are not able to capture the thermal performance of these EMMCs due to the complexity of the flow delivery between the manifold and the cold plate channels inside the cooler [5]. Therefore, numerical methods based on CFD are preferred for predicting the thermofluidic behavior of the EMMC heatsink [6,7]. Most studies have opted to break the overall cooler into smaller repeating pieces and perform simulations on these reduced-order models, for example, unit cell analysis [810], single-channel manifold simulations [11], and porous medium model approach [9]. Zhou et al. [8] utilized a unit cell model to analyze the sensitivity of the heat sink performance to material selection and geometric parameters. Moreover, deterministic and probabilistic optimization is obtained by employing the unit-cell model simulation [9]. Mandel et al. [10] developed a reduced-order “2.5-D” CFD modeling approach for single-phase flow and heat transfer in manifold-microchannel heat exchangers, indicating an order-of-magnitude reduced computational cost compared to a full 3D model. However, it is shown that this method is above 30%–40% mean error when the dimensionless length is smaller than 0.1. Several studies [812] utilize a two-level modeling approach using the porous media model and system manifold-level analysis. Sharma et al. [11] used a single channel model with half inlet and outlet manifold slot nozzles for developing hot-spot-targeted embedded microchannel structure optimization methodology. The microchannels were modeled as a fluid-saturated porous medium in the single-channel model. However, the accuracy of these models is highly determined by the simplification of the porous medium model and also the detailed information regarding the flow behavior is missing in this hybrid model approach [8], such as the flow distributions across the cold plate microchannels along a single inlet manifold and the pressure/velocity relations between the manifold level and microchannel levels. Most importantly, the unit cell model and single-channel manifold model are all based on the assumption of the uniform flow distribution across the individual inlet manifold and microchannel. This assumption can result in large errors when the flow nonuniformity across the unit cell and single-channel model is pronounced.

Full-scale CFD analysis including the geometric parameters of the manifold and the entire microchannel heat sink can provide us better understanding of the flow patterns inside the microchannels and also the flow interactions between the manifold level and microchannel level. The local level and system level pressure drop can also be extracted from the full-scale model. Additionally, the full-scale CFD simulation can depict how the device responds to a nonuniform heat flux map—this is especially important in scenarios where several heat-producing chips are packed together tightly to be cooled by a single aggressive cooler. Lastly, the full-scale model can be also used to capture the effect of dimensional inaccuracies and variations that are associated with the etching-based fabrication of such large area coolers in silicon [13]. Jung et al. did extensive numerical investigations on the thermal and hydraulic performance of a small-scale EMMC cooler with a dimension of 25 mm2 [4]. Piazza et al. further investigated the broad-level effects of scaling up of the device design, with up to 400 mm2, and reported a shift of the chip temperature hot spot at high flow rates [14]. However, full-scale large device footprint size (>400 mm2) with small cold plate channels (10–100 μm) results in an exponential increase of the number of the meshing element, since minimal meshing element across the minimal feature size is required to achieve sufficient resolution to capture the flow and thermal physics. It is reported that the meshing elements for a full device level EMMC model with 5 × 5 mm2 are around 50 × 106 mesh elements, with a minimal microchannel feature size of 20 μm [15].

In this paper, we investigate a larger cooler with an area of 600 mm2 by simulating the full geometry with improved computation capabilities. We will perform the one-quarter CFD simulations for a large-scale chip size of 24 × 24 mm2 and 3D manifold of 0.7 mm height, with the minimal feature size of 90 μm. The CFD model geometry with its modeling and meshing methodology will be introduced in section “Large Footprint (24 × 24 mm2) EMMCs Model.” The temperature distributions and pressure/flow patterns variations across the CP microchannel and MF will be systematically analyzed in section “Thermal and Flow Pattern Analysis.” In section “Multi-Level Pressure Analysis,” the pressure distributions across the microchannel and manifold level will be investigated. To provide guidance for the unit cell model and single-channel manifold model, the flow uniformity analysis is conducted in section “Manifold Flow Uniformity Analysis.” In section “Impact of Manifold Thickness,” the impact of the manifold thickness is also investigated, thus providing a guideline for the selection of manifold manufacturing techniques.

Large Footprint (24 × 24 mm2) EMMCs Model

Geometry and Parameters.

The geometry of the EMMC consists of two parts, CP and 3D manifold, as illustrated in Fig. 2(a). For the cold plate, there are inlet and outlet channels alternating with a silicon wall width WCP,wall of 90 μm between them. The cold plate microchannel width WCP,channel is chosen as 90 μm. The height of the cold plate microchannel HCP,channel is fixed as 300 μm, leaving the base thickness of the cold plate at 200 μm. The heat source is directly located on the back of the cold plate, where the microchannels are on the front side. The 3D manifold as shown in Fig. 2(b) provides multiple inlets and outlets that are perpendicular with the embedded cold plate microchannels, which can shorten the fluid flow path and further reduce the pressure drop. The manifold height HMF,in used in this study is 0.7 mm and the base of the manifold is kept as 0.3 mm. There are 28 inlet manifold channels and 29 outlet manifold channels in total across the 3D manifold. This study targeted the numerical investigation of the large device footprint area with around 24 × 24 mm2. The geometry parameters used in this study are summarized in Table 1.

Fig. 2
(a) Schematic of the EMMC and fluid routing through the manifold and cold plate and (b) cross section of the manifold with in/out fluid flow. A–A: cold plate microchannel view and B–B: manifold dimensions.
Fig. 2
(a) Schematic of the EMMC and fluid routing through the manifold and cold plate and (b) cross section of the manifold with in/out fluid flow. A–A: cold plate microchannel view and B–B: manifold dimensions.
Close modal
Table 1

Simulated EMMC design dimensions

Manifold designHCP(μm)WCP,channel (μm)WCP,wall (μm)HCP,channel (μm)HMF(μm)HMF,in (μm)WMF,in (μm)WMF,wall (μm)WMF,out (μm)
D150090903001000700250200250
Manifold designHCP(μm)WCP,channel (μm)WCP,wall (μm)HCP,channel (μm)HMF(μm)HMF,in (μm)WMF,in (μm)WMF,wall (μm)WMF,out (μm)
D150090903001000700250200250

CP = cold plate and MF = manifold; dimensions are defined in Fig. 2.

Modeling Methodology.

As illustrated in Fig. 3(a), the two inlet flows coming from the x and x direction, and the inlet flows are convergent into the 28 inlet manifold channels. The outlet flows from 29 outlet manifold channels are finally collected normally (y direction) through a larger outlet plenum located on top of the outlet manifold channels. It is shown that the full EMMC device is symmetric in the X and Z directions. Therefore, a quarter geometry with symmetric boundary conditions is considered for the simulations as shown in Fig. 3(b). The total size dimension of the quarter cut model including silicon substrate is 35 mm × 35 mm, with a hot-spot size of 12 mm × 12 mm in the middle of the device. The fluid directions are indicated in the single manifold channel shown in Figs. 3(c) and 3(d).

Fig. 3
Large footprint 24 × 24 mm2 EMMC CFD model: (a) top view of the EMMC; (b) quarter cut of the EMMC schematic; (c) side view of the single manifold channel schematic with fluid directions; and (d) top view of the single manifold channel
Fig. 3
Large footprint 24 × 24 mm2 EMMC CFD model: (a) top view of the EMMC; (b) quarter cut of the EMMC schematic; (c) side view of the single manifold channel schematic with fluid directions; and (d) top view of the single manifold channel
Close modal

As the fluid flow passes from inlet manifolds with larger hydraulic diameter to cold plate microchannels with smaller hydraulic diameter, and then back to the outlet manifolds, the flow encounters different flow regimes, including turbulent flow, transition flow, and laminar flow based on the hydraulic diameter. It is reported that the kω shear stress transport (SST) model can capture the transition flow with less than 25% prediction error for complicated channel geometry [16]. Therefore, for the laminar flow with Reynolds number below 2000, the laminar model is used in the simulation. As the Reynolds number is higher than 2000, the kω SST model is implemented. The quarter model is meshed in the meshing module in ansys. “Cutcell” meshing with hexahedral elements is used due to the smaller meshed elements number, faster and well-behaved convergence behavior during CFD simulation [17]. The mesh size independence study with five different numbers of mesh elements, from 4.3 × 106 to 24 × 106 is performed, as shown in the Appendix. Based on the meshing sensitivity study, the minimal bulk mesh size is chosen as 0.01 mm to capture the thermal/hydraulic performance of the EMMC cooler. The details of the optimal meshing are shown in Fig. 4, including the manifold channel and cold plate microchannels. The total generated mesh element number is about 24 × 106. The mesh sizing “growth rate” which controls the maximum allowable rate of growth between adjacent mesh elements is set as 1.05. Moreover, the meshing parameter “Number of cells across the gap” controls the size of the minimum cell and simulation resolution across the thinnest and difficult mesh geometry sections, which is the cold plate channel in our case [17]. This value is varied between 3 and 10. Most importantly, the discretization and meshing of the fluid domain at the boundary layer are very crucial to achieve high accuracy for the flow and thermal conjugate modeling at the fluid/solid interface. As illustrated in Fig. 4(c), the boundary layer is meshed with very fine mesh size and inflated until fully covering the boundary layer thickness [17].

Fig. 4
Meshing details of the EMMC quarter model: (a) top view of the meshing; (b) meshing details of the manifold inlet/outlet fluid domain (solid domain of silicon manifold is not shown in this meshed model); and (c) cold plate microchannel fluid and solid domain
Fig. 4
Meshing details of the EMMC quarter model: (a) top view of the meshing; (b) meshing details of the manifold inlet/outlet fluid domain (solid domain of silicon manifold is not shown in this meshed model); and (c) cold plate microchannel fluid and solid domain
Close modal

The numerical analysis is based on the steady-state CFD modeling in ansysfluent v16.0 and v19.0 software. Constant heat flux with 800 W/cm2 was applied at the top surface of the gold heater. Deionized water at an inlet temperature of 300 K was used as the working fluid. All other solid walls are defined as adiabatic. No-slip and coupled wall boundary conditions are held at the interfaces between the solid domain and the fluid domain. For the one-quarter model, symmetric boundary conditions are applied at the middle of the manifold. Constant inlet velocity boundary condition is applied at the fluid inlet, and zero-gauge pressure is applied at the fluid outlet. The thermal and hydraulic boundary conditions are summarized in Table 2. The material dimensions and properties are listed in Table 3.

Table 2

Thermal/hydraulic boundary conditions

ConditionsValue
Inlet velocity0.13 m/s to 0.36 m/s
Heat fluxes800 W/cm2
Viscous modelLaminar/SST
Fluid inlet temperature300 K
Outlet pressure (gauge)0 Pa
ConditionsValue
Inlet velocity0.13 m/s to 0.36 m/s
Heat fluxes800 W/cm2
Viscous modelLaminar/SST
Fluid inlet temperature300 K
Outlet pressure (gauge)0 Pa
Table 3

Dimensions and material properties

MaterialLayer thicknessThermal conductivity
Au0.2 mm297.73 W/m K
Si (cold plate base)0.2 mm130 W/m K
Water thermal conductivityPiecewise-linear, temperature dependent
MaterialLayer thicknessThermal conductivity
Au0.2 mm297.73 W/m K
Si (cold plate base)0.2 mm130 W/m K
Water thermal conductivityPiecewise-linear, temperature dependent

For the numerical scheme of the CFD modeling, SIMPLE (Semi-Implicit Method for Pressure Linked Equations) scheme is used for the pressure–velocity coupling. The numerical discretization with the least square cell-based method is used for the gradient spatial discretization. QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme is used for the momentum and energy discretization and the standard scheme is used for pressure terms in the governing continue, momentum and energy equation. The investigated inlet velocity varies from 0.13 m/s to 0.36 m/s. For the thermal/hydraulic analysis, the inlet velocity is fixed as 0.13 m/s, and the Reynolds number based on the cold-plate microchannel hydraulic diameter is 1125.

Boundary Conditions.

For the large footprint EMMC CFD quarter model, including the silicon manifold solid domain will increase the computation costs dramatically. As shown in Figs. 5(a) and 5(b), the thermal resistance network is illustrated, which indicates that there is a heat conduction path going through the heater to the cold plate substrate to the silicon manifold. In this section, two full-scale quarter models with silicon manifold and without silicon manifold (adiabatic interface) are compared.

Fig. 5
Temperature comparison of the quarter model with silicon manifold and without silicon manifold (adiabatic manifold surface): (a) thermal resistance network model; (b) heat conduction path through the heater to manifold; (c) local temperature profile comparison; and (d) average temperature comparison under different inlet velocity
Fig. 5
Temperature comparison of the quarter model with silicon manifold and without silicon manifold (adiabatic manifold surface): (a) thermal resistance network model; (b) heat conduction path through the heater to manifold; (c) local temperature profile comparison; and (d) average temperature comparison under different inlet velocity
Close modal

As shown in Fig. 5(c), the heat conduction through the Si manifold reduces the overall heat surface temperature. The local temperature profile with boundary condition of the adiabatic manifold surface shows around 10% temperature difference compared with the model with silicon manifold. This is because the total liquid wetting area will be reduced when the manifold surface is considered adiabatic. Moreover, the average heater temperature between these two models is also plotted and compared as the inlet velocity changes from 0.13 m/s to 0.36 m/s. The comparison in Fig. 5(d) shows the maximum difference of 2.5%. Therefore, the model can be simplified by using adiabatic boundary conditions on the manifold–fluid interface, which can highly reduce the computation cost for full-scale model.

Thermal/Flow Pattern Modeling Analysis

After the quarter model simulation was finished, the temperature/velocity/pressure data were then processed in CFD-post. The residuals for momentum, continuity, and energy are kept below 1 × 10−6 for (x, y, z-velocity), 1 ×10−5 for (continuity), and 1 × 10−8 for (energy) as the convergence criteria.

Thermal and Flow Pattern Analysis.

As illustrated in the temperature distribution map in Fig. 6(a), the hot spot is observed after the entrance regime of the device with the least effective cooling. The maximum temperature in that region is 124 °C under the laminar flow condition. Our previous work reported that the location of the hottest region will shift away from the center of the EMMC (x=0mm) to the entrance region (x=12mm) as the flowrate increases [14]. This is due to the upper level inlet manifold flow expanding to the entrance region of the cold plate microchannels, resulting in a higher velocity and low-temperature region at the entrance. Then, the flow velocity reduces along the inlet manifold channel until reaching the middle of the manifold (x=0mm). Finally, the liquid flow will be forced into the cold plate microchannel, resulting in an impingement effect with a high heat transfer coefficient across the cold plate microchannel. Therefore, the hottest region is observed around the position between x=10mm and x=4mm. As the flowrate increase, the inlet manifold flow will bypass the entrance of the bottom level cold plate microchannels with lower flow volume and high temperature at the entrance region. Flow pattern behaviors inside the manifold and cold plate microchannels are investigated in this paper to further understand this phenomenon of temperature nonuniformity.

Fig. 6
One-quarter modeling result with (a) heat source temperature distribution map and (b) temperature profiles along single manifold (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Fig. 6
One-quarter modeling result with (a) heat source temperature distribution map and (b) temperature profiles along single manifold (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Close modal

As shown in Fig. 6(b), the temperature profiles along with single manifold numbers 1 to 11 almost overlaps identically, indicating that the corresponding flow velocity profile from line 1 to line 11 are similar to each other. This attributes to the superiority of the plenum that is able to uniformly distribute the flow across the different inlet manifold channels. The lower temperatures along line 12 and line 13 are observed due to the heat being dissipated by conduction by the silicon near the edge of the cooler. Figure 7(a) illustrates the relationship between the temperature and channel velocity. The curve shows that the lower channel velocity in the middle of the quarter model (x=10mm to x=4mm) is consistent with the high temperature in this region, while larger velocities seen near the entrance (x=12mm to x=10mm) and device center (x=0) are characterized by a higher cooling effect and thus, lower temperature. Figure 7(b) illustrates the relationship between the in-plane velocity along with the inlet manifold level VMx and out-of-plane velocity in the cold plate microchannel VCPz. It can be seen that whenever the VMx is high, the fluid flow primarily happens along the manifold and less fluid is pushed through the underlying cold plate channels as shown by lowered VCPz in these sections. As the flow approaches, the device center stagnation point VMx drops to zero and the rest of flow is forced down to the central cold plate microchannels. This can be observed in Fig. 7(b), where VCPz is seen to gradually rise as we move to the device center from the hot spot (x=4mm).

Fig. 7
One-quarter modeling result with (a) velocity and temperature profiles along single manifold line 1 and (b) velocity profile along inlet manifold and cold plate microchannel (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Fig. 7
One-quarter modeling result with (a) velocity and temperature profiles along single manifold line 1 and (b) velocity profile along inlet manifold and cold plate microchannel (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Close modal

In order to further understand the flow behaviors inside the manifold and cold plate channels, the flow patterns with velocity vectors are analyzed, as shown in Fig. 8. Three typical regions are identified, which are entrance region (x=12mm), hot-spot region (10mm<x=4mm) and manifold middle region (x=0mm). At the entrance region shown in Fig. 8(a), because of the sudden availability of a short length path to the pressure outlet, a significant amount of fluid tries to attain 0 pressure (outlet through the first set of microchannels). So, the flow velocity is high here. This phenomenon is defined as “expansion flow” at the entrance region. The expansion flow hits the entrance of the cold plate microchannel corner and turns into the microchannel, which results in a higher velocity and recirculated flow at the entrance region. At the hot-spot region, just after the entrance, this expansion flow effect isn't present anymore. In this section, the fluid preferably flows in the wider manifold section which provides a lower pressure drop pathway than pressure drop during flow through the comparatively smaller cold plate microchannels. Therefore, the inlet flow moves forward and bypasses the underlying cold plate channels, resulting in a low velocity inside the microchannels, as shown in Fig. 8(b). An accompanying rise in chip temperature is observed as the flow finally reaches the middle portion of the manifold (x=0mm). As the flow approaches the symmetry surface in the center of the inlet manifold (x=0mm), the velocity along the MF channel drops to 0 and the rest of the flow is forced down inside the CP microchannels, resulting in a high velocity in this region. The confined flow pattern is shown in Fig. 8(c). The flow patterns with velocity vectors are directly linked with the temperature distribution of the chip.

Fig. 8
One-quarter modeling result with flow pattern summary with three typical regions (a) entrance region; (b) hot-spot region; and (c) center of the manifold region (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Fig. 8
One-quarter modeling result with flow pattern summary with three typical regions (a) entrance region; (b) hot-spot region; and (c) center of the manifold region (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Close modal

Multi-Level Pressure Analysis.

In Fig. 9(a), the pressure distributions inside the manifold and microchannel are investigated which dictates the velocity distribution. The pressure profiles along the base of the cold plate channel and inlet manifold level for 14 single manifolds in Fig. 9(b) both show small deviations. At the flow entrance region, the pressure difference between Manifold and cold plate is high resulting in the large flow velocities in the cold plate. An increase in pressure is observed at the center of the manifold (x=0mm) on the manifold side, indicating the formation of a stagnation pressure point. As shown in Figs. 9(b) and 9(c), ΔPMC(defined as the pressure drop from the manifold to the CP base) is between 300 and 400 Pa in the hot-spot region. Based on the Darcy–Weisbach equation, the pressure drop of the cold plate is shown as below:
ΔPCP=fHCP/2DcpρVCPz22
(1)
Dcp=4WCPLCP2(WCP+LCP)
(2)
LCP=WMF
(3)
Fig. 9
One-quarter modeling result with (a) pressure drop distribution map; (b) pressure profiles with different cross section along single manifold; and (c) schematics of the pressure drop inside the cooling system (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Fig. 9
One-quarter modeling result with (a) pressure drop distribution map; (b) pressure profiles with different cross section along single manifold; and (c) schematics of the pressure drop inside the cooling system (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125)
Close modal

Based on the equation, ΔPCP varies linearly with the square of the cold plate feeding velocity VCPz, shown as ΔPCPVCPz2.

Manifold Flow Uniformity Analysis.

In this section, the flow uniformity across different manifold channels will be analyzed by using ansyscfd postanalysis. As illustrated in Fig. 10, the velocity distributions for the cross section of 14 inlet manifold channels and 15 outlet manifold channels in the one-quarter model are plotted. Three different locations along the liquid flow directions are investigated. Plane (1) is the flow entrance region along the manifold (x=12mm), while plane (3) is the manifold middle region (x=0mm). Plane (2) is chosen in the hot-spot region (x=6mm). The results show that the highest velocity is located at the Manifold in the plane (1), with less velocity in the cold plate microchannel shown in Fig. 10(a). As illustrated in Fig. 10(b), the inlet manifold velocity Vmanifold continues to decrease while the velocity inside the microchannel Vcpchannel increases. In the symmetric interface plane (3), the velocity distribution across the cold plate microchannel is relatively higher than in other locations. In the interface plane (3), the impingement jet cooling behaviors can be observed, resulting in high cooling performance on the microchannel bottom surface.

Fig. 10
Velocity and pressure distribution inside the manifold and microchannels: (a) flow patterns for the three interface planes (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125); (b) indications of the three different analyzed interface planes; and (c) flow uniformity distribution: average velocity across the single inlet manifold for different Reynolds number ranging from 1125 to 3115
Fig. 10
Velocity and pressure distribution inside the manifold and microchannels: (a) flow patterns for the three interface planes (Vin = 0.13 m/s, flow rate = 2.2 LPM, Re = 1125); (b) indications of the three different analyzed interface planes; and (c) flow uniformity distribution: average velocity across the single inlet manifold for different Reynolds number ranging from 1125 to 3115
Close modal

The average velocity is taken from the inlet flow entrance from plane 1, as defined in Fig. 10(a). It shows good uniformity of the flow distribution across the individual single-channel manifolds for different flow rates, with Reynolds numbers from 1125 to 3115. Therefore, good temperature uniformity between different cooling manifolds can be achieved due to the uniform flow distributions. This observation indicates that the plenum in the full cooler model indeed distributes the fluid almost equally among the manifold channels; however, each of the manifold channels nonuniformly distributes the fluid among the CP channels. This further means that the full model can be easily represented by a much simpler reduced-order single manifold channel model [18].

Impact of Manifold Thickness.

Three-dimensional manifold EMMC can lower the flow path within the EMMC device by as much as 50% as compared to their conventional microchannels, therefore enabling further reduction in pressure drop and increase in coefficient of performance (COP), at similar levels of thermal performance. Second, the manifold side pressure distribution is very important in determining nonuniformity of flow distribution among CP channels and thus change in temperature map of the chip. In this session, two different 3D manifold inlet channel heights based on the quarter model are compared.

As shown in Fig. 11(b), the manifold pressure drop ΔPm can be reduced by changing the manifold conduit height or width. In this section, the impact of the manifold conduit height based on the one-quarter model will be investigated. Since the critical minimal parameters of the cold plate microchannel are the same, the final meshing element of 0.7 mm and 1.5 mm manifold thicknesses are both approximately 24 × 106.

Fig. 11
Tradeoff chart based on the thermal performance and pressure drop of EMMCs with cooling footprints of 24 × 24 mm2 for two different manifold conduit height: 0.7 mm and 1.5 mm (based on the one-quarter model simulation results)
Fig. 11
Tradeoff chart based on the thermal performance and pressure drop of EMMCs with cooling footprints of 24 × 24 mm2 for two different manifold conduit height: 0.7 mm and 1.5 mm (based on the one-quarter model simulation results)
Close modal
In order to evaluate the thermal performance of the EMMC cooler, normalized thermal resistance is defined as below for the performance benchmarking:
R*=(TavgchipTin)q
(4)

where Tavgchip is the average temperature of the heater surface, Tin is the inlet liquid temperature, q is the heat flux of the heater.

The metric COP can be used to include both aspects of the thermal resistance and pumping power, which is defined as the ratio between the cooling power and the required pumping power. Therefore, a higher COP indicates a better cooler design. The definition is shown as follows:
COP=PowerprovidedtoheaterRequiredpumppower=Maxallowedtempincrease/RthRequiredpumppower
(5)

As shown in Fig. 11, for a fixed flowrate of 17 L/min, the pressure drop can be reduced by a factor of 4 while keeping the same thermal resistance of 0.04 K cm2/W, by increasing the manifold conduit height. The normalized thermal resistance of the silicon base substrate and gold heater is also indicated in the figure, with a total value of 0.021 K cm2/W. This is because the increase of the manifold height can reduce the pressure drop across the manifold.

For a fixed pressure drop of 80 kPa, the average heater temperature can be reduced by a factor of 1.2 as the manifold increases from 0.7 mm to 1.5 mm. For a constant pressure drop situation, designs with a taller manifold height require to be able to accommodate a larger coolant mass flux than its shorter manifold counterpart. Therefore, for a fixed pump power, more liquid can be pushed through the taller manifold cooler, thereby enabling much better thermal performance and superior COP. Alternatively, for a fixed maximum chip temperature condition, the taller manifold would be able to achieve the same cooling performance at a lower pressure drop, often allowing us to get rid of bulky pump components that would otherwise be needed. Current silicon microfabrication methods put a limitation on the maximum inlet manifold height of around 0.8 mm; in order to facilitate larger heights, other techniques like 3D printing can be employed. However, a remaining drawback is that the manifold material would be different Si and thus would have a drastically different thermal expansion coefficient (CTE). Thus, debonding of the MF and CP might occur at elevated working temperatures.

Conclusions

In this paper, we present a thermal–fluidic numerical analysis for silicon embedded microchannel cooling. We develop a conjugate fluid and thermal CFD model for a large footprint (24 × 24 mm2) device with embedded microchannels alongside a 3D manifold. Three regions including the entrance region, hot-spot region and manifold middle region are identified to have a big impact on the temperature distribution across the entire footprint. The expansion flow at the entrance region is found because of the sudden availability of a short length path to the pressure outlet. The expansion flow hits the entrance of the cold plate microchannel corner and turns into the microchannel, which results in a higher velocity and recirculated flow at the entrance region. At the hot-spot region, just after the entrance, small velocity inside the CP microchannels is observed since the expansion flow effect is not present anymore and inlet flow moves forward and bypasses the underlying cold plate channels. A high velocity is observed as the flow approaches the middle of the inlet manifold because the velocity along the MF channel drops to 0 and the rest of the flow is forced down inside the CP microchannels.

The pressure distributions inside the manifold and microchannel are also investigated, which in turn dictates the velocity distribution. The simulations show that the pressure profiles along the base of the cold plate channel and inlet manifold level for single manifolds show small deviations. At the flow entrance region, pressure difference between Manifold and cold plate is high resulting in the large flow velocities in the cold plate. An increase in pressure is observed at the center of the manifold in the manifold side, indicating the formation of a stagnation pressure point. Moreover, this study investigates the degree of flow maldistribution, first between the manifold channels caused by the plenum and then between the cold plate channels caused by the individual MF channels. It shows good uniformity of the flow distribution across the individual single-channel manifolds for different flow rates, with Reynolds numbers from 1125 to 3115. This observation indicates that the plenum in the full cooler model indeed distributes the fluid almost equally among the manifold channels; however, each of the manifold channels nonuniformly distributes the fluid among the CP channels. This further means that the full model can be easily represented by a much simpler reduced-order single manifold channel model.

In the final section, two different 3D manifold inlet channel heights based on the quarter model are compared. This revealed that the manifold with 1.5 mm height can reduce the pressure drop by a factor of 4 while maintaining the same thermal resistance of 0.04 K cm2/W, thus indicating a potential four times increase in COP.

Acknowledgment

The authors acknowledge Ford Motor Company and the National Science Foundation Engineering Research Center for Power Optimization of Electro-Thermal Systems (POETS) for support of this work.

Funding Data

  • This project is supported by Ford Motor Company and the National Science Foundation Engineering Research Center for Power Optimization of Electro-Thermal Systems (POETS) with cooperative Agreement No. EEC-1449548.

Nomenclature

COP =

coefficient of performance

Dh =

hydraulic diameter

HCP,channel =

microchannel height

HCP =

cold plate total height

HMF,in =

inlet manifold height

HMF =

manifold height

L =

fluid travel length

q =

heat flux

R* =

normalized thermal resistance

Re =

Reynolds number

Tavgchip =

average temperature of the heater surface

Tin =

inlet temperature

Vin =

inlet velocity

WCP,channel =

microchannel width

WCP,wall =

wall width between the channels (fin width)

WMF,in =

manifold width (inlet)

WMF,out =

manifold width (outlet)

ΔP =

pressure drop between the inlet and outlet

Appendix: One-Quarter Model: Mesh Sensitivity Study

In this Appendix, the meshing sensitivity analysis of the one-quarter mode is first conducted. Jung et al. [19] performed a mesh size independence study in multiple EMMC designs with various number of mesh elements, while the footprint size is 5 × 5 mm2. Since the number of the mesh elements for the independence highly relies on the size of small features such as the hydraulic diameter of CP microchannels, the investigation range of the minimal mesh size will be referred to the previous work [19]. The relevant meshing setting parameters are listed in Table 4. As suggested by the previous work, cut cell mesh is chosen as the assembly meshing method. Two different sizing growth rates with 1.2 and 1.05 are used in this study.

Table 4

Meshing parameters comparisons

No. 1No. 2No. 3No. 4No. 5
Total mesh number4.3 × 1067.7 × 10616 × 10620 × 10624 × 106
Assembly meshing methodCut cellCut cellCut cellCut cellCut cell
Proximity min size0.01 mm0.01 mm0.01 mm0.01 mm0.01 mm
Sizing growth rate1.21.21.21.21.05
Proximity size function sourcesEdgesEdges + facesEdgesEdges + facesEdges + faces
Inflation optionSmooth transitionSmooth transitionSmooth transition
Inflation max layers555
No. 1No. 2No. 3No. 4No. 5
Total mesh number4.3 × 1067.7 × 10616 × 10620 × 10624 × 106
Assembly meshing methodCut cellCut cellCut cellCut cellCut cell
Proximity min size0.01 mm0.01 mm0.01 mm0.01 mm0.01 mm
Sizing growth rate1.21.21.21.21.05
Proximity size function sourcesEdgesEdges + facesEdgesEdges + facesEdges + faces
Inflation optionSmooth transitionSmooth transitionSmooth transition
Inflation max layers555

The results of mesh size independence study with five different number of mesh elements, from 4.3 × 106 to 24 × 106, are plotted in Fig. 12. There is less than 2% change in the average heater temperature and 0.3% change in the total pressure drop, as the number of the mesh element is higher than 20 M. Based on the tradeoff of the accuracy and the computation cost, the inflation layer mesh with growth ratio of 1.2 with maximum five boundary layers is chosen in this study. For the cross section view of the manifold fluid domain meshing (see Fig. 13), the boundary layer mesh can be observed, where the finer mesh is generated at the bottom of the manifold fluid domain. For the meshing details of the cold plate, there are four mesh cells across the cold plate microchannel with 90 μm width, and also thin boundary mesh is observed between the fluid and solid interface.

Fig. 12
Plot of mesh element independence analysis at Tin = 300 K, Vin = 0.13 m/s, heat flux q = 800 W/cm2 conditions. The number of mesh elements used in the study are from 4.3 M to 24 M.
Fig. 12
Plot of mesh element independence analysis at Tin = 300 K, Vin = 0.13 m/s, heat flux q = 800 W/cm2 conditions. The number of mesh elements used in the study are from 4.3 M to 24 M.
Close modal
Fig. 13
Meshing parameters comparisons
Fig. 13
Meshing parameters comparisons
Close modal

References

1.
Agostini
,
B.
,
Fabbri
,
M.
,
Park
,
J. E.
,
Wojtan
,
L.
,
Thome
,
J. R.
, and
Michel
,
B.
,
2007
, “
State of the Art of High Heat Flux Cooling Technologies
,”
Heat Transfer Eng.
,
28
(
4
), pp.
258
281
.10.1080/01457630601117799
2.
Wei
,
T.
,
2020
, “
All-in-One Design Integrates Microfluidic Cooling Into Electronic Chips
,”
Nature
,
585
(
7824
), pp.
188
189
.10.1038/d41586-020-02503-1
3.
Tuckerman
,
D. B.
, and
Pease
,
R. F. W.
,
1981
, “
High-Performance Heat Sinking for VLSI
,”
IEEE Electron Device Lett.
,
2
(
5
), pp.
126
129
.10.1109/EDL.1981.25367
4.
Jung
,
K. W.
,
Zhou
,
F.
,
Asheghi
,
M.
,
Dede
,
E. M.
, and
Goodson
,
K. E.
,
2019
, “
Experimental Study of Single-Phase Cooling With DI Water in an Embedded Microchannels-3D Manifold Cooler
,” IEEE 21st Electronics Packaging Technology Conference (
EPTC
), Singapore, Dec. 4–6, pp.
164
166
.10.1109/EPTC47984.2019.9026600
5.
Jung
,
K. W.
,
Kharangate
,
C. R.
,
Lee
,
H.
,
Palko
,
J.
,
Zhou
,
F.
,
Asheghi
,
M.
,
Dede
,
E. M.
, and
Goodson
,
K. E.
,
2017
, “
Microchannel Cooling Strategies for High Heat Flux (1 kW/cm2) Power Electronic Applications
,” 16th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (
ITherm
), Orlando, FL, July 27, pp.
98
105
.10.1109/ITHERM.2017.7992457
6.
Ryu
,
J. H.
,
Choi
,
D. H.
, and
Kim
,
S. J.
,
2003
, “
Three-Dimensional Numerical Optimization of a Manifold Microchannel Heat Sink
,”
Int. J. Heat Mass Transfer
,
46
(
9
), pp.
1553
1562
.10.1016/S0017-9310(02)00443-X
7.
Boteler
,
L.
,
Jankowski
,
N.
,
McCluskey
,
P.
, and
Morgan
,
B.
,
2012
, “
Numerical Investigation and Sensitivity Analysis of Manifold Microchannel Coolers
,”
Int. J. Heat Mass Transfer
,
55
(
25–26
), pp.
7698
7708
.10.1016/j.ijheatmasstransfer.2012.07.073
8.
Zhou
,
F.
,
Dede
,
E. M.
, and
Joshi
,
S. N.
,
2015
, “
A Novel Design of Hybrid Slot Jet and Mini-Channel Cold Plate for Electronics Cooling
,” 31st Thermal Measurement, Modeling and Management Symposium (
SEMI-THERM
), San Jose, CA, Mar. 15–19, pp.
60
67
.10.1109/SEMITHERM.2015.7100141
9.
Sarangi
,
S.
,
Bodla
,
K. K.
,
Garimella
,
S. V.
, and
Murthy
,
J. Y.
,
2014
, “
Manifold Microchannel Heat Sink Design Using Optimization Under Uncertainty
,”
Int. J. Heat Mass Transfer
,
69
, pp.
92
105
.10.1016/j.ijheatmasstransfer.2013.09.067
10.
Mandel
,
R.
,
Shooshtari
,
A.
, and
Ohadi
,
M.
,
2018
, “
A ‘2.5-D’ Modeling Approach for Single-Phase Flow and Heat Transfer in Manifold Microchannels
,”
Int. J. Heat Mass Transfer
,
126
(Part A), pp.
317
330
.10.1016/j.ijheatmasstransfer.2018.04.145
11.
Sharma
,
C. S.
,
Tiwari
,
M. K.
,
Zimmermann
,
S.
,
Brunschwiler
,
T.
,
Schlottig
,
G.
,
Michel
,
B.
, and
Poulikakos
,
D.
,
2015
, “
Energy Efficient Hot-Spot-Targeted Embedded Liquid Cooling of Electronics
,”
Appl. Energy
,
138
, pp.
414
422
.10.1016/j.apenergy.2014.10.068
12.
Escher
,
W.
,
Michel
,
B.
, and
Poulikakos
,
D.
,
2010
, “
A Novel High Performance, Ultra Thin Heat Sink for Electronics
,”
Int. J. Heat Fluid Flow
,
31
(
4
), pp.
586
598
.10.1016/j.ijheatfluidflow.2010.03.001
13.
Hazra
,
S.
,
Piazza
,
A.
,
Jung
,
K. W.
,
Asheghi
,
M.
,
Gupta
,
M. P.
,
Jih
,
E.
,
Degner
,
M.
, and
Goodson
,
K. E.
,
2020
, “
Microfabrication Challenges for Silicon-Based Large Area (>500 mm2) 3D-Manifolded Embedded Microcooler Devices for High Heat Flux Removal
,” 19th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (
ITherm
), Orlando, FL, July 21–23,
pp.
83
90
.10.1109/ITherm45881.2020.9190541
14.
Piazza
,
A.
,
Hazra
,
S.
,
Jung
,
K. W.
,
Degner
,
M.
,
Gupta
,
M. P.
,
Jih
,
E.
,
Asheghi
,
M.
, and
Goodson
,
K. E.
,
2020
, “
Considerations and Challenges for Large Area Embedded Micro-Channels With 3D Manifold in High Heat Flux Power Electronics Applications
,” 19th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (
ITherm
), Orlando, FL, July 21–23, pp.
77
82
.10.1109/ITherm45881.2020.9190179
15.
Jung
,
K. W.
,
Hazra
,
S.
,
Kwon
,
H.
,
Piazza
,
A.
,
Jih
,
E.
,
Asheghi
,
M.
,
Gupta
,
M. P.
,
Degner
,
M.
, and
Goodson
,
K. E.
,
2019
, “
Parametric Study of Silicon-Based Embedded Microchannels With 3D Manifold Coolers (EMMC) for High Heat Flux (∼1 kW/cm2) Power Electronics Cooling
,”
ASME
Paper No. IPACK2019-6472.10.1115/IPACK2019-6472
16.
Wei
,
T.
,
Oprins
,
H.
,
Cherman
,
V.
,
Beyne
,
E.
, and
Baelmans
,
M.
,
2019
, “
Conjugate Heat Transfer and Fluid Flow Modeling for Liquid Microjet Impingement Cooling With Alternating Feeding and Draining Channels
,”
Fluids
,
4
(
3
), p.
145
.10.3390/fluids4030145
17.
ANSYS Tutorial,
2022, “
Near-Wall Mesh Guidelines
,” ANSYS Tutorial, Canonsburg, PA, accessed Sept. 13, 2022, https://www.afs.enea.it/project/neptunius/docs/fluent/html/ug/node410.htm
18.
Hazra
,
S.
,
Wei
,
T. W.
,
Lin
,
Y. J.
,
Asheghi
,
M.
,
Goodson
,
K. E.
,
Gupta
,
M. P.
, and
Degner
,
M.
, “
Parametric Design Analysis of a Multi-Level 3D Manifolded Microchannel Cooler Via Reduced Order Numerical Modeling
,”
Int. J. Heat Mass Transfer
, 197, p. 123356.10.1016/j.ijheatmasstransfer.2022.123356
19.
Jung
,
K. W.
,
Kharangate
,
C. R.
,
Lee
,
H.
,
Palko
,
J.
,
Zhou
,
F.
,
Asheghi
,
M.
,
Dede
,
E. M.
, and
Goodson
,
K. E.
,
2019
, “
Embedded Cooling With 3D Manifold for Vehicle Power Electronics Application: Single-Phase Thermal-Fluid Performance
,”
Int. J. Heat Mass Transfer
,
130
, pp.
1108
1119
.10.1016/j.ijheatmasstransfer.2018.10.108