Numerical Study of Large Footprint (24 × 24 mm²) Silicon-Based Embedded Microchannel Three-Dimensional Manifold Coolers

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/cm² [1–3]. Microchannel liquid cooling pioneered by Tuckerman and Pease [3] shows that it can remove heat flux of 790 W/cm² over a 10 $×$ 10 mm² 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 D_h, but at a cost of higher pressure drop $ΔP$ (h ∼ 1/D_h and $ΔP$ ∼ L/D_h⁴). 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.

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 [8–10], 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 [8–12] 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 mm² [4]. Piazza et al. further investigated the broad-level effects of scaling up of the device design, with up to 400 mm², and reported a shift of the chip temperature hot spot at high flow rates [14]. However, full-scale large device footprint size (>400 mm²) 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 mm² are around 50 × 10⁶ 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 mm² 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 mm² 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 mm²) 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.

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 mm². The geometry parameters used in this study are summarized in Table 1.

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).

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 × 10⁶ to 24 × 10⁶ 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 × 10⁶. 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].

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/cm² 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.

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.

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.

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.

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⁻⁶ for (x, y, z-velocity), 1 ×10⁻⁵ for (continuity), and 1 × 10⁻⁸ for (energy) as the convergence criteria.

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=0 mm$⁠) to the entrance region (⁠$x=−12 mm$⁠) 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=0 mm$⁠). 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=−10 mm$ and $x=−4 mm$⁠. 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.

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=−10 mm$ to $x=−4 mm$⁠) is consistent with the high temperature in this region, while larger velocities seen near the entrance (⁠$x=−12 mm$ to $x=−10 mm$⁠) 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=−4 mm)$⁠.

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=−12 mm$⁠), hot-spot region (⁠$−10 mm<x=−4 mm$⁠) and manifold middle region (⁠$x=0 mm$⁠). 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=0 mm$⁠). As the flow approaches the symmetry surface in the center of the inlet manifold (⁠$x=0 mm$⁠), 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.

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

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=−12 mm$⁠), while plane (3) is the manifold middle region (⁠$x=0 mm$⁠). Plane (2) is chosen in the hot-spot region (⁠$x=−6 mm$⁠). 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 $Vcp−channel$ 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.

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].

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 × 10⁶.

where $Tavg−chip$ is the average temperature of the heater surface, $Tin$ is the inlet liquid temperature, $q$ is the heat flux of the heater.

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 cm²/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 cm²/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.

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 mm²) 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 cm²/W, thus indicating a potential four times increase in COP.

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.

• 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.

COP =

coefficient of performance

D_h =

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

$Tavg−chip$ =

average temperature of the heater surface

$Tin$ =

inlet temperature

V_in =

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

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 mm². 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.

The results of mesh size independence study with five different number of mesh elements, from 4.3 × 10⁶ to 24 × 10⁶, 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.