Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
Heat Sink
Author
Christophe Prud'homme
Baptiste Morin
Date
2011-06-28




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

complete_cooler.png
Mesh of a classical CPU cooler

We are here going to describe how it is theorically working and how it is impleted with Feel++.

Description

Domain

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 :

sink_geom.png
Geometry of heat sink

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 :

mesh_2d.png
mesh_3d.png
2D mesh
3D mesh

Inputs

The implementation of thoses parameters is described in the section Application Parameters.

Material

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 thermal conductivity \(\kappa_i\)
  • the material's density \(\rho_i\)
  • the heat capacity of the material \(C_i\)

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\) :

MaterialThermal 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

Physical

  • Depth
    This parameter is only to take into account for the 3D simulation. It represents the depth of the caracteristical heat sink and is called depth in the application.
  • Length
    You can also parameterize the length of the fin. This one is called L in the application's parameters, its dimension is the meter.
  • Width
    Typically, this parameter is linked with constructor's standards. This parameter is called width in the application's implementation.

Thermal

  • Heat flux
    It represents the heat flux brought by the electronic component at the bottom of the base. Here it's typically the heat brought by the processor.
  • Thermal coefficient
    The thermal coefficient \(h\) named \(therm_coeff\) in the application is representative of the heat transfer between the fin and the air flow.
  • Ambien temperature}
    This parameter called \(T_{amb}\) represents the temperature around the heat sink at the beginning. That means the ambient temperature before the computer is turned on.

Summary table

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}\)

Theory

Figure

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 :

figure_2d.png
figure_3d.png
2D geometry details
3D geometry details

Equations

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 .

Boundary Conditions

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 :

  • Homogeneous Neumann condition ( \((2)\) ) and ( \((3)\) ) : it represents the fact that the heat flux is only vertical (for \(\Gamma_6\) and \(\Gamma_7\)) or the fact that the heat flux is only provided by \(\Gamma_4\) (for \(\Gamma_2\) and \(\Gamma_5\)).
  • Homogeneous Neumann condition ( \((4)\) ) : it imposes that the heat flux is brought by this surface (it mathematically represents that the heat sink is placed on the heat source).
  • Non-homogeneous Neumann condition ( \((5)\) ) : this boundary condition represents the transient state for the heat transfer calculation.
  • Temperature continuity ( \((6)\) ) : it imposes that the temperature is continute at the interface between the two materials (if there are two materials, we can also have the same one for the two pieces).
  • Heat flux continuity ( \((7)\) ) : it represents that the heat flux is continute at the interface between the two materials. Literally, it means that the two flows offset each other.

Theses conditions have been coded as explained in the section Equations.

Finite Element Method

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.

Implementation

Application Parameters

The parameters of the application are implemented such as

inline
Feel::po::options_description
makeOptions()
{
Feel::po::options_description heatsinkoptions( "heatsink options" );
heatsinkoptions.add_options()
// mesh parameters
( "L", Feel::po::value<double>()->default_value( 20 ),
"dimensional length of the sink (in meters)" )
( "width", Feel::po::value<double>()->default_value( 0.5 ),
"dimensional width of the fin (in meters)" )
// 3D parameter
( "depth", Feel::po::value<double>()->default_value( 0 ),
"depth of the mesh (in meters) only in 3D simulation" )
// thermal conductivities parameters
( "kappa_s", Feel::po::value<double>()->default_value( 386e-3 ),
"thermal conductivity of the base spreader in SI unit W.mm^{-1}.K^{-1}" )
( "kappa_f", Feel::po::value<double>()->default_value( 386e-3 ),
"thermal conductivity of the fin in SI unit W.mm^{-1}.K^{-1}" )
// density parameter
( "rho_s", Feel::po::value<double>()->default_value( 8940e-9 ),
"density of the spreader's material in SI unit kg.mm^{-3}" )
( "rho_f", Feel::po::value<double>()->default_value( 8940e-9 ),
"density of the fin's material in SI unit kg.mm^{-3}" )
// heat capacities parameter
( "c_s", Feel::po::value<double>()->default_value( 385 ),
"heat capacity of the spreader's material in SI unit J.kg^{-1}.K^{-1}" )
( "c_f", Feel::po::value<double>()->default_value( 385 ),
"heat capacity of the fin's material in SI unit J.kg^{-1}.K^{-1}" )
// physical coeff
( "therm_coeff", Feel::po::value<double>()->default_value( 1e3 ),
"thermal coefficient" )
( "Tamb", Feel::po::value<double>()->default_value( 300 ),
"ambiant temperature" )
( "heat_flux", Feel::po::value<double>()->default_value( 1e6 ),
"heat flux generated by CPU" )
( "steady", Feel::po::value<bool>()->default_value( false ),
"if true : steady else unsteady" )
// export
( "export-matlab", "export matrix and vectors in matlab" );
return heatsinkoptions.add( Feel::feel_options() );
}

Surfaces

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 :

/*
* Calculate the two surfaces used for averages calculation
*/
surface_base = measure( _range= markedfaces( mesh, "gamma4" ) );
surface_fin = measure( _range= markedfaces( mesh, "gamma1" ) );

Equations

First we start by calculate the non-steady state which means that we integrate all the time-independant terms, which is done with :

/*
* Right hand side construction (steady state)
*/
auto l = form1( _test=Xh );
auto lt = form1( _test=Xh );
l = integrate( _range= markedfaces( mesh, "gamma1" ), _expr= therm_coeff*Tamb*id( v ) );
/*
* Left hand side construction (steady state)
*/
auto a = form2( Xh, Xh );
a = integrate( _range= markedelements( mesh,"spreader" ), _expr= kappa_s*gradt( T )*trans( grad( v ) ) );
a += integrate( _range= markedelements( mesh,"fin" ), _expr= kappa_f*gradt( T )*trans( grad( v ) ) );
a += integrate( _range= markedfaces( mesh, "gamma1" ), _expr= therm_coeff*idt( T )*id( v ) );
M_bdf->start();
//from now if the option "steady" is set to true then M_bdf->setSteady will set time-step=time-final=1e30
if ( steady )
{
M_bdf->setSteady();
}
a +=
integrate( _range=markedelements( mesh, "spreader" ), _expr=rho_s*c_s*idt( T )*id( v )*M_bdf->polyDerivCoefficient( 0 ) )
+ integrate( _range=markedelements( mesh, "fin" ), _expr=rho_f*c_f*idt( T )*id( v )*M_bdf->polyDerivCoefficient( 0 ) );

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 :

for ( M_bdf->start(); M_bdf->isFinished()==false; M_bdf->next(T) )
{
lt.zero();
// update right hand side with time dependent terms
auto bdf_poly = M_bdf->polyDeriv();
lt =
integrate( _range=markedelements( mesh, "spreader" ), _expr=rho_s*c_s*idv( bdf_poly )*id( v ) ) +
integrate( _range=markedelements( mesh, "fin" ), _expr=rho_f*c_f*idv( bdf_poly )*id( v ) );
lt +=
integrate( _range= markedfaces( mesh,"gamma4" ), _expr= heat_flux*( 1-math::exp( -M_bdf->time() ) )*id( v ) );
lt += l;
a.solve( _solution=T, _rhs=lt );
Tavg = integrate( _range=markedfaces( mesh,"gamma4" ), _expr=( 1/surface_base )*idv( T ) ).evaluate()( 0,0 );
Tgamma1 = integrate( _range=markedfaces( mesh,"gamma1" ), _expr=( 1/surface_fin )*idv( T ) ).evaluate()( 0,0 );
// export results
out << M_bdf->time() << " " << Tavg << " " << Tgamma1 << "\n";
M_exporter->step( M_bdf->time() )->addScalar( "Tavg", Tavg );
M_exporter->step( M_bdf->time() )->addScalar( "Tgamma1", Tgamma1 );
M_exporter->step( M_bdf->time() )->add( "Temperature", T );
M_exporter->save();
}

Results

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.
    If it is installed, open a command line and go to
    ~/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.

How to use it ?

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 :

# file heatsink_1.cfg
# spreader and fin in copper
# 2D simulation
hsize=1e-4
kappa_s=386 # W/m/K
c_s=385
rho_s=8940
kappa_f=386 # W/m/K
c_f=385 #J/kg/K
rho_f=8940
L=15e-3
width=5e-4
therm_coeff=1000 #W/(m2K)
heat_flux=1e6
[bdf]
order=2
time-step=0.05
time-final=100
steady=0
[exporter]
format=gmsh
# file heatsink_3.cfg
# spreader in copper
# fin in aluminium
# 3D simulation
hsize=3e-4
kappa_s=386 # W/m/K
c_s=385
rho_s=8940
kappa_f=386 # W/m/K
c_f=385 #J/kg/K
ho_f=8940
L=15e-3
width=5e-4
deep=4e-2
therm_coeff=1000 #W/(m2K)
heat_flux=1e6
[bdf]
order=2
time-step=0.05
time-final=100
steady=0
[exporter]
format=gmsh

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).

Figures

2D cases

Here are some results of the 2D simulations considering different configurations files. The figures have been extracted thanks to Gmsh and Octave :

heatsink_1.png
heatsink_1.cfg : steady state, spreader and fin in copper, \(Q=1e6\) and \(h=1e3\)
heatsink_1_gamma4.png
heatsink_1_gamma1.png
heatsink_1.cfg : transient state on \(\Gamma_4\)
heatsink_1.cfg : transient state on \(\Gamma_1\)

3D cases

Here is the result of 3D simulations considering the following configurations :

heatsink_3.png
heatsink_3.cfg : spreader and fin in copper, \(Q=1e6\) and \(h=1e3\)
heatsink_3_gamma4.png
heatsink_3_gamma1.png
heatsink_3.cfg : transient state on \(\Gamma_4\)
heatsink_3.cfg : transient state on \(\Gamma_1\)