## Abstract

Three-dimensional (3D) stacked integrated circuit (SIC) chips are one of the most promising technologies to achieve compact, high-performance, and energy-efficient architectures. However, they face a heat dissipation bottleneck due to the increased volumetric heat generation and reduced surface area. Previous work demonstrated that pin-fin enhanced microgap cooling, which provides fluidic cooling between layers could potentially address the heat dissipation challenge. In this paper, a compact multitier pin-fin single-phase liquid cooling model has been established for both steady-state and transient conditions. The model considers heat transfer between layers via pin-fins, as well as the convective heat removal in each tier. Spatially and temporally varying heat flux distribution, or power map, in each tier can be modeled. The cooling fluid can have different pumping power and directions for each tier. The model predictions are compared with detailed simulations using computational fluid dynamics/heat transfer (CFD/HT). The compact model is found to run 120–600 times faster than the CFD/HT model, while providing acceptable accuracy. Actual leakage power estimation is performed in this codesign model, which is an important contribution for codesign of 3D-SICs. For the simulated cases, temperatures could decrease 3% when leakage power estimation is adopted. This model could be used as electrical-thermal codesign tool to optimize thermal management and reduce leakage power.

## 1 Introduction

Micro-electronics have become smaller, and transistor density has doubled every three years since 2013 following the celebrated Moore's law [1]. This rapid development contributes to the heat dissipation bottleneck of two-dimensional stacked integrated circuits (2D-SICs). Thermal challenges in 2D structural ICs appear for several reasons. The size of transistors cannot be miniaturized indefinitely, since it will eventually reach the fundamental physical barrier [2]. That is, when the thickness of a gate oxide layer is less than four layers of silicon atoms, current will penetrate the gate oxide, leading to the failure of the transistor because of large leakage. Fin field-effect transistor (FinFet) has been widely investigated and has contributed to scaling since 2001 [3]. Although 16 nm and 22 nm transistor gates have been accomplished with this technology, the physical limitation of atomic layer thickness still exists, which explains the prediction that transistor size will still decrease with time but will eventually reach a limit [4–6].

In addition to the scaling problem of transistors, interconnects have also become a major bottleneck for high-performance computing systems [7]. As ICs contain more electronic components in the planar direction, such as transistors, the interconnects of these components become longer, which increases both latency and required switching energy. For example, the latency of a 1 mm long interconnect will be 100 times as large as, and its binary switching energy 30 times as large as, the latency and energy of a transistor with 35 nm generation technology [7]. Higher energy consumption from the longer interconnection length contributes to intensifying Joule heating. To sum up, all the above factors contribute to increases in latency and power consumption, and reductions in the expected benefits from scaling.

Overcoming the above problems associated with 2D-SICs, the vertical integration of chips, or three-dimensional (3D)-SICs, could reduce the wire interconnection length by as much as 50% [8]. The shorter wire length reduces the global RC delay and increases the wire-limited clock frequency by 3.9 times [9]. Thus, remarkable power reduction in the interconnection can be expected. In addition, 3D structures also achieve wide bandwidth buses between functional blocks in the chips [10], which could contribute to improved signal quality [11,12]. More importantly, 3D integration realizes a combination of various heterogeneous chips within a single package, including memory, logic, micro-electromechanical systems, optoelectronics, etc. which would increase the functionality of the package with reduced size.

Thermal challenges arising with the 3D-SICs can potentially be handled with active interlayer liquid cooling [13–15]. Pin-fin enhanced microgap cooling has received increasing interest as chip layers are vertically connected by through silicon vias (TSVs) for die-to-die signal communication [16–18]. The pin-fins in the microgap also increase the liquid contact surface area for improved cooling capacity [19–22]. For multiple layers in 3D ICs, simulations have been employed [23,24].

Arvind et al. modeled the fluidic layers, as a porous medium, to solve the inline and staggered pin fin configurations [25,26]. This flexible model is compatible with thermal computer-aided design tools for ICs and offers higher speed-up over commercial computational fluid dynamic (CFD) simulators. However, it does not resolve the effects of pin fin dimensions, and transversal, and longitudinal spacing, on the heat transfer coefficient. The thermal and hydraulic characteristics for different dimensions are not expected to be the same.

Wan et al. [27,28] performed computational fluid dynamic/heat transfer (CFD/HT) modeling for various pin-fin configurations to predict accurate heat transfer coefficient and pressure drops. Correlations were established for the friction factor and Colburn *j* factor for dense arrays of micropins (22< Re <357, pitch/diameter ratios of 1.5–2.25, and height/diameter ratios of 1.5–2.25), based on CFD/HT modeling. The optimized 3D-SICs packaging structure was selected from four alternative packaging organizations by using this codesign model. However, this model is only limited to steady-state simulation and to two chip layers.

High speed thermal modeling is a required tool for electrical-thermal codesign. Serafy et al. [29] codesigned 3D-SICs architectures to optimize the electrical performance and cooling capacity. It was demonstrated that, the codesign approach using compact thermal simulations achieves reduced electrical leakage, better performance and energy efficiency, than optimizing only the cooling capacity.

In this work, a compact transient thermal model for 3D-SICs with pin-fin enhanced microgap is presented. The model considers spatially and temporally varying heat flux distributions in each tier, and various coolant characteristics in each microgap. In addition to accurate thermal simulations based on pin-fin configurations, the rapid execution feature enables the compact model as a tool for electrical and thermal codesign. Section 2 describes the methodology of this thermal compact model and its CFD/HT verification. Applications of the compact model are described in Sec. 3, followed by performance analysis of the results in Sec. 4.

## 2 Compact Thermal Model

Figure 1(a) shows a schematic of 3D-IC microfluidic cooling for four layers, each of which contains several components with different power. Heat is dissipated by a coolant that runs through the microgap between layers connected by cylindrical pin-fins, which transfer both heat and electronic signals, as the pins act as TSVs.

### 2.1 Thermal Analysis.

*j*factor have been utilized to calculate the heat transfer coefficient and pressure drop [27], as indicated in Eqs. (1) and (2) and Table 1.

C | $\alpha 1$ | $\alpha 2$ | $\alpha 3$ | m | ||
---|---|---|---|---|---|---|

f | Re $<$ 100 | 3.1335 | −0.4485 | −0.4965 | −0.5553 | −0.6292 |

Re $>$ 100 | 1.246 | −0.3362 | −0.4478 | −0.4615 | −0.4393 | |

j | Re $<$ 100 | 0.5885 | 0.0072 | −0.1432 | −0.1289 | −0.5697 |

Re $>$ 100 | 0.4481 | −0.1285 | −0.1707 | 0.0804 | −0.4864 |

C | $\alpha 1$ | $\alpha 2$ | $\alpha 3$ | m | ||
---|---|---|---|---|---|---|

f | Re $<$ 100 | 3.1335 | −0.4485 | −0.4965 | −0.5553 | −0.6292 |

Re $>$ 100 | 1.246 | −0.3362 | −0.4478 | −0.4615 | −0.4393 | |

j | Re $<$ 100 | 0.5885 | 0.0072 | −0.1432 | −0.1289 | −0.5697 |

Re $>$ 100 | 0.4481 | −0.1285 | −0.1707 | 0.0804 | −0.4864 |

These correlations have been developed based on detailed CFD/HT simulations within certain ranges (22< Re < 357, 1.5 < ratios of the pitch to diameter <2.25, and 1.5< ratios of height to diameter <2.25), and have been validated [27,28]. The heat transfer coefficient inside the pin-fin enhanced microgap would be in the certain range (23,045 $<$h $<$ $56,731$, if $22<$ Re$<$ 100; 45,686$<$h $<$ $101,690$, if $100<$ Re $<$ 357). In addition to the forced convection in the microgap, the model includes an effective heat transfer coefficient from the 3D-SICs to the ambient, including the net effect of natural convection and radiation.

### 2.2 Verification of the Compact Model.

To validate the compact microfluidic cooling model, a four-layer CFD/HT model was built. The CFD/HT model, Fig. 2(a), simulates one stripe along the way because of its symmetric structure. Symmetric boundary condition has been applied to two sides of the simulated strip, as shown in Fig. 2(a). At top surface of the strip, a heat transfer coefficient of 10 W/m^{2}K has been applied. At the bottom of the strip, a heat transfer coefficient of 550 W/m^{2}K has been adopted as an effective heat transfer coefficient from the chip package to the circuit board via underfills [28]. Silicon is used as the chip material and water is used as the coolant at a constant inlet velocity of 0.3 m/s, opening condition with zero-gage pressure at the outlet. Transient simulation has been performed for 1 min with 1 s time-step. Radiation heat transfer, on the external surfaces, is negligible. Model parameters are listed in Table 2.

T_{amb} | 300 K | T_{in} | 300 K |
---|---|---|---|

h_{top} | 10 W/m^{2 }K | T_{ini} | 300 K |

L | 9400 × 10^{−6} m | h_{bot} | 550 W/m^{2}K |

W | 200 × 10^{−6} m | H_{pins} | 300 × 10^{−6} m |

D_{p} | 100 × 10^{−6} m | S_{lo} | 200 × 10^{−6} m |

q | 2.5 × 10^{10} W/m^{3} | S_{tr} | 200 × 10^{−6} m |

T_{amb} | 300 K | T_{in} | 300 K |
---|---|---|---|

h_{top} | 10 W/m^{2 }K | T_{ini} | 300 K |

L | 9400 × 10^{−6} m | h_{bot} | 550 W/m^{2}K |

W | 200 × 10^{−6} m | H_{pins} | 300 × 10^{−6} m |

D_{p} | 100 × 10^{−6} m | S_{lo} | 200 × 10^{−6} m |

q | 2.5 × 10^{10} W/m^{3} | S_{tr} | 200 × 10^{−6} m |

The dimensions of the layers and boundary conditions from the compact thermal model have been applied to the CFD/HT model with 3.5 × 10^{6} nodes. Temperature variations stay below 0.2 K after nodes reach 3.3 × 10^{6}, which meets the requirement of mesh independence for reliable temperature results. Sweep meshing method [30] has been adopted for the four pin-fin layers, which enables higher quality meshing. Inflation meshing option [30] has been adopted for pin-fin liquid interfaces and set as six maximum layers with 1.1 growth rate. This meshing option increases meshes around the pin-fin to accurately resolve fluid boundary layers. Performance of this inflation meshing option can be seen in partial enlargement of meshing in the side view, Fig. 2(c), which shows that the maximum temperatures from these two models closely agree. In addition to its accuracy, the compact model is also computationally efficient, taking minutes to execute. The CFD/HT model, when running on a 2.5 GHz CPU and 8 GB memory desktop takes around 20 h.

### 2.3 Thermal-Electrical Codesign.

The thermal compact model presented here has been used for leakage prediction. For the Penryn processor investigated in this paper, 80% of chip power is considered as dynamic, and the remaining 20% as leakage by implementing McPat power simulator [31]. Leakage power is correlated with chip temperature, as shown in Fig. 3 by implementing 45 nm high-performance model with 50 inverter circuits in HSPICE [32]. The higher the temperature, the higher the leakage power. The leakage power is estimated based on a 100 °C chip temperature, which can be an over-estimation depending on the actual temperate of the chip [33]. This over-estimated leakage power is corrected by using the realistic power ratio and temperature correlation, as shown in Fig. 3. Realistic power is acquired by deducting the over-estimated leakage power. Based on this updated power, chip temperatures are updated in the compact model by considering both thermal and electrical aspects. The overall calculation procedure is shown in Fig. 4.

## 3 Configurations and Power Maps of the Four Layers Three-Dimensional-Stacked Integrated Circuits Compact Model

The 3D-SICs model with pin-fin enhanced microgaps explored here is shown in Fig. 5. It is a four active layers model, with each layer having a pin-fin enhanced microgap on top. The thickness of microgaps is the height of the staggered pin-fins, as all four layers are stacked together and connected through pin-fins. De-ionized water as a coolant is individually pumped through the four microgaps. Parameters set in the model are listed in Table 3.

T_{amb} | 293 K | T_{in} | 298 K |
---|---|---|---|

h_{top} | 10 W/m^{2 }K | T_{ini} | 298 K |

L | 12000 × 10^{−6} m | h_{bot} | 550 W/m^{2}K |

d_{Si} | 100 × 10^{−6} m | W | 16000 × 10^{−6} m |

d_{SiO2} | 10 × 10^{−6} m | H_{pins} | 300 × 10^{−6} m |

D_{p} | 100 × 10^{−6} m | S_{lo} | 200 × 10^{−6} m |

Cp_{Si} | 707 J/kg/K | S_{tr} | 200 × 10^{−6} m |

k_{Si} | 149 W/mK | ρ_{Si} | 2330 kg/m^{3} |

Cp_{SiO2} | 730 J/kg/K | k_{SiO2} | 1.4 W/mK |

ρ_{SiO2} | 2200 kg/m^{3} | q_{pump/layer} | 1 W |

T_{amb} | 293 K | T_{in} | 298 K |
---|---|---|---|

h_{top} | 10 W/m^{2 }K | T_{ini} | 298 K |

L | 12000 × 10^{−6} m | h_{bot} | 550 W/m^{2}K |

d_{Si} | 100 × 10^{−6} m | W | 16000 × 10^{−6} m |

d_{SiO2} | 10 × 10^{−6} m | H_{pins} | 300 × 10^{−6} m |

D_{p} | 100 × 10^{−6} m | S_{lo} | 200 × 10^{−6} m |

Cp_{Si} | 707 J/kg/K | S_{tr} | 200 × 10^{−6} m |

k_{Si} | 149 W/mK | ρ_{Si} | 2330 kg/m^{3} |

Cp_{SiO2} | 730 J/kg/K | k_{SiO2} | 1.4 W/mK |

ρ_{SiO2} | 2200 kg/m^{3} | q_{pump/layer} | 1 W |

In this paper, Penryn [34], a 45 nm Intel Core 2 Duo processor architecture with doubled power to achieve 95 W for each core and 8 W for each cache is employed. The total power of each layer with two cores and two caches is 206 W. Each layer contains two identical cores and two identical caches but are located in reverse arrangements. Two cores in the first and third layers are located on one side. However, cores in the second and fourth layers are located on the opposite side. The power maps of the four layers are shown in Fig. 6.

## 4 Performance Analysis of Three-Dimensional-Stacked Integrated Circuits Compact Model

Figures 7(a) and 7(b) show the temperature field and flow directions for the layers. Because of the different flow directions in Figs. 7(a) and 7(b), the temperature fields of the layers with the same power distribution are also different. Table 4 shows the maximum temperature difference for all four layers in the two cases. In the opposite flow direction, the maximum temperature decreases 5 °C for the second layer, and 3.6 °C for the fourth layer. By contrast, for the same flow direction case, the maximum temperatures of the other two layers are relatively constant. These results suggest that, the flow direction between layers affects the cooling performance for nonuniform power maps. When the hotspots of the nonuniform power map are located upstream, the temperature of layers decreases. Therefore, this function can be used to achieve the best cooling effect for each layer by changing the flow direction based on the power map distribution.

Maximum temperature and leaking power | |||||
---|---|---|---|---|---|

First layer | Second layer | Third layer | Fourth layer | ||

Case 1 (Fig. 7(a)): uniform flow directions | Thermal model only | 62.3 (°C) | 60.0 (°C) | 52.6 (°C) | 55.3 (°C) |

With thermal-electrical codesign | 61.8 (°C) | 59.1 (°C) | 52.1 (°C) | 54.4 (°C) | |

Electrical leakage | 31.2 (W) | 30.9 (W) | 30.9 (W) | 30.8 (W) | |

Case 2 (Fig. 7(b)): various flow directions | Thermal model only | 62.6 (°C) | 55.1 (°C) | 54.8 (°C) | 51.6 (°C) |

With thermal-electrical codesign | 62.2 (°C) | 54.6 (°C) | 54.3 (°C) | 51.2 (°C) | |

Electrical leakage | 31.3 (W) | 31.0 (W) | 31.0 (W) | 30.9 (W) |

Maximum temperature and leaking power | |||||
---|---|---|---|---|---|

First layer | Second layer | Third layer | Fourth layer | ||

Case 1 (Fig. 7(a)): uniform flow directions | Thermal model only | 62.3 (°C) | 60.0 (°C) | 52.6 (°C) | 55.3 (°C) |

With thermal-electrical codesign | 61.8 (°C) | 59.1 (°C) | 52.1 (°C) | 54.4 (°C) | |

Electrical leakage | 31.2 (W) | 30.9 (W) | 30.9 (W) | 30.8 (W) | |

Case 2 (Fig. 7(b)): various flow directions | Thermal model only | 62.6 (°C) | 55.1 (°C) | 54.8 (°C) | 51.6 (°C) |

With thermal-electrical codesign | 62.2 (°C) | 54.6 (°C) | 54.3 (°C) | 51.2 (°C) | |

Electrical leakage | 31.3 (W) | 31.0 (W) | 31.0 (W) | 30.9 (W) |

Table 4 illustrates the temperature differences with and without thermal-electrical codesign in steady-state. Codesign based chip temperatures are normally 0.5–1 °C lower than the simulation results from thermal only compact thermal model, because of updated power map by implementing maximum temperatures. Table 4 also includes the leakage power of each layer based on the realistic average chip temperature. Although case 2 with various flow directions has lower maximum chip temperatures, it has higher leakage power considering the entire chip.

### 4.1 Transient Feature of the Compact Model.

One of the features of the compact model developed in this study is that it can simulate the transient temperature field. An important parameter setting for evaluating the thermal performance of the transient model is the power map, or the spatial heat flux distribution. The parameters required in the compact model appear in Table 3. Figure 8 shows the transient temperature field of all four layers at time 20, 40, and 80 s. At *t* = 20 s, the temperatures of a majority of the area of the chip are below 40 °C, except the two cores, which have reached 50 °C. After 20 s, because of heat conduction within each layer, the maximum temperatures of each layer increase, and the areas of temperatures higher than 50 °C also increase. Besides the heat conduction within each layer, heat conduction through pin fins between the layers in the vertical direction also has significant effects. As indicated at 80 s, nonuniform temperatures of the two caches with uniform power maps occur because of the vertical heat conduction through the cylindrical pin fins. As summarized in Table 5, all the realistic maximum chip temperatures with electrical codesign in the 3D-SICs at each time-step are normally 0.5–1.5 °C lower than the simulation results from only compact thermal model, because of updated power map by implementing maximum temperature of each layer into electrical codesign function. The leakage power of each layer based on the realistic average chip temperature also have been presented for each time-step in Table 5.

Maximum temperature and leakage power | |||||
---|---|---|---|---|---|

First layer | Second layer | Third layer | Fourth layer | ||

t = 20 s | Thermal model only | 55.5 (°C) | 51.4 (°C) | 49.7 (°C) | 49.8 (°C) |

With thermal-electrical codesign | 53.7 (°C) | 49.7 (°C) | 48.3 (°C) | 48.1 (°C) | |

Electrical leakage | 33.7 (W) | 30.6 (W) | 30.6 (W) | 30.6 (W) | |

t = 40 s | Thermal model only | 60.9 (°C) | 56.7 (°C) | 52.2 (°C) | 53.7 (°C) |

With thermal-electrical codesign | 59.8 (°C) | 55.4 (°C) | 51.5 (°C) | 52.5 (°C) | |

Electrical leakage | 31.0 (W) | 30.8 (W) | 30.8 (W) | 30.7 (W) | |

t = 80 s | Thermal model only | 62.2 (°C) | 59.6 (°C) | 52.6 (°C) | 55.1 (°C) |

With thermal-electrical codesign | 61.7 (°C) | 58.6 (°C) | 52.1 (°C) | 54.2 (°C) | |

Electrical leakage | 31.1 (W) | 30.9 (W) | 30.8 (W) | 30.7 (W) |

Maximum temperature and leakage power | |||||
---|---|---|---|---|---|

First layer | Second layer | Third layer | Fourth layer | ||

t = 20 s | Thermal model only | 55.5 (°C) | 51.4 (°C) | 49.7 (°C) | 49.8 (°C) |

With thermal-electrical codesign | 53.7 (°C) | 49.7 (°C) | 48.3 (°C) | 48.1 (°C) | |

Electrical leakage | 33.7 (W) | 30.6 (W) | 30.6 (W) | 30.6 (W) | |

t = 40 s | Thermal model only | 60.9 (°C) | 56.7 (°C) | 52.2 (°C) | 53.7 (°C) |

With thermal-electrical codesign | 59.8 (°C) | 55.4 (°C) | 51.5 (°C) | 52.5 (°C) | |

Electrical leakage | 31.0 (W) | 30.8 (W) | 30.8 (W) | 30.7 (W) | |

t = 80 s | Thermal model only | 62.2 (°C) | 59.6 (°C) | 52.6 (°C) | 55.1 (°C) |

With thermal-electrical codesign | 61.7 (°C) | 58.6 (°C) | 52.1 (°C) | 54.2 (°C) | |

Electrical leakage | 31.1 (W) | 30.9 (W) | 30.8 (W) | 30.7 (W) |

### 4.2 Various Heating Periods Feature of Compact Model.

The power map distribution of each layer can differ and change with time within the modeling framework. In this paper, the four-layer power map, including two cores and two caches for each layer, is changed to simulate real-time workloads. In addition to the power map parameters, the pumping powers for fluidic cooling for all layers remain constant at 1 W for each layer.

Figure 9 shows the temperature fields of the first layer for a simulation with two heating periods: period 1 from 0 to 60 s, and period 2 from 60 to 120 s. Initially, both cores in this layer are working at full capacity. Then the left core turns off at the beginning of the first period. After 60 s, at the beginning of the second period, the active cores are switched. As only one core and one cache are active in both periods, the total power for each layer is always 103 W. As the active core changes from one period to the next, Fig. 9 shows the temperature transition at the beginning of the second period. The resolution of the transition can be determined by user requirements.

## 5 Conclusion

This paper introduced a high-speed compact model used in the codesign of thermal and electrical aspects related to leakage. The model accounts for both spatially and temporally varying heat-flux distributions with high flexibility in each layer and the fluid pumping power with varying flow directions between layers. The temperature fields of 3D-stacked multilayer chips with pin-fin-enhanced microgap cooling have been simulated in a broad situation for both steady and transient states. The model can rapidly simulate the steady-state within 2 min, and transient state within 15 min, for steady and time-varying power maps. With thermal-electrical codesign included, the compact model can simulate the actual leakage power and provide accurate, updated temperatures. The model can be used for complex electronic operating process simulations, as an electrical-thermal codesign tool.

## Acknowledgment

This work was funded by the Defense Advanced Research Projects Agency (DARPA) ICECool Fundamentals (FUN) and Applications (APPS) programs, under the contracts HR0011-13-2-0008 and HR0011-14-1-0002 with Dr. Avram Bar-Cohen as Program Director.

## Funding Data

Defense Advanced Research Projects Agency (Grant Nos. HR0011-13-2-0008 and HR0011-14-1-0002; Funder ID: 10.13039/100000185).

## Nomenclature

*C*=coefficient in correlations of friction factor and Colburn

*j*factor*Cp*=specific heat at constant pressure, J/kg/K

*D*=diameter of fins, m

*h*=heat transfer coefficient, W/m

^{2}K*H*=pin-fin height, m

*k*=thermal conductivity, W/m K

*L*=chip length, m

*m*=mass-flow rate, kg/s

*q*=heat flux, W/m

^{2}*Q*=flow rate, L/hr

- Re =
Reynolds number

*S*=fin pitch, m

*t*=time, s

*T*=temperature, K

*W*=chip width, m

*α*=coefficient in correlations of friction factor and Colburn

*j*factor*δ*=thickness, m

*ρ*=density, kg/m

^{3}*Δ*=difference

## References

**143**(3), p.

*Handbook of 3D Integration: Design, Test, and Thermal Management*, Wiley-VCH, Weinheim, Germany, pp.