This problem considers the performance of a heat sink designed for the thermal management of high-density electronic components. The heat sink is comprised of a base/spreader which in turn supports a number of plate fins exposed to flowing air. We model the flowing air through a simple convection heat transfer coefficient. From the engineering point of view, this problem illustrates the application of conduction analysis to an important class of cooling problems: electronic components and systems.
Our interest is in the conduction temperature distribution at the base of the spreader. The target is to study how the heat transfer occures with different parameters on our heat sink. The heat generated by high-density electronic components is such that it's very expensive to cool large structures (data center). The cooling optimization is consequent in the run for decreasing operating costs.
A classical thermal CPU cooler looks like this
![]() |
We are here going to describe how it is theorically working and how it is impleted with Feel++.
We consider here a classical "radiator" which is a CPU heat sink. Those types of coolers are composed with a certain number of plate fins exposed to flowing air or exposed to a ventilator. Regarding the periodicity and geometry of our concern, we can make our study on a characteristic element of the problem : a half cell of the heat sink single thermal fin with its spreader at the basis. Let's take a look at the geometry of our problem :
![]() |
Our study is avaible in 2 or 3 dimensions, depending on the application's parameters. You'll see later how to work with it. Let's see on which meshes we are working on :
![]() | ![]() |
The implementation of thoses parameters is described in the section Application Parameters.
Here the material parameter can be described with furthers parameters. We have, with \(i=1\) for the fin and \(i=2\) for the base :
The term \(\rho_i C_i\) corresponds to the heat volumetric capacity. In that way, we make possible the construction of a heat sink with \(2\) different materials. Here is a list of the well-known ones, \(\rho\) and \(C\) are gave at \(298 K\) :
Material | Thermal conductivity \(\kappa$ in $W.m^{-1}.K^{-1}\) | Density \(\rho$ in $kg.m^{-3}\) | Heat Capacity \(C\) in \(J.kg^{-1}.K^{-1}\) |
---|---|---|---|
Aluminium | 180 (alloys) or 290 (pure) | 2700 | 897 |
Copper | 386 | 8940 | 385 |
Gold | 314 | 19320 | 129 |
Silver | 406 | 10500 | 233 |
depth
in the application. L
in the application's parameters, its dimension is the meter. width
in the application's implementation.The following table displays the various fixed and variables parameters of this application.
Name | Description | Nominal Value | Range | Units |
---|---|---|---|---|
BDF parameters | ||||
\(time-initial\) | begining | \(0\) | ||
\(time-final\) | end | \(50\) | \(]0, 1500]\) | |
\(time-step\) | time step | \(0.1\) | \(]0,1[\) | |
\(steady\) | steady state | \(0\) | \(\{0,1\}\) | |
\(order\) | order | \(2\) | \([0, 4]\) | |
Physical parameters | ||||
\(L\) | fin's length | \(2\cdot 10^{-2}\) | \([0.02, 0.05]\) | \(m\) |
\(width\) | fin's width | \(5\cdot 10^{-4}\) | \([10^{-5}, 10^{-4}]\) | \(m\) |
\(deep\) | heat sink depth | \(0\) | \([0, 7\cdot 10^{-2}]\) | \(m\) |
Mesh parameter | ||||
\(hsize\) | mesh's size | \(10^{-4}\) | \([10^{-5},10^{-3} ]\) | |
Fin Parameters | ||||
\(\kappa_f\) | thermal conductivity | \(386\) | \([100,500]\) | \(W \cdot m^{-1} \cdot K^{-1}\) |
\(\rho_f\) | material density | \(8940\) | \([10^3,12\cdot 10^3 ]\) | \(kg\cdot m^{-3}\) |
\(C_f\) | heat capacity | \(385\) | \([10^2,10^3]\) | \(J\cdot kg^{-1} \cdot K^{-1}\) |
Base/spreader Parameters | ||||
\(\kappa_s\) | thermal conductivity | \(386\) | \([100,500]\) | \(W \cdot m^{-1} \cdot K^{-1}\) |
\(\rho_s\) | material density | \(8940\) | \([10^3,12\cdot 10^3 ]\) | \(kg\cdot m^{-3}\) |
\(C_s\) | heat capacity | \(385\) | \([10^2,10^3]\) | \(J\cdot kg^{-1} \cdot K^{-1}\) |
Heat Parameters | ||||
\(T_{amb}\) | ambient temperature | \(300\) | \([300,310]\) | \(K\) |
\(heat\_flux\) | heat flux \(Q\) | \(10^6\) | \([0 ,10^{6}]\) | \(W \cdot m^{-3}\) |
\(therm\_coeff\) | thermal coefficient \(h\) | \(10^3\) | \([0,10^3]\) | \(W\cdot m^{-2} \cdot K^{-1}\) |
The global domain is \(\varOmega = \varOmega_1 \cup \varOmega_2 \) where \(\varOmega_1\) is the fin's domain and \(\varOmega_2\) the spreader's domain. We note \(\partial\varOmega\) the border of the domain \(\varOmega\). The physical lines we are using will be noted as \(\Gamma_i\) such as described above. The following figure describes the parameters and the geometry we are using in the equations to solve our 3D issue : The following figures describe the parameters and the geometry we are using in the equations to solve our 2D or 3D issue :
![]() | ![]() |
Our concern satisfies the heat equation which reads
\begin{eqnarray} \sum_{i=1}^{2} \kappa_i \Delta T - \rho_i C_i \frac{ \partial T}{\partial t} & =& 0& \\ \kappa_1 \frac{\partial T}{\partial n} &= & 0 & \text{on } \Gamma_2 \text{ and } \Gamma_6\quad\quad \\ \kappa_2 \frac{\partial T}{\partial n} &= & 0 & \text{on } \Gamma_5 , \Gamma_7 \text{ and } \Gamma_8\quad\quad \\ \kappa_1 \frac{\partial T}{\partial n} &= &- h( T - T_{amb}) & \text{on } \Gamma_1\quad\quad \\ \kappa_2 \frac{\partial T}{\partial n} &= & Q(1-e^{-t}) & \text{on } \Gamma_4\quad\quad \\ T_{|\varOmega_1} &= &T_{| \varOmega_2} & \text{on } \Gamma_3\quad\quad \\ \kappa_1 \nabla T \cdot n &=& \kappa_2 \nabla T \cdot n & \text{on } \Gamma_3\quad\quad \end{eqnarray}
with \(i=1\) for the fin and \(i=2\) for the base and where \(\kappa_i\) is the thermal conductivity, \(\rho_i\) is the material's density ( \(kg.m^{-3}\) in the SI unit), \(C_i\) the heat capacity and \(T\) the temperature at a precise point (in 2D or 3D). To see how it has been coded, you can read Implemtation .
The problem requieres that the temperature and heat flux are continute on \(\Gamma_3\). Considering the problem's geometry, we also impose zero heat flux on the vertical surfaces of the spreader. Let's detail the conditions we have imposed :
Theses conditions have been coded as explained in the section Equations.
Let's apply the method to our concern, we introduce the test function \(v\) and we integrate the main equation, which reads now as :
\begin{equation*} \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} - \kappa_i \int_{\varOmega_i} v\Delta T & = 0 \end{equation*}
We integrate by parts, which leads to :
\begin{equation*} \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} - \kappa_i \int_{\partial \varOmega_i} {(\nabla T \cdot n) v} & = 0 \end{equation*}
then, by decomposing the borders \(\partial\varOmega_i\), we obtain :
\begin{multline*} \displaystyle{- \kappa_1 \int_{\Gamma_1}{(\nabla T \cdot n) v} - \kappa_2 \int_{\Gamma_4} {(\nabla T \cdot n) v} - \kappa_1 \int_{\Gamma_{2,6}}{(\nabla T \cdot n) v} - \kappa_2 \int_{\Gamma_{5,7,8}}{(\nabla T \cdot n) v} } \quad + \\ \displaystyle{ \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} - \kappa_i \int_{\partial\varOmega_i \cap \Gamma_3}{(\nabla T \cdot n) v} } = 0 \end{multline*}
Now, we apply the conditions \((2)\), \((3)\), \((4)\) and \((5)\) which brings us to :
\begin{equation*} \int_{\Gamma_1}{hv(T-T_{amb})} - \int_{\Gamma_4} {vQ(1-e^{-t})} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} - \underbrace{\kappa_i \int_{\partial\varOmega_i \cap \Gamma_3}{(\nabla T \cdot n) v}}_{\text{=0}} & = 0 \end{equation*}
Now we apply the boundary conditions ( \((7)\) ) which results in :
\begin{equation*} \begin{split} \displaystyle{ h \int_{\Gamma_1}{v(T-T_{amb})} - \int_{\Gamma_4} {vQ(1-e^{-t})} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} } & = 0 \end{split} \end{equation*}
We can now start to transform the equation by puting in the right hand the known terms :
\begin{equation*} \begin{split} \displaystyle{ h \int_{\Gamma_1}{v T} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{ \partial T}{\partial t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} } & = \int_{\Gamma_4} {vQ(1-e^{-t})} + hT_{amb}\int_{\Gamma_1}{v} \end{split} \end{equation*}
We discretize \(\displaystyle{\frac{\partial T}{\partial t}}\) where \(\delta t\) is the time step, such as:
\begin{equation*} \begin{split} \displaystyle{ h \int_{\Gamma_1}{v T} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{T^{n+1} - T^n}{\delta t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T} } & = \int_{\Gamma_4} {vQ(1-e^{-t})} + hT_{amb}\int_{\Gamma_1}{v} \end{split} \end{equation*}
Finally we obtain :
\begin{equation*} \color{red} \begin{split} \displaystyle{ h \int_{\Gamma_1}{v T} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v\frac{T^{n+1}}{\delta t} + \kappa_i \int_{\varOmega_i} {\nabla v \cdot \nabla T}} & = \displaystyle{ \int_{\Gamma_4} {vQ(1-e^{-t})} + hT_{amb}\int_{\Gamma_1}{v} + \sum_{i=1}^{2} \rho_i C_i \int_{\varOmega_i} v \frac{T^n}{\delta t} } \end{split} \end{equation*}
This is that equation which is implemented in the application feel_heatsink
.
The parameters of the application are implemented such as
To be able to calculate the surfaces in further dimension without changing the code, we have given the same names for the faces we were interested in. In 2D \(\Gamma_i\) represents a line whereas in 3D it represents a surface. The calculation of those surfaces which makes possible the calculation of averages temperature is as follow :
First we start by calculate the non-steady state which means that we integrate all the time-independant terms, which is done with :
Then, to compute the transient state, which means time dependant terms, you have to initialize the temperature (which is initialized as \(T_{amb}\) on \(X_h\) space) and create a new vector \(F_t\) which corresponds to the time dependent term. The code is as follow :
As you can see in the equation's implementation above, there are two ouputs :
Gmsh
format : this file contains the entire mesh and the temperatures associated to each degrees of freedom of the mesh. To open it, you juste have to do as you always do with Gmsh
: gmsh
heatsink-1_0.msh. You will obtain the figure with the different temperatures, you are now able to click on "play" with its significative logo and admire the evolution averages
file : this file is completed at each time step, each line contains the current time, the average temperature on \(\Gamma_4\) (surface where is the contact between the heat sink and the heat source) and the average temperature on \(\Gamma_1\). To analyze this file, we recommend you to work with Octave
which is an open-source software similar to Matlab
.~/feel/heatsink/Simplex_'.'.'/0.000'/and try :
> octave octave:1> M=load('averages'); octave:2> plot(M(:,1),M(:,2)) octave:3> plot(M(:,1),M(:,3)) octave:4> plot(M(1:70,1),M(1:70,2)) octave:5> plot(M(1:70,1),M(1:70,3))The \( 4^{th}\) and \( 5^{th}\) lines are here to observe the transient state.
To make easier the use of this application, we recommand you to use the configurations files. This is the fastest way : to do it, you juste have to create the file heatsink.cfg
and place it in the same directory that your application's executable.
We have created \(3\) typical cfg
files such as :
This file is the only modification you will have to bring to the application, in that way you won't have to compile each time the files (except for heatsink.cpp
if you want to increase the order and/or the dimension, in that case you'ill have to modify this parameter at then end of the file in the main
method).
Here are some results of the 2D simulations considering different configurations files. The figures have been extracted thanks to Gmsh
and Octave
:
![]() |
heatsink_1.cfg : steady state, spreader and fin in copper, \(Q=1e6\) and \(h=1e3\) |
![]() | ![]() |
heatsink_1.cfg : transient state on \(\Gamma_4\) | heatsink_1.cfg : transient state on \(\Gamma_1\) |
Here is the result of 3D simulations considering the following configurations :
![]() |
heatsink_3.cfg : spreader and fin in copper, \(Q=1e6\) and \(h=1e3\) |
![]() | ![]() |
heatsink_3.cfg : transient state on \(\Gamma_4\) | heatsink_3.cfg : transient state on \(\Gamma_1\) |