## Abstract

The industry shift to multicore microprocessor architecture will likely cause higher temperature nonuniformity on chip surfaces, exacerbating the problem of chip reliability and lifespan. While advanced cooling technologies like two phase embedded cooling exist, the technological risks of such solutions make conventional cooling technologies more desirable. One such solution is remote cooling with heatsinks with sequential conduction resistance from chip to module. The objective of this work is to numerically demonstrate a novel concept to remotely cool chips with hotspots and maximize chip temperature uniformity using an optimized flow distribution under constrained geometric parameters for the heatsink. The optimally distributed flow conditions presented here are intended to maximize the heat transfer from a nonuniform chip power map by actively directing flow to a hotspot region. The hotspot-targeted parallel microchannel liquid cooling design is evaluated against a baseline uniform flow conventional liquid cooling design for the industry pressure drop limit of approximately 20 kPa. For an average steady-state heat flux of 145 $W/cm2$ on core areas (hotspots) and 18 $W/cm2$ on the remaining chip area (background), the chip temperature uniformity is improved by 10%. Moreover, the heatsink design has improved chip temperature uniformity without a need for any additional system level complexity, which also reduces reliability risks.

## 1 Introduction

On-chip hotspots, which are areas of elevated temperature that occur on the chip level, have intense adverse effects on electronic device performance and reliability. New techniques and remedies have been innovated by thermal researchers to alleviate these effects [1]. Much research has focused on microfabricated thermoelectric coolers (TEC), orthotropic thermal interface materials (TIMs), and anisotropic heat spreaders for hotspot remediation [2–6].

The hotspot remedy considered in this study is a modification of traditional microchannel liquid-cooled heatsinks, which adapts them to cooling nonuniform heat fluxes. This modification to heatsink architecture (channel network and/or manifolds) results in conducting more flow toward the hotspots and/or increasing convective heat transfer surface area on the top of the hotspots [7–11]. The former (conducting more flow toward the hotspots) is the core idea of this study. To the best of our knowledge, the work presented in this paper is a precursive study on the design and optimization of liquid-cooled heatsinks that prioritizes hotspots over background for the flow distribution. A few designs were found in literature that target areas of high heat flux by coolant flow rate [12–14]. A microgap heatsink that operated with two separate unmixed fluids was studied by Green et al. [12]. In their design, one fluid is acted as a fluidic spreader dedicated to cooling the hotspot, while the second fluid served as both a coolant for the background heat flux and an on-chip regenerator for the hotspot fluid. Sharma et al. [13] designed a manifold microchannel liquid-cooled heatsink to target on-chip hotspots by conducting the inlet flow to the top of the hotspots. In their design, coolant impinged on the top of the hotspots after leaving the inlet manifold and was collected from the background to the outlet manifold. Kumar and Singh [14] managed chip temperature nonuniformity by taking advantage of flow maldistribution created by different types of heatsink manifolds. Using numerical modeling, they showed that different manifolds (I-Type, Z-Type, C-Type, and N-Type) could reduce the thermal resistance of heatsinks cooling chips with center, right, left, and dual streamlines heating.

A scientific optimization methodology based on design of experiments (DOE) was implemented to optimize this study's design. DOE is a statistical technique that examines the relationship between design and response parameters [15]. When combined with regression analysis and optimization algorithms, DOE is used extensively for optimizing the thermohydraulic performance of cooling devices such as gas and liquid-cooled heatsinks. Rao et al. [16] used DOE response surface method (RSM) and computational fluid dynamics (CFD) to calculate the regression functions for thermal resistance and pump power of a microchannel heatsink. Then, using an iterative JAYA algorithm, they calculated the optimal design parameters. In order to minimize the thermal resistance and pump power of a double-layered microchannel heatsink, Kulkarni et al. [17] performed a multi-objective optimization using DOE RSM and numerical simulation. Naqiuddin et al. [18] employed DOE Taguchi–Grey method to maximize the heat transfer performance and minimize the pressure drop and temperature gradient in a microchannel heatsink. In another study, RSM was utilized for the shape optimization of a microchannel heatsink employed to cool a multidie module (discrete nonuniform heat flux boundary condition), which maximized the heatsink thermal performance and minimized its pressure drop [19]. Using DOE, microchannel heatsinks with different wall surfaces including uncoated, hydrophobic, and super-hydrophobic, were optimized geometrically to minimize thermal resistance and pump power [20]. Nakhchi [21] optimized the heat transfer and pressure drop inside a liquid-cooled sinusoidal wavy channel. The author used DOE response surface method to calculate the optimal geometries that possessed minimum pressure drop and maximum Nusselt number. Chen et al. [22] used DOE Taguchi method to perform a screening experimental design that determined the crucial parameters in the design of a parallel-plain fin heatsink. Using DOE response surface method, they optimized the thermohydraulic performance of the heatsink.

Modern remedies for hotspots, such as TECs, anisotropic graphite heat spreaders, and orthotropic TIMs benefit high-tech designs in material science. However, there are serious limitations and drawbacks associated with them. TECs suffer from high cost and low thermal efficiency. Graphite heat spreaders offer low mechanical strength and flexibility. Orthotropic TIMs have operational temperature limits that make them inappropriate for high heat flux applications. The disadvantages of modern hotspot solutions were the motivation behind this study.

The objective of this work is to numerically demonstrate a novel concept to remotely cool chips with hotspots and to maximize chip temperature uniformity using an optimized flow distribution under constrained geometric parameters for the heatsink and constant total coolant flowrate. A system of manifolds was designed to concentrate the coolant flowrate to the channels over a hotspot with a heat flux eight times that of background. The hotspot, located on the center of the chip, yielded a heat flux of 145.5 $W/cm2$ compared to the background heat flux of 18.1 $W/cm2$. The effect of flow concentration toward the hotspot on the thermal and heat spreading resistivity of the heatsink, chip temperature nonuniformity, and pump power was investigated. Next, using DOE with a full factorial design (FFD), the distribution was numerically optimized to achieve a minimum chip temperature nonuniformity, heatsink thermal resistivity, and pump power. The obtained optimal designs were compared with a baseline uniform flow conventional liquid cooling design.

The studies available in the literatures on hotspot mitigation deal with various techniques to direct the coldest fluid to the hottest region at high flow rates or the cooling solution is embedded on to the die targeting the hotspot. Lee et al. [23,24] presented a modified fin structure that could be the closest study we could observe on modified microchannel design. We have sincerely made a literature search and to the best of our knowledge we would like to affirm that this is the foremost study on split flow multizone microchannels for hotspot mitigation.

## 2 Problem Description

An exploded schematic view of the heatsink design is shown in Fig. 1. As shown in the figure, the chip sits at the bottom of the heatsink with a hotspot located at its center. A thin layer of TIM is placed between the chip and the heatsink. There are two inlet flows of coolant at the inlet manifold of the heatsink ($QI$ and $QII$). Each of them is supplied by separate pumps. That is because the channels in this heatsink design are classified as middle class (on the hotspot) (class I) and side class (class II) channels. This design enables different flows of the coolant to be channeled to each of these channel classes. In this case, $QI$ is the flow supplied to the middle channels over the hotspot and $QII$ is the flow directed towards the side channels. It is worth mentioning that both $QI$ and $QII$ flow rates and the number of channels in each class can be controlled based on requirements. Similarly, two ports in the outlet manifold direct the coolant out of the system. Lastly, a cover sits atop the entire heatsink.

Limitations of the design are categorized as geometrical and operational. The server designers who determine the available space in a server and a heatsink location dictate the geometrical limitations of the heatsink design, which are considered fixed and unchangeable during the optimization process. The geometrical parameters of the module are shown in Fig. 2 and provided in Table 1. Additionally, server designers determine the operational conditions that must be considered in the heatsink design and optimization process. The operational parameters are provided in Table 2. As seen in the table, the hotspot area is less than 3% of the chip area, while the hotspot heat flux is almost eight times that of the background. It is worth mentioning that, in this study, the values of the geometric and operational fixed parameters (fin size, coolant flow rate, hotspot to background heat flux ratio, etc.) are selected based on the previous researches on liquid-cooled heatsinks used for high power applications [16,25–27].

Symbol | Definition | Value (mm) |
---|---|---|

$hB$ | Base height | 3.00 |

$hF$ | Fin height | 4.00 |

$hT$ | TIM height | 0.10 |

$lChip$ | Chip length | 43.40 |

$lF$ | Fin length | 61.80 |

$lH$ | Hotspot length | 6.50 |

$tF$ | Fin thickness | 0.18 |

$wCN$ | Channel network width | 38.20 |

$wC$ | Channel width | 0.20 |

$wChip$ | Chip width | 26.80 |

$wH$ | Hotspot width | 4.76 |

Symbol | Definition | Value (mm) |
---|---|---|

$hB$ | Base height | 3.00 |

$hF$ | Fin height | 4.00 |

$hT$ | TIM height | 0.10 |

$lChip$ | Chip length | 43.40 |

$lF$ | Fin length | 61.80 |

$lH$ | Hotspot length | 6.50 |

$tF$ | Fin thickness | 0.18 |

$wCN$ | Channel network width | 38.20 |

$wC$ | Channel width | 0.20 |

$wChip$ | Chip width | 26.80 |

$wH$ | Hotspot width | 4.76 |

Symbol | Definition | Value |
---|---|---|

T_{In} | Coolant inlet temperature | 300 ($K$) |

$q\u02d9Chip$ | Chip total power | 250 ($W$) |

$q\u02d9HS$ | Hotspot power | 45 ($W$) |

$q\u02d9BG$ | Background power | 205 ($W$) |

$AChip$ | Chip area | 1163.12 ($mm2$) |

$AHS$ | Hotspot area | 30.94 ($mm2$) |

$AHS/AChip$ | Area ratio | 2.66% |

$qHS\u2033$ | Hotspot heat flux | 145.5 ($W/cm2$) |

$qBG\u2033$ | Background heat flux | 18.1 ($W/cm2$) |

$qHS\u2033/qBG\u2033$ | Heat flux ratio | 8.03 |

$Q$ | Coolant flow rate | 33.33 ($cm2/s$) |

Symbol | Definition | Value |
---|---|---|

T_{In} | Coolant inlet temperature | 300 ($K$) |

$q\u02d9Chip$ | Chip total power | 250 ($W$) |

$q\u02d9HS$ | Hotspot power | 45 ($W$) |

$q\u02d9BG$ | Background power | 205 ($W$) |

$AChip$ | Chip area | 1163.12 ($mm2$) |

$AHS$ | Hotspot area | 30.94 ($mm2$) |

$AHS/AChip$ | Area ratio | 2.66% |

$qHS\u2033$ | Hotspot heat flux | 145.5 ($W/cm2$) |

$qBG\u2033$ | Background heat flux | 18.1 ($W/cm2$) |

$qHS\u2033/qBG\u2033$ | Heat flux ratio | 8.03 |

$Q$ | Coolant flow rate | 33.33 ($cm2/s$) |

The hotspot creates a drastically nonuniform temperature distribution on the chip and heatsink base, which reduces their reliability and lifespan.

## 3 Optimization Design Parameters

in which, $Q$ is the total heatsink coolant flow rate (a constant value of $2\u2009lit/min$) and $QI$ is the portion of this flow rate that travels through the $NI$ channels. The number of channels in class II and the coolant flowrate through them depend on the dimensionless values of $\lambda I$ and $\omega I$. Values of $\lambda I$ range from 0.2 to 0.5 ($0.2\u2264\lambda I\u22640.5$), and $\omega I$ can take values ranging between 0.5 and 0.8 ($0.5\u2264\omega I\u22640.8$). A uniform flow is achieved when $\lambda I$ is at its maximum ($\lambda I=0.5$) and $\omega I$ is at its minimum ($\omega I=0.5$), whereas a maximum flow concentration on the hotspot is achieved when $\lambda I$ is at its minimum ($\lambda I=0.2$) and $\omega I$ is at its maximum ($\omega I=0.8$).

## 4 Response Parameters

where $TBAve$ denotes the heatsink average base temperature. With the presence of hotspots in a module, removing the excessive heat and temperature distribution nonuniformity are of higher priority.

in which, $TChipMax$ and $TChipMin$ represent the maximum and minimum chip temperatures, respectively. To achieve a close-to-uniform temperature distribution profile, it is desirable that $\psi $ be close to unity. A perfectly uniform temperature distribution is achieved when $\psi =1$, which means a completely uniform on-chip temperature distribution is obtained (due to $TChipMax=TChipMin=TChipAve$).

where $W\u02d9P,I$ and $W\u02d9P,II$ are the powers for pump I and II, respectively, and $\Delta PI$ and $\Delta PII$ are the average pressure drops across channels in class I and class II, respectively.

## 5 Numerical Analysis

### 5.1 Basic Assumptions.

The following assumptions are considered in the numerical modeling:

steady, incompressible, three-dimensional flow;

fully laminar flow in all channels ($ReC<1000$);

gravity is negligible;

constant thermophysical properties of both solid and fluid phases due to minute temperature variations;

negligible viscous heat generation;

negligible radiation heat transfer.

### 5.2 Governing Equations.

The following are the governing equations for a three-dimensional conduction–convection heat transfer problem in an incompressible, steady laminar flow regime [28,29].

where $V$ is velocity vector field.

In Eq. (11), $\rho f$ is the density of the coolant fluid, $P$ is the pressure field, and $\mu f$ represents the absolute viscosity of the fluid.

In which $kf$ is the conductivity, and $Tf$ is the temperature of the fluid.$\u2009cp,f$ denotes the fluid specific heat capacity at constant pressure.

where $Ts$ is the temperature in the solid phase.

### 5.3 Numerical Domain and Boundary Conditions.

Previous studies were able to reduce the computational domain to one-half of a unit cell of the mini/microchannel array, given the simplicity of a uniform heat flux boundary condition [16,17,30,31]. In this study, the presence of a hotspot introduces a nonuniform heat flux boundary condition. With the geometric symmetry of both the heatsink and the heat flux boundary condition in mind, the only way to conserve computational resources is to use a half symmetry model (Fig. 1).

A schematic of the computational domain and the corresponding boundary conditions are illustrated in Fig. 3. It is important to mention that this figure is not drawn to scale but is drawn with exaggerated dimensions to provide a clear view. The thickness of the fin and width of the channel are significantly larger in the figure than their actual dimensions in the model. In Fig. 3, the coolant enters the channels from the “velocity inlet” and leaves the heatsink from the “pressure constant” surfaces. In order to employ symmetry in the model, zero-gradient boundary conditions are imposed on the symmetry plane, shown as the violet dashed area in Fig. 3. Specific heat fluxes are applied to the chip background and the hotspot. The heat flux value at the hotspot is eight times greater than that of the background to enforce the nonuniform temperature distribution that is typical in such cases. No-slip hydraulic and thermally coupled boundary conditions are applied to the surfaces of the channel walls that are in contact with the flow. An adiabatic boundary condition is applied to all the surrounding surfaces.

### 5.4 Numerical Method.

The system of equations for motion and energy in the computational domain is discretized on a staggered grid. It is performed using the finite volume method (FVM). Second-order upwind and central difference schemes are applied to the spatial discretization of the convective and diffusion terms, respectively. The governing equations are solved by ansysfluent and the field functions (pressure, velocity, and temperature) calculated using the iterative SIMPLEC algorithm [32].

Compared to other discretization methods (finite difference and finite element), the finite volume method is more compatible with the transport equations governing fluid flow and heat transfer. FVM discretizes the governing equations in an integral form by dividing the physical space into several arbitrary polyhedral control volumes (cells). Since FVM is based on the direct discretization of the conservation equations, mass, momentum, and energy are conserved by the numerical scheme that results in high robustness and calculation efficiency [33]. FVM is straightforward to implement on nonuniform/unstructured grids and flexible enough for complex geometries.

A mapped grid, made of hexahedral cells with fine divisions at the TIM layer and near the channel walls, is applied to the computational domain. Regions with finer cell dimensions correspond to regions of higher gradients.

A grid sensitivity analysis is performed based on both temperature nonuniformity of the chip and pump power given by Eqs. (5) and (6), respectively. Figure 4 shows the results of this analysis. As seen in the figure, when the number of cells is greater than 9 × 10^{6}, the pump power and chip temperature nonuniformity variations are less than 1%. As a result, the model is discretized to almost 9 × 10^{6} cells in this study.

### 5.5 Coolant, Metal, and Thermal Interface Material Properties.

There are two distinct solid phases in the numerical simulation, which possess different thermophysical properties. These two phases are TIM and copper. The thermophysical properties of copper at room temperature (300 K) are used. As for the TIM, the material properties of ShinEtsuMicroSi, X23-7783D, which is a commonly used TIM in electronic packaging, is incorporated into the model. Water is chosen as the coolant, with an inlet temperature set to 300 K. It is assumed that the thermophysical properties of the coolant are constant. The thermophysical properties of the coolant and the solid phases are provided in Table 3.

## 6 Hydraulic Verification of the Numerical Simulation

Using the analytical approach below, the hydraulic pump power results across the DOE design points from the numerical simulation are verified. It is important to mention that the analytical solution is obtained for an asymptotic (including both hydrodynamically developing and fully developed) laminar flow in rectangular channels [34].

where $\u03f5$ is the aspect ratio of the channel, which is the ratio of the longer to shorter edges of the channel cross section.

As seen in Fig. 5, the simulation results match very well with the analytical solution for all the cases that are studied in the analysis. Design model numbers (Table 4) in this figure refer to the pump power for each case in ascending order. The maximum and average discrepancies are less than 1.6% and 2.2%, respectively.

$No.$ | $\omega I$ | $\lambda I$ | $NI$ | $NII$ | $QI$$(lit/min)$ | $QII\u2009(lit/min)$ | $ReDhIReDhII$ | $W\u02d9P\u2009$$(W)$ | $RTh\u2033$$(K.cm2/W)$ | $\psi $ |
---|---|---|---|---|---|---|---|---|---|---|

1 | 0.5 | 0.5 | 51 | 50 | 1.0 | 1.0 | 1.0 | 0.22 | 0.625 | 9.202 |

2 | 0.6 | 0.5 | 51 | 50 | 1.2 | 0.8 | 1.5 | 0.23 | 0.614 | 9.113 |

3 | 0.5 | 0.4 | 41 | 60 | 1.0 | 1.0 | 1.5 | 0.23 | 0.613 | 9.105 |

4 | 0.7 | 0.5 | 51 | 50 | 1.4 | 0.6 | 2.3 | 0.26 | 0.608 | 9.007 |

5 | 0.6 | 0.4 | 41 | 60 | 1.2 | 0.8 | 2.2 | 0.26 | 0.605 | 9.002 |

6 | 0.5 | 0.3 | 31 | 70 | 1.0 | 1.0 | 2.3 | 0.26 | 0.604 | 9.001 |

7 | 0.8 | 0.5 | 51 | 50 | 1.6 | 0.4 | 3.9 | 0.30 | 0.602 | 8.845 |

8 | 0.7 | 0.4 | 41 | 60 | 1.4 | 0.6 | 3.4 | 0.30 | 0.600 | 8.870 |

9 | 0.6 | 0.3 | 31 | 70 | 1.2 | 0.8 | 3.4 | 0.32 | 0.600 | 8.890 |

10 | 0.5 | 0.2 | 21 | 80 | 1.0 | 1.0 | 3.8 | 0.34 | 0.602 | 8.917 |

11 | 0.8 | 0.4 | 41 | 60 | 1.6 | 0.4 | 5.9 | 0.37 | 0.599 | 8.653 |

12 | 0.7 | 0.3 | 31 | 70 | 1.4 | 0.6 | 5.3 | 0.39 | 0.599 | 8.738 |

13 | 0.6 | 0.2 | 21 | 80 | 1.2 | 0.8 | 5.7 | 0.44 | 0.602 | 8.806 |

14 | 0.8 | 0.3 | 31 | 70 | 1.6 | 0.4 | 9.0 | 0.49 | 0.602 | 8.472 |

15 | 0.7 | 0.2 | 21 | 80 | 1.4 | 0.6 | 8.9 | 0.56 | 0.607 | 8.644 |

16 | 0.8 | 0.2 | 21 | 80 | 1.6 | 0.4 | 15.2 | 0.72 | 0.619 | 8.339 |

$No.$ | $\omega I$ | $\lambda I$ | $NI$ | $NII$ | $QI$$(lit/min)$ | $QII\u2009(lit/min)$ | $ReDhIReDhII$ | $W\u02d9P\u2009$$(W)$ | $RTh\u2033$$(K.cm2/W)$ | $\psi $ |
---|---|---|---|---|---|---|---|---|---|---|

1 | 0.5 | 0.5 | 51 | 50 | 1.0 | 1.0 | 1.0 | 0.22 | 0.625 | 9.202 |

2 | 0.6 | 0.5 | 51 | 50 | 1.2 | 0.8 | 1.5 | 0.23 | 0.614 | 9.113 |

3 | 0.5 | 0.4 | 41 | 60 | 1.0 | 1.0 | 1.5 | 0.23 | 0.613 | 9.105 |

4 | 0.7 | 0.5 | 51 | 50 | 1.4 | 0.6 | 2.3 | 0.26 | 0.608 | 9.007 |

5 | 0.6 | 0.4 | 41 | 60 | 1.2 | 0.8 | 2.2 | 0.26 | 0.605 | 9.002 |

6 | 0.5 | 0.3 | 31 | 70 | 1.0 | 1.0 | 2.3 | 0.26 | 0.604 | 9.001 |

7 | 0.8 | 0.5 | 51 | 50 | 1.6 | 0.4 | 3.9 | 0.30 | 0.602 | 8.845 |

8 | 0.7 | 0.4 | 41 | 60 | 1.4 | 0.6 | 3.4 | 0.30 | 0.600 | 8.870 |

9 | 0.6 | 0.3 | 31 | 70 | 1.2 | 0.8 | 3.4 | 0.32 | 0.600 | 8.890 |

10 | 0.5 | 0.2 | 21 | 80 | 1.0 | 1.0 | 3.8 | 0.34 | 0.602 | 8.917 |

11 | 0.8 | 0.4 | 41 | 60 | 1.6 | 0.4 | 5.9 | 0.37 | 0.599 | 8.653 |

12 | 0.7 | 0.3 | 31 | 70 | 1.4 | 0.6 | 5.3 | 0.39 | 0.599 | 8.738 |

13 | 0.6 | 0.2 | 21 | 80 | 1.2 | 0.8 | 5.7 | 0.44 | 0.602 | 8.806 |

14 | 0.8 | 0.3 | 31 | 70 | 1.6 | 0.4 | 9.0 | 0.49 | 0.602 | 8.472 |

15 | 0.7 | 0.2 | 21 | 80 | 1.4 | 0.6 | 8.9 | 0.56 | 0.607 | 8.644 |

16 | 0.8 | 0.2 | 21 | 80 | 1.6 | 0.4 | 15.2 | 0.72 | 0.619 | 8.339 |

## 7 Interpretation of Computational Fluid Dynamics Results

The FFD with four levels is applied to the design parameters to generate a total of 16 points in the design space. The total pump power, thermal resistivity, and temperature nonuniformity of the chip are calculated using CFD at the design points.

Table 4 shows the response parameters ($W\u02d9P,\psi ,$ and $RTh\u2033$) and $ReDhI/ReDhII$ at the FFD points. For model No. 1 with a uniform flow distribution ($ReDhI/ReDhII=1$), the minimum pump power is achieved, but both thermal resistivity and temperature nonuniformity are at their peaks. The minimum value of $\psi $ is obtained when pump power and the Reynolds number ratio ($ReDhI/ReDhII=15$) are at their maximum. This corresponds to the design point of model No. 16 in which $80%$ of the flowrate is pumped to $20%$ of the channels on the hotspot (class I). Increasing $ReDhI/ReDhII$ causes the pump power to increase and the temperature nonuniformity to decrease.

Although the maximum thermal resistivity coincides with the maximum temperature nonuniformity and minimum pump power (model No. 1), it reaches a minimum value when pump power and temperature nonuniformity are at intermediate values (model No. 11). When 80% of the flowrate travels through 40% of the channels in the middle (class I), as is the case for model No. 11, the maximum chip/base temperature is at its lowest.

By concentrating the flowrate to the channels above the hotspot (increasing $W\u02d9P$) and deviating from a uniform flow distribution, the heat spreading resistivity ($RSp\u2033$), at the right-hand side of Eq. (21) decreases, which results in a reduction in $RTh\u2033$. This reduction in $RTh\u2033$ initiates at uniform flow distribution of model No. 1 and continues to model No. 11 where $80%$ of the flow travels through $40%$ of the channels on the hotspot. At this point, $RTh\u2033$ reaches its minimum. Concentrating more flow to the channels on the hotspot increases the $RTh\u2033$ and changes its trend from descending to ascending (Fig. 6(b)). The reason is that concentrating most of the flow to the middle channels attenuates the thermal performance of the side channels. This results in an increase in the heatsink $RConv\u2033$. Figure 7 illustrates how the base of the heatsink spreads the heat generated by the chip hotspot and background to the channels with low flowrate at the sides. Using the temperature profile on the plane in the middle of the channels perpendicular to the flow direction, Fig. 8 compares the spread of the heat for two design models, No. 1 and No. 16, corresponding to maximum and minimum $RSp\u2033$ (or $\psi $), respectively. Comparing Figs. 8(a) and 8(b), it is concluded that concentrating the flow on the hotspot improves the spread of the heat through the heatsink base and fins, which, in turn, results in more uniform chip temperature. Figure 9 shows the corresponding temperature profiles of the models No. 1 and No. 16 on a plane close to the root of the fins. These temperature profiles confirm there is an enhanced spread of heat in model No. 16 compared to that of model No. 1.

## 8 Parametric Study

### 8.1 Pump Power.

Figure 10 shows the variation of pump power with the design parameters $\lambda I$ and $\omega I$. As expected, design model No. 1 with a uniform flow distribution possesses the minimum pump power. Increasing $\omega I$ or decreasing $\lambda I$ concentrates the flow to the channels on the hotspot, which results in a higher pressure drop and total pump power.

### 8.2 Chip Temperature Nonuniformity.

Figure 11 illustrates the effect of on-chip flow concentration on the chip temperature nonuniformity and heatsink heat spreading resistivity. The figure shows that both chip temperature nonuniformity and heat spreading resistivity follow the same trend with $\omega I$ at different values of $\lambda I$. Focusing the flow to the channels on the hotspot enhances heat spreading from it and results in a more uniform chip temperature profile.

### 8.3 Thermal Resistivity.

The effect of $\omega I$ and $\lambda I$ on the total heatsink thermal resistivity is demonstrated in Fig. 12. The graph for $\lambda I=0.5$ and $\lambda I=0.4$, which correspond to 50% and 40% of the channels on the hotspot, shows that $RTh$ is reduced by increasing $\omega I$. In these cases ($\lambda I=0.5$ and $\lambda I=0.4$), the low number of side channels (class II), makes it so that the concentration of the flow on the middle channels does not increase heatsink convective thermal resistivity remarkably. Consequently, by increasing $\omega I$ and decreasing $RSp\u2033$ (Eq. (21)), the heatsink thermal resistivity is reduced. In other words, for $\lambda I=0.5$ and $\lambda I=0.4$, $RSp\u2033$ is the only term on the right side of Eq. (21) that depends on $\omega I$. As Fig. 12 shows for $\lambda I=0.3$, the trend of thermal resistivity is almost flat. For $\lambda I=0.3$, when 60% of the channels are placed on the sides, $RSp\u2033$ and$\u2009RConv\u2033$ in Eq. (21) depend on $\omega I$. Increasing $\omega I$ decreases $RSp\u2033$, on the one hand, and increases $RConv\u2033$ on the other hand, which results in nearly constant $RTh\u2033$. This is the reason that when $\lambda I=0.3$, thermal resistvity ($RTh\u2033$) is a weak function of $\omega I$. An ascending trend for thermal resistivity is recognizable at $\lambda I=0.2$. As mentioned in Sec. 7, due to the large number of channels located on the sides (class II), increasing $\omega I$ weakens convective heat transfer and flow bulk motion in these channels which increases heatsink $RConv\u2033$ and results in $RTh\u2033$ growth.

## 9 Optimization

The calculated regression functions are verified both statistically, considering coefficients of determination ($R2$s), and numerically by modeling a few intermediate design points.

The obtained values for $R2$ and $R2(Adjusted)$ for $W\u02d9P$ are 99.8% and 99.7%, respectively. Similar $R2$ values of 99.4% and 99.2% are calculated for $\psi $. Finally, the values of 96.8% and 94.6% are acquired for $R2$ and $R2(Adjusted)$ of $RTh\u2033$.

In addition to the statistical validation, the accuracy of the regression functions in Eqs. (29)–(31) is confirmed numerically. Three intermediate designs are chosen from the available range of design parameters. Table 5 compares the results of the numerical modeling with regression function predictions. As Table 5 shows the discrepancy is less than 1% at all intermediate points for all regression functions.

Design Parameters | $W\u02d9P\u2009(W)$ | $\psi $ | $\u2009RTh\u2033$ ($K.\u2009cm2/W$) | |||||||
---|---|---|---|---|---|---|---|---|---|---|

$\omega I$ | $\lambda I$ | Simulation | Regression | Discrepancy | Simulation | Regression | Discrepancy | Simulation | Regression | Discrepancy |

0.55 | 0.44 | 0.2333 | 0.2335 | 0.09% | 9.0981 | 9.1053 | 0.08% | 0.6125 | 0.6115 | 0.16% |

0.65 | 0.34 | 0.3160 | 0.3157 | 0.09% | 8.8676 | 8.8804 | 0.14% | 0.5996 | 0.5995 | 0.02% |

0.75 | 0.24 | 0.5350 | 0.5372 | 0.40% | 8.5560 | 8.5401 | 0.19% | 0.6049 | 0.6067 | 0.35% |

Design Parameters | $W\u02d9P\u2009(W)$ | $\psi $ | $\u2009RTh\u2033$ ($K.\u2009cm2/W$) | |||||||
---|---|---|---|---|---|---|---|---|---|---|

$\omega I$ | $\lambda I$ | Simulation | Regression | Discrepancy | Simulation | Regression | Discrepancy | Simulation | Regression | Discrepancy |

0.55 | 0.44 | 0.2333 | 0.2335 | 0.09% | 9.0981 | 9.1053 | 0.08% | 0.6125 | 0.6115 | 0.16% |

0.65 | 0.34 | 0.3160 | 0.3157 | 0.09% | 8.8676 | 8.8804 | 0.14% | 0.5996 | 0.5995 | 0.02% |

0.75 | 0.24 | 0.5350 | 0.5372 | 0.40% | 8.5560 | 8.5401 | 0.19% | 0.6049 | 0.6067 | 0.35% |

This study also tried linear regression functions (quadratic forms) for all the response parameters. Still, in the case of pump power and thermal resistivity, they failed to satisfy linearity and homoscedasticity assumptions. To fix the problem, a logarithmic transformation is used for pump power and thermal resistivity (Eqs. (29) and (31)) [35]. A regression function is valid only if it satisfies the standard regression assumptions [36] expressed as:

*Linearity assumption*: The relation between response and design parameters can be estimated by standard multiregression when it is linear in nature.*Normality assumption*: The regression function errors should be normally distributed.*Homoscedasticity (constant variance) assumption*: The variance of the regression function errors should be constant and independent of their mean.

In addition to the above assumptions, the regression functions should be checked for outliers and influential points. A more detailed analysis of regression diagnostics can be found in Ref. [36].

Although logarithmic transformation increases the normality and decreases the effects of outliers, it is mostly used to enhance data linearity and remove heteroscedasticity. A regression diagnostics analysis is performed for the new nonlinear regression functions (Eqs. (29) and (31)) that shows they satisfy standard regression assumptions thoroughly.

In Eqs. (32) and (33), $W\u02d9P$, $\psi $, and $RTh\u2033$ are substituted by Eqs. (29)–(31), respectively. $W\u02d9PMin$, $\psi Min$, and $RTh\u2033\u2009Min$ are the minimum values of the regression functions of Eqs. (29)–(31).$\u2009\xi W\u02d9P$, $\xi \psi $, and $\xi RTh\u2033$ are the weight factors assigned to $W\u02d9P$, $\psi $, and $RTh\u2033$, respectively. Weight factors are chosen according to the designer's priority and the importance they place on the design parameters. The iterative algorithm of JAYA [37] is used for minimizing the objective functions expressed in Eqs. (32) and (33). Figure 13 shows the Pareto curve for the optimization problem of $W\u02d9P\u2212\psi $ associated with the values of thermal resistivity ($RTh\u2033$) at these points. Table 6 indicates the design parameters of these optimal designs. The Pareto curve for the optimization problem of $W\u02d9P\u2212RTh\u2033$ together with the values of $\psi $ at these points are shown in Fig. 14. Table 7 displays $\omega I$ and $\lambda I$ at the obtained optimal points.

$\omega I$ | $\lambda I$ | $W\u02d9P\u2009(W)$ | $\psi $ | $RTh\u2033$ ($K.\u2009cm2/W$) | |
---|---|---|---|---|---|

1 | 0.5 | 0.50 | 0.22 | 9.18 | 0.626 |

2 | 0.5 | 0.40 | 0.23 | 9.10 | 0.612 |

3 | 0.8 | 0.50 | 0.30 | 8.82 | 0.602 |

4 | 0.8 | 0.46 | 0.32 | 8.72 | 0.600 |

5 | 0.8 | 0.42 | 0.35 | 8.69 | 0.599 |

6 | 0.8 | 0.38 | 0.38 | 8.63 | 0.600 |

7 | 0.8 | 0.34 | 0.43 | 8.56 | 0.602 |

8 | 0.8 | 0.31 | 0.48 | 8.51 | 0.604 |

9 | 0.8 | 0.27 | 0.54 | 8.46 | 0.607 |

10 | 0.8 | 0.23 | 0.65 | 8.39 | 0.612 |

11 | 0.8 | 0.20 | 0.73 | 8.34 | 0.616 |

$\omega I$ | $\lambda I$ | $W\u02d9P\u2009(W)$ | $\psi $ | $RTh\u2033$ ($K.\u2009cm2/W$) | |
---|---|---|---|---|---|

1 | 0.5 | 0.50 | 0.22 | 9.18 | 0.626 |

2 | 0.5 | 0.40 | 0.23 | 9.10 | 0.612 |

3 | 0.8 | 0.50 | 0.30 | 8.82 | 0.602 |

4 | 0.8 | 0.46 | 0.32 | 8.72 | 0.600 |

5 | 0.8 | 0.42 | 0.35 | 8.69 | 0.599 |

6 | 0.8 | 0.38 | 0.38 | 8.63 | 0.600 |

7 | 0.8 | 0.34 | 0.43 | 8.56 | 0.602 |

8 | 0.8 | 0.31 | 0.48 | 8.51 | 0.604 |

9 | 0.8 | 0.27 | 0.54 | 8.46 | 0.607 |

10 | 0.8 | 0.23 | 0.65 | 8.39 | 0.612 |

11 | 0.8 | 0.20 | 0.73 | 8.34 | 0.616 |

$\omega I$ | $\lambda I$ | $W\u02d9P\u2009(W)$ | $\psi $ | $RTh\u2033$ ($K.\u2009cm2/W$) | |
---|---|---|---|---|---|

1 | 0.5 | 0.50 | 0.219 | 9.18 | 0.626 |

2 | 0.5 | 0.48 | 0.219 | 9.16 | 0.623 |

3 | 0.5 | 0.44 | 0.222 | 9.13 | 0.617 |

4 | 0.5 | 0.42 | 0.226 | 9.11 | 0.613 |

5 | 0.5 | 0.37 | 0.238 | 9.07 | 0.608 |

6 | 0.5 | 0.36 | 0.243 | 9.06 | 0.607 |

7 | 0.54 | 0.36 | 0.254 | 9.04 | 0.604 |

8 | 0.59 | 0.37 | 0.270 | 8.99 | 0.602 |

9 | 0.65 | 0.37 | 0.294 | 8.92 | 0.600 |

10 | 0.69 | 0.38 | 0.312 | 8.86 | 0.599 |

11 | 0.72 | 0.20 | 0.337 | 8.78 | 0.599 |

$\omega I$ | $\lambda I$ | $W\u02d9P\u2009(W)$ | $\psi $ | $RTh\u2033$ ($K.\u2009cm2/W$) | |
---|---|---|---|---|---|

1 | 0.5 | 0.50 | 0.219 | 9.18 | 0.626 |

2 | 0.5 | 0.48 | 0.219 | 9.16 | 0.623 |

3 | 0.5 | 0.44 | 0.222 | 9.13 | 0.617 |

4 | 0.5 | 0.42 | 0.226 | 9.11 | 0.613 |

5 | 0.5 | 0.37 | 0.238 | 9.07 | 0.608 |

6 | 0.5 | 0.36 | 0.243 | 9.06 | 0.607 |

7 | 0.54 | 0.36 | 0.254 | 9.04 | 0.604 |

8 | 0.59 | 0.37 | 0.270 | 8.99 | 0.602 |

9 | 0.65 | 0.37 | 0.294 | 8.92 | 0.600 |

10 | 0.69 | 0.38 | 0.312 | 8.86 | 0.599 |

11 | 0.72 | 0.20 | 0.337 | 8.78 | 0.599 |

## 10 Sensitivity Analysis

Sensitivity analysis is a mathematical technique used to evaluate how the response parameters are influenced by uncertainties in the design parameters. The analysis can be performed both by taking partial derivatives of regression functions (Eqs. (29)–(31)) with respect to independent (design) parameters [38] or the results of F-values obtained from the analysis of variance (ANOVA) [39,40]. The *F*-value is a non-negative real number defined as the ratio of the mean regression sum of squares to mean error sum of squares. The higher *F*-value that a term has in the analysis of variance, the more significant its effect on the response parameter. Using the results of *F*-values calculated through the analysis of variance for generating regression functions (Eqs. (29)–(31)), a sensitivity analysis was done to estimate the influence of design parameters and their combinations (regression functions terms) on the response parameters. Figure 15 shows the results of this sensitivity analysis. As seen in the figure, quadratic and interaction terms have greater impacts on the response parameters than linear terms. The interaction term $\omega 1\lambda 1$ is the most influential term for this problem. It is attributed to the fact that the flow concentration on the high heat flux area includes the contribution of both design parameters simultaneously. In other words, to create a concentrated flow distribution, one must specify both the number of channels and the flowrate in class I.

## 11 Conclusion

In this study, the authors numerically demonstrated a novel concept for energy efficient, hotspot-targeted liquid cooling of multicore microprocessors. The multizone heatsink consisted of varyingly distributed channels and flow over the hotspots and the background as a potential design to maximize temperature uniformity over the chip surface. Uniform flow distribution over the entire chip surface was taken as the baseline for comparison. Detailed three dimensional CFD simulations were performed for a series of designs and then a surrogate model representing the actual physics was developed. The surrogate model was used for optimization. The optimization study results made it evident that concentrating the flow on the small areas of elevated heat flux improves the spread of heat from these regions and results in a more uniform chip temperature. An improvement of 10% was achieved in the chip temperature uniformity compared to the baseline case. However, the results showed that concentrating the flow on the hotspot decreased the system thermal resistivity up to some point (where $RSp\u2033$ and $RConv\u2033$ are minimum) after which, its effect followed an adversely reverse course. In other words, there was a relative minimum for thermal resistivity with respect to flow concentration on the hotspot.

## Funding Data

NSF Industry/University Cooperative Research Center (Award No. IIP-1738793; Funder ID: 10.13039/100000001).

## Nomenclature

- $A$ =
area ($m2$)

- $cp$ =
specific heat capacity ($J/kg\u2009K$)

- $Dh$ =
hydraulic diameter ($m$)

- $f$ =
- $F$ =
- $h$ =
height ($m$)

- $k$ =
thermal conductivity ($J/m\u2009K$)

- $l$ =
length (m)

- $m\u02d9$ =
mass flow rate ($kg/s$)

- $p$ =
perimeter ($m$)

- $P$ =
pressure ($Pa$)

- $Q$ =
coolant volumetric flow rate ($m3/s$)

- $q\u02d9$ =
power ($W$)

- $q\u2033$ =
heat Flux ($W/m2$)

- $R\u2033$ =
thermal resistivity ($K\u2009\u2009m2/W$)

- $R2$ =
coefficient of determination

- $Re$ =
Reynolds number

- $t$ =
thickness ($m$)

- $T$ =
temperature ($K$)

- $V$ =
velocity ($m/s$)

- $w$ =
width (m)

- $W\u02d9$ =
pump power ($W$)