The rest of the dissertation is organized as follows. In chapter 2, the problem formulation, the simulation algorithm and the experimental results of the proposed deterministic thermal simulation method are detailed. In chapter 3, the problem formulation, the models of leakage powers, the simulation algorithm and the experimental results of the proposed statistical thermal simulation method are detailed. In chapter 4, the problem formulation, the simulation algorithm and the experimental results of the proposed look-up table based thermal simulation method for early design stages of 3-D ICs are detailed. Finally, the conclusion of this dissertation is presented in chapter 5.
Chapter 2
Simulation Method I – Full-Chip Thermal Analysis for Early Design Stages via
Generalized Integral Transforms
As addressed by [40–42, 56–59], temperature-aware design should be brought to early design stages such as thermal-aware floor-planning and placement. Therefore, the proposed simulation method in this chapter is on efficiently proving the on-chip temperature profile for early design stages. This chapter is organized as follows. First, the thermal model for early design stages is presented in section 2.1. The generalized integral transform (GIT) based computational formula of the on-chip temperature profile and the proposed evaluating algorithms are described in sec-tion 2.2. After that, in secsec-tion 2.3, a method with a hybrid scheme of the GIT based analytical and the FDM based numerical formulations will be addressed for simulating the temperature profile of the stacked dies and package structure. Finally, the experimental results are given in sections 2.5.
2.1 Thermal Modeling for Early Design Stages
To estimate temperature with a reasonable accuracy and a little computational effort, instead of employing a detail thermal model for the post-layout thermal verification, a compact thermal model is essential for fast temperature simulation in early design stages. A popular compact thermal model of the chip for early design stages is illustrated in Figure 2.1 [56–59]. This model consists of three portions: the primary heat flow path, the secondary heat flow path, and the heat transfer characteristic of each macro/block on the silicon die. The primary heat flow path is composed of thermal interface material, heat spreader and heat sink. The secondary
heat flow path contains interconnect layers, I/O pads and the print circuit board (PCB). The functional blocks are modeled as many power generating sources attached to a thin layer close to the top surface of die with the thickness being equal to the junction depth1. The major concerns
Heat Sink
Thermal Interface Material
Ambient Air
Interconnect Layers
I/O Pads & PCB
h s
Secondary Heat Transfer Path Ambient Air
Primary Heat Transfer Path
Thin layer with
junction depth thickness Heat Spreader
z=-j
dFigure 2.1: Compact thermal model for early design stages.
of early-stage temperature-aware optimization procedure are to reduce the temperature or the thermal gradient of die. Thus, we focus on estimating the temperature profile of die.
According to energy conservation law, the changing rate of energy in a unit volume of sub-strate equals to the conduction heat through the unit volume [65]. Figure 2.2 illusub-strates this heat conduction mechanism. In Figure 2.2, dE/dt is the energy change rate for the unit volume and is equal to σ∆x∆y∆z∂T/∂t. The conduction heat flowing into the unit volume is equal to the sum of q|x0 = −κ∆y∆z ∂T/∂x|x0, q|y0 = −κ∆x∆z ∂T/∂y|y0 and q|z0 = −κ∆x∆y ∂T/∂z|z0. The conduction heat flowing outward the unit volume is the sum of q|x0+∆x = −κ∆y∆z ∂T/∂x|x0+∆x, q|y0+∆y = −κ∆x∆z ∂T/∂y|y0+∆yand q|z0+∆z = −κ∆x ∆y∂T/∂z |z0+∆z. p∆x∆y∆z is the energy gen-eration rate of that unit volume. κ and σ are the thermal conductivity, and the product of the
material density and the heat in the unit volume, respectively. p is the density of power
A Unit Volume of Die Energy Conservation and Heat Flow Mechanism
y
0q dE q y
0+∆ y
dt
Heat Conduction of the Unit Volume
(
0 0) (
0 0) (
0 0)
Figure 2.2: Energy conservation law and the heat conduction equation.
Based on the heat conduction mechanism, the temperature Td(r, t) of die can be governed by the following heat transfer equations [51–59]
σ(Td)∂Td(r, t) arbitrary function on bs, and ∂/∂nbs is the differentiation along the outward direction which is normalized to bs.
To provide reasonable accuracy of the temperature estimation with a little computational effort during executing early-stage temperature-aware optimization procedures, heat transfer
coefficients on the boundary surfaces of are suggested to be appropriately modeled [56–59].
Based on the model proposed by [56–59], the thermal model of the primary heat flow path can be modeled as an effective heat transfer coefficient hpby combining the effect of each compo-nent on the primary heat flow path. Since the detailed layout of interconnects is not available in early design stages, interconnect layer is modeled as an equivalent thermal resistance based on the densities and the regularity structure assumption of metal and dielectric material [56–59].
Furthermore, the I/O pads and print circuit board (PCB) can also be modeled as an effective thermal resistance by using the technique proposed by [68]. Then, the equivalent heat transfer coefficient hs of these successively connected thermal resistors can be calculated by the tech-nique shown in [49]. After hpand hshave been obtained, fbs(r)’s for the top and bottom surfaces are set to hsTaand −hpTa, respectively [5, 49, 51–59]. Here, Ta is the ambient air temperature.
Because of the chip and package structures, the area of vertical surface is strictly less than the area of horizontal surface, and the thermal conductivity of air is much less than the thermal conductivities of primary and secondary heat flow paths. Therefore, the boundary condition of each vertical surface can be reasonably set to be adiabatic [5].
Generally, the values of κ(Td) and σ(Td) are temperature dependent. The difference of peak temperature is about 5◦C between the result with temperature-dependent thermal parameters and the result with constant thermal parameters at 25◦C [55]. In current VLSI design, the on-die temperature can be in the degree of 100◦C. Under this situation, this difference may lead to about 5% error for the peak temperature of die. However, the effort to amend this error is relatively high because several iterations of the thermal simulation have to be executed for correcting the difference caused by the temperature dependences of κ(Td) and σ(Td). For practical purposes, these thermal parameters are usually treated as appropriate constants while performing temperature-aware floor-planning and placement [40–42].
The value of each thermal parameter can be found by applying a 1-D thermal circuit shown in Figure 2.3 to estimate the roughly average steady state temperature of the die. In Figure 2.3, the values of thermal resistors are Rs = 1/Adzhs, Rp = 1/Adzhp and Rdie = DT/κAdz. Tavg(z) is the average steady state temperature on the lateral planes at arbitrary z position of die. Here, Rdiecan be viewed as a variable resistor when obtaining Tavg(z) at certain z position. PT is the
( )
T
avgz Ambient Air
Ambient Air
R s
R die
R p
Die
Primary Heat Flow Path
A dz
A dz
D
TSubstrate
R s
R die
R p
P T
T
aT
aSecondary Heat Flow Path
Figure 2.3: The 1-D thermal model for estimating the roughly steady state average temperature of die. The modeled thermal resistance network is shown in the right hand side.
total average steady state power consumption of die. Adz is the cross area of die normal to the z-direction, and DT is the thickness of die. The computation flow is processed as follows. In the beginning, Rdieis calculated using the thermal conductivity of die at the room temperature.
After thermal resistances Rs, Rpand Rdieare obtained, the average rising temperature Tavgof die is equal to
Tavg = Tavg(0)+ Tavg(−Lz)
2 , (2.3)
where Tavg(0) and Tavg(−Lz) can be obtained by solving the temperatures in the 1-D thermal circuit shown in Figure 2.3.
Once Tavg is calculated, Rdie is re-calculated using the thermal conductivity of die at Tavg. This calculating procedure is repeated until Tavg converges. After that, the thermal parame-ters are calculated at the average temperature. With these estimated thermal parameparame-ters, the error of the peak temperature between the simulated temperature profiles for with and without considering the temperature dependence of the thermal parameters can be reduced.
With the above models, the heat diffusion equations for the rising temperature profile of die,
T(r, t)= Td(r, t) − Ta, in early design stages can be written as
Here, κ and σ are the thermal conductivity, and the product of the material density and the specific heat of die got by using the roughly steady state average temperature, respectively, and the initial condition T (r, 0)= 0.
As shown in Figure 2.1, after dividing the layer with power generating sources into MN grids with M and N are numbers of divisions in x- and y- directions, respectively, the power density profile p(r, t) stated in equation (2.4) can be written as
p(r, t)= waveform of grid cell (m, n) in the thin layer with thickness jd.
For the transient (dynamic) thermal simulation, pmn(t) is a time-interval function with the magnitude of each interval being equal to the average power density of each time interval. We should note that the thermal time-constant of heat conduction is much larger than the clock period of circuit [51, 56]. As indicated by [56], the temperature takes at least 100K cycles to rise 0.1◦C. Practically, the time interval specified by the user can be much larger than the clock period of circuit. For the steady state thermal simulation, the input power profile is usually set to the steady power profile (the average power profile for a very long time period estimation) [51, 56]. Therefore, pmn(t) can be reasonably viewed as a step function with the magnitude being equal to its average power density for a long time period.
With the above discussion and governing equations (2.4)–(2.7), our goal is to get the rising
Input
•The chip geometry and chip package configuration
•The power density of grid cells
Computational formulas construction
Fast temperature evaluating algorithm
Perform the proposed 2D-LTS-FFT to the power density map of grid cells
Perform the proposed 2D-STL-FFT to get the average temperatures of grid cells
Output
The temperature profile of grid cells Construct a set of
appropriate ortho-normal bases for approximating temperature distribution
Construct a system to calculate time varying coefficients by using Galerkin’s scheme
Get computational formulas for average temperatures of grid cells
Figure 2.4: The executing flow of the proposed GIT based thermal simulation method.