Thermal simulations for integrated circuits

At macroscopic level, Fourier's law states a relation between the transmitted heat per unit time per unit area and the gradient of temperature:

${\displaystyle q=-\kappa \nabla T}$

Where ${\displaystyle \kappa }$ is the thermal conductivity, [W·m⁻¹ K⁻¹].

Electronic systems work based on current and voltage signals. Current is the flow of charged particles through the material and these particles (electrons or holes), interact with the lattice of the crystal losing its energy which is released in form of heat. Joule Heating is a predominant mechanism for heat generation in integrated circuits[1] and is an undesired effect in most of the cases. For an ohmic material, it has the form:

${\displaystyle Q=j^{2}\rho }$

Where ${\displaystyle j}$ is the current density in [A·m⁻²], ${\displaystyle \rho }$ is the specific electric resistivity in [${\displaystyle {\Omega }}$·m] and ${\displaystyle Q}$ is the generated heat per unit volume in [W·m⁻³].[1]

The governing equation of the physics of the heat transfer problem relates the flux of heat in space, its variation in time and the generation of power by the following expression:

${\displaystyle \nabla \left(\kappa \left(T\right)\nabla T\right)+g=\rho C{\frac {\partial T}{\partial t}}}$

Where ${\displaystyle \kappa }$ is the thermal conductivity, ${\displaystyle \rho }$ is the density of the medium, ${\displaystyle C}$ is the specific heat, ${\displaystyle \alpha ={\frac {\kappa }{\rho C}}}$, the thermal diffusivity and ${\displaystyle g}$ is the rate of heat generation per unit volume. Heat diffuses from the source following the above equation and solution in an homogeneous medium follows a Gaussian distribution.

To get rid of the temperature dependence of ${\displaystyle \kappa }$, Kirchhoff transformation can be performed [2]

${\displaystyle \theta =T_{s}+{\frac {1}{\kappa _{s}}}\int _{T_{s}}^{T}\kappa (T)dT}$

where ${\displaystyle \kappa _{s}=\kappa \left(T_{s}\right)}$ and ${\displaystyle T_{s}}$ is the heat sink temperature. When applying this transformation, the heat equation becomes:

${\displaystyle \alpha \nabla ^{2}\theta +{\frac {\alpha }{\kappa _{s}}}g={\frac {\partial \theta }{\partial t}}}$

where ${\displaystyle \alpha ={\frac {\kappa }{\rho C}}}$ is called the diffusivity,[2] which also depends on the temperature. To completely linearize the equation, a second transformation is employed:

${\displaystyle \alpha _{s}\tau =\int _{0}^{t}\alpha (\theta )dt}$

yielding the expression:

${\displaystyle \nabla ^{2}\theta -{\frac {1}{\alpha _{s}}}{\frac {\partial \theta }{\partial \tau }}=-{\frac {g}{\kappa _{s}}}}$

Simple, direct application of this equation requires approximation. Additional terms arising in the transformed Laplacian are dropped, leaving the Laplacian in its conventional form.[2]

Although analytical solutions can only be found for specific and simple cases, they give a good insight to deal with more complex situations. Analytical solutions for regular subsystems can also be combined to provide detailed descriptions of complex structures. In Prof. Batty's work,[2] a Fourier series expansion to the temperature in the Laplace domain is introduced to find the solution to the linearized heat equation.

This procedure can be applied to a simple but nontrivial case: an homogeneous cube die made out of GaAs, L=300 um. The goal is to find the temperature distribution on the top surface. The top surface is discretized into smaller squares with index i=1...N. One of them is considered to be the source.

Taking the Laplace transform to the heat equation:

${\displaystyle \nabla ^{2}{\bar {\Theta }}-{\frac {s}{k_{s}}}{\bar {\Theta }}=0}$

where ${\displaystyle {\overline {\Theta }}=s\theta -\theta \left(\tau =0\right)}$

Function ${\displaystyle {\overline {\Theta }}}$ is expanded in terms of cosine functions for the ${\displaystyle x}$ and ${\displaystyle y}$ variables and in terms of hyperbolic cosines and sines for ${\displaystyle z}$ variable. Next, by applying adiabatic boundary conditions at the lateral walls and fix temperature at the bottom (heat sink temperature), thermal impedance matrix equation is derived:

${\displaystyle \Delta \theta _{i}=\sum _{j=1}^{N}R_{TH_{ij}}(t)P_{j}(t)}$

Where the index ${\displaystyle j}$ accounts for the power sources, while the index ${\displaystyle i}$ refers to each small area.

For more details about the derivation, please see Prof. Batty's paper,.[2] The below figure shows the steady state temperature distribution of this analytical method for a cubic die, with dimensions 300 um. A constant power source of 0.3W is applied over a central surface of dimension 0.1L x 0.1L. As expected, the distribution decays as it approaches to the boundaries, its maximum is located at the center and almost reaches 400K

Numerical solutions use a mesh of the structure to perform the simulation. The most popular methods are: Finite difference time-domain (FDTD) method, Finite element method (FEM) and method of moments (MoM).

The finite-difference time-domain (FDTD) method is a robust and popular technique that consists in solving differential equations numerically as well as certain boundary conditions defined by the problem. This is done by discretizing the space and time, and using finite differencing formulas, thus the partial differential equations that describe the physics of the problem can be solved numerically by computer programs.

The FEM is also a numerical scheme employed to solve engineering and mathematical problems described by differential equations as well as boundary conditions. It discretizes the space into smaller elements for which basis functions are assigned to their nodes or edges. Basis functions are linear or higher order polynomials. Applying the differential equation and the boundary conditions of the problem to the basis functions, a system of equations is formulated using either the Ritz or Galerkin method. Finally, a direct or iterative method is employed to solve the system of linear equations.[3] For the thermal case, FEM method is more suitable due to the nonlinearity nature of the thermal properties.

The previous example can be solved with a numerical method. For this case, the cube can by discretized into rectangular elements. Its basis functions can be chosen to be a first order approximation (linear):

${\displaystyle N_{i}^{e}={\frac {1}{2}}\xi \left(1\mp \zeta \right),\qquad i=1,4}$

${\displaystyle N_{i}^{e}={\frac {1}{2}}\eta \left(1\mp \zeta \right),\qquad i=2,5}$

${\displaystyle N_{i}^{e}={\frac {1}{2}}\left(1-\xi -\eta \right)\left(1\mp \zeta \right),\qquad i=3,6}$

where ${\displaystyle \zeta =2(z-z_{c})/h_{z}}$. If ${\displaystyle z_{c}=0}$, then ${\displaystyle \zeta =2z/h_{z}}$.

Using this basis functions and after applying Galerkin's method to the heat transfer equation, a matrix equation is obtained:

${\displaystyle \left[S\right]\left\{\theta \right\}+\left[R\right]{\frac {d}{dt}}\left\{\theta \right\}=\left\{B\right\}}$

where,

${\displaystyle R_{ij}=\int _{v}N_{j}N_{i}dV}$

${\displaystyle S_{ij}=k\int _{v}\nabla N_{j}.\nabla N_{i}dV}$

${\displaystyle B_{i}={\frac {k}{\kappa _{s}}}\int _{\Omega _{1}}N_{i}p(x,y)d\Omega +{\frac {k}{\kappa _{s}}}\int _{v}N_{i}gdV-kT_{o}\sum _{j=0}^{N_{D}}\int _{v}\nabla N_{j}^{D}.\nabla N_{i}dV}$.

This expressions can be evaluated by using a simple FEM code. For more details, please see.[3] The figure below shows the temperature distribution for the numerical solution case. This solution shows very good agreement with the analytical case, its peak also reaches 390 K at the center. The apparent lack of smoothness of the distribution comes from the first order approximation of the basis functions and this can be solved by using higher order basis functions. Also, better results might be obtained by employing a denser mesh of the structure; however, for very dense meshes the computation time increases a lot, making the simulation non-practical.

The next figure shows a comparison of the peak temperature as a function of time for both methods. The system reaches steady state in approximately ${\displaystyle 1ms}$.

The numerical methods such as FEM or FDM derive a matrix equation as shown in the previous section. To solve this equation faster, a method called Model order reduction can be employed to find an approximation of lower order. This method is based on the fact that a high-dimensional state vector belongs to a low-dimensional subspace .

Figure below shows the concept of the MOR approximation: finding matrix V, the dimension of the system can be reduced to solve a simplified system.

Therefore, the original system of equation:

${\displaystyle C\left\{x\right\}'+K\left\{x\right\}=F\left\{u\right\}}$

becomes:

${\displaystyle V^{T}CV\left\{z\right\}'+V^{T}KV\left\{z\right\}=V^{T}F\left\{u\right\}}$

Whose order is much lower than the original making the computation much less expensive. Once the solution is obtained, the original vector is found by taking the product with V.