High Performance Heat-spreader for Heatsink

Using heat sink specification provided by manufacturers requires understanding of thermal design. Most design engineer not trained in heat transfer can erroneously rely on the published specifications for heat-sink selection. The heat flux present at the semiconductor transfer surface can reduce the published performance specification due to thermal resistance of the material at the mating surface. This reduction in heatsink performance can easily be mitigated by careful selection of  heat-spreader.

Matlab generated thermal image. Thermal system without optimization.
Matlab generated thermal image. Thermal system without optimization.

Matlab generated thermal image. With optimization
Matlab generated thermal image. Thermal System with optimization

Finite element or more accurately, the finite difference method provides a powerful and flexible method of optimizing the performance of the entire heatsink. Most CAD provide powerful simulation features with intuitive user interface to simulate the performance of thermal system. However, optimizing these thermal system in CAD may require tens or even hundreds of simulations. By modeling the system in Matlab, engineers can easily run these simulations by creating a loop and varying parameters of the design. The basics of finite difference modeling is introduces and the complete model is developed and simulated to show the performance increase by adding a heat spreader to the system.

Shiraz Macuff, Nesh Basnet, Subash Shrestha, Anil Rajpatei, Francisco Santos,

Adrian Serrano,Lhakpa Dhondup, Lobsang Lekshey

Example heat sink

 Assem1Cu  Assem1Cu - Sheet1

Figure 1 and figure 2 above shows an example of two configurations of the same heatsink. The TO-247 transistor in figure 1 and the  TO-220 in figure 2 both dissipate the same amount of power. However the absence of a high performance heat-spreader on the heatsink results in the TO-220 transistor in figure 2 running at a higher temperature due to the higher heat-flux at the mating surface.

Two versions of the heat-sinks are examined. One without a heat spreader and one with a high performance heat-spreader made of phase change system. The problem is simplified for two dimensional analysis  by looking at the heat-sink from the top face. Only thermal conduction is considered since we are only varying the performance of the heat-spreader which contributes to conduction performance. Mounting screws and mounting flange are omitted since they will be held constant while comparing the two models. Because of the limitation of analytic method,  finite difference method is used instead to build a foundation for other irregular geometries. References 1 and 2 are highly recommended for thermal studies.

Heat flux

The heat flux present at the mating surface of the transistor is given by \frac{q}{A} where  A is the mating area given by  A=H_1*W_1. From this we see that given the same power dissipation in both the TO-220 and TO-247 transistor, the heat flux will be much more on the TO-220 transistor. The thermal resistance of the heat-spreader material will produce a thermal gradient parallel to L1  and parallel to W1 at the base  of the transistor. This thermal gradient reduces the efficiency of the entire system.

2D view of the problem

GriddedHeatsinkReducing the problem to 2D requires that we examine the top face of the assembly as shown in figure 3. The heat spreader and the heat-sink are considered to be the same material in this case aluminum. The superimposed grid allow us to formalize the equations for analysis. The thermal gradient from nodes   (m,n) to  (m+1,n) is given by \frac{\partial T}{\partial x}\Big|_{m+\frac{1}{2},n}\approx \frac{T_{m+1,n}-T_{m,n}}{\Delta x}.

Similarly the gradient between nodes surrounding node (m,n) are given by:

 (m,n) to  (m-1,n)\frac{\partial T}{\partial x}\Big|_{m-\frac{1}{2},n}\approx \frac{T_{m,n}-T_{m-1,n}}{\Delta x}

 (m,n) to  (m,n+1)\frac{\partial T}{\partial x}\Big|_{m,n+\frac{1}{2}}\approx \frac{T_{m,n+1}-T_{m,n}}{\Delta x}

 (m,n) to  (m,n-1)\frac{\partial T}{\partial x}\Big|_{m,n-\frac{1}{2}}\approx \frac{T_{m,n}-T_{m,n-1}}{\Delta x}.

Taking the second derivative

\frac{\partial^2 T}{\partial x^2}\Big|_{m,n}\approx \frac{\frac{\partial T}{\partial x}\Big|_{m+\frac{1}{2},n}-\frac{\partial T}{\partial x}\Big|_{m-\frac{1}{2},n}}{(\Delta x)^2}=\frac{T_{m+1,n}+T_{m-1,n}-2T_{m,n}}{(\Delta x)^2} \frac{\partial^2 T}{\partial y^2}\Big|_{m,n}\approx \frac{\frac{\partial T}{\partial y}\Big|_{m,n+\frac{1}{2}}-\frac{\partial T}{\partial y}\Big|_{m,n-\frac{1}{2}}}{(\Delta y)^2}=\frac{T_{m,n+1}+T_{m,n-1}-2T_{m,n}}{(\Delta y)^2} Two dimensional heat conduction equation is given below with the finite difference approximation \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial x^2}=0\approx\frac{T_{m+1,n}+T_{m-1,n}-2T_{m,n}}{(\Delta x)^2}+\frac{T_{m,n+1}+T_{m,n-1}-2T_{m,n}}{(\Delta y)^2} Making \Delta y=\Delta x results in

1.0) T_{m+1,n}+T_{m-1,n}+T_{m,n+1}+T_{m,n-1}-4T_{m,n}=0

For the four nodes surrounding  T_{m,n}, the four remaining equations are:

1.1) T_{m,n}+T_{m-1,n}+T_{m,n+1}+T_{m,n-1}-4T_{m+1,n}=0

1.2) T_{m,n}+T_{m+1,n}+T_{m,n+1}+T_{m,n-1}-4T_{m-1,n}=0

1.3) T_{m,n}+T_{m+1,n}+T_{m-1,n}+T_{m,n-1}-4T_{m,n+1}=0

1.4) T_{m,n}+T_{m+1,n}+T_{m-1,n}+T_{m,n+1}-4T_{m,n-1}=0

Using the method shown above, a system of equation is derived and solved for the nodal temperatures. Note that this does not account for heat generation nor does account for convection at the edges.


% This is an example of two dimensional heat distribution
% on a square plate. For illustration,
% the number of nodes are limited 12.Reference [1] and [4] used.
% The figure below shows a square plate with the
% nodes at the boundaries held at specific temperatures
% temperatures. The temperatures at internal nodes 1
% trough 4 is solved for in the code below.
%
%
%                     T=500C
%         o---------o---------o--------o
%         |         |         |        |
%         |         |         |        |
%         |         |         |        |
%         o--------o1--------o2--------o
%         |         |         |        |
%  T=100C |         |         |        |  T=100C
%         |         |         |        |
%         o--------o3--------o4--------o
%         |         |         |        |
%         |         |         |        |
%         |         |         |        |
%         o---------o---------o--------o
%                     T=100C
%
% the equations for the four nodes are given below
% Node1) -4*T1+T2+T3+0*T4=-600
% Node2) T1-4*T2+0*T3+T4=-600
% Node3) T1+0*T2-4*T3+T4=-200
% Node4) 0*T1+T2+T3-4*T4=-200
% clear all previous variables
clear all
%clear the screen
clc
% load matrix with coefienct of the equations
M= [-4 1 1 0;
    1 -4 0 1;
    1 0 -4 1;
    0 1 1 -4];
% load vector B with the right hand side of the equations
B=[-600;-600;-200;-200];
% multiply the inverse of M with B to solve
%for variables T1 though T4
X=inv(M)*B;
Results:</pre>
250 250 150 150
<pre>

The example above does not include convection and heat generation

The 2D problem with convective heat transfer and heat generation at the boundaries.

To include convective heat transfer and heat generation at the boundaries we need to modify the nodal equation above. griddedHeatsink2CroppedFor the three nodes surrounding node (m,n)  , where the right side is exposed to convection, the energy balance for node m,n is given below. -K \Delta y \frac{T_{m,n}-T_{m-1,n}}{\Delta x}-K \frac{\Delta x}{2} \frac{T_{m,n}-T_{m,n+1}}{\Delta y}-K \frac{\Delta x}{2} \frac{T_{m,n}-T_{m,n-1}}{\Delta y}= h\Delta y (T_{m,n} -T_{\infty}) Setting \Delta y=\Delta x  and simplifying the equation leads to 3.0) T_{m,n} (\frac{h \Delta x}{k}+2)-\frac{h \Delta x}{k} T_{\infty}-\frac{1}{2}(2T_{m-1,n}+T_{m,n+1}+T_{m,n-1})=0

griddedHeatsink3Cropped

Convective heat transfer at the corner node is slightly different from the equation derived above. We move the node to the corner in this figure. -K \frac{ \Delta y}{2} \frac{T_{m,n}-T_{m-1,n}}{\Delta x}-K \frac{ \Delta x}{2} \frac{T_{m,n}-T_{m,n-1}}{\Delta x}= h\Delta y (T_{m,n} -T_{\infty})+h\Delta x (T_{m,n} -T_{\infty}) setting \Delta x =\Delta y and simplifying the above equation results in 4.0) 2T_{m,n}(\frac{h \Delta x}{k}+1)-2 \frac{h \Delta x}{k}T_{\infty}-(T_{m-1,n}+T_{m,n-1})=0

Heat generation at the boundary

griddedHeatsink4Cropped

The heat generated by the transistor \dot{q} at node (m,n) is in \frac{power}{length} unit. The heat flux at the mating surface of the transistor with area  W_{1} *H_{1} is in units  \frac{watt}{{meter}^2}. This figure shows the grid super imposed on the image of the system to be analyzed. Nine nodes (1,1) through (3,3) are used for brevity. The energy balance at node  (m,n) is given by q_{(m-1,n)->(m,n)}+q_{(m+1,n)->(m,n)}+q_{(m,n+1)->(m,n)}+q_{trans}+\frac{dq_{m,n}}{dt}=0 At steady state \frac{dq_{m,n}}{dt}=0 and the heat transferred to the system by the transistor is given by q_{trans}=\frac{q_{tr}*\Delta x}{k}. The heat transfer equations for q_{(m-1,n)->(m,n)} and q_{(m+1,n)->(m,n)} are q_{(m-1,n)->(m,n)}=k\frac{\Delta y}{2}*\frac{T_{m-1,n}-T_{m,n}}{\Delta x} q_{(m+1,n)->(m,n)}=k\frac{\Delta y}{2}*\frac{T_{m+1,n}-T_{m,n}}{\Delta x} For node q_{(m,n+1)->(m,n), the heat transfer equation is q_{(m+1,n)->(m,n)}=k \Delta x *\frac{T_{m,n+1}-T_{m,n}}{\Delta y} substituting these equations into the energy balance equation above and simplifying the equation results in \frac{q\Delta x}{k}-2k T_{m,n}+kT_{m,n+1}+\frac{1}{2}k T_{m+1,n}+\frac{1}{2}T_{m-1,n}=0

# maple code for algebraic manipulation</pre>
# restart; # NULL; eq1 := q[1]+q[2]+q[3]+q[4]+q[5] = 0;
q[1] := (1/2)*k*dy*(T[m-1, n]-T[m, n])/dx;
q[2] := (1/2)*k*dy*(T[m+1, n]-T[m, n])/dx;
q[3] := k*(dx*1)*(T[m-1, n]-T[m, n])/dy;
q[4] := q[gen]*dx/k; q[5] := 0; dx := dy;
2*collect(eq1, T[m, n]); 

Thermal model for finite difference method using conduction , convection and heat generation equations we develop a thermal model for finite difference method solution below.

% This is an example of two dimensional heat distribution
% on a square plate. For illustration,
% the number of nodes are limited 9.Reference [1] and [4] used.
% The figure below shows a square plate with the
% nodes at sides and top exposed to convective heat transfer.
% Heat enters the system at node T_3,2
%
%
%                     T_inf
%           o(T_1,1)----o(T_1,2)----o(T_1,3)
%           |           |           |
%           |           |           |
%           |           |           |
%           |           |           |
%           o(T_2,1)----o(T_2,2)----o(T_2,3)
%           |           |           |
%           |           |           |
%           |           |           |
%           |           |           |
%           o(T_3,1)----o(T_3,2)----o(T_3,3)
%                     q in
% Nodel equations for heat transfer of the system above
%
%
% clear all previous variables
clear all
Tinf=20;
h=10;k=200;dx=0.025
g=((h*dx)/k)+2;
g1=2*(((h*dx)/k)+1);
g3=((h*dx)/k)*Tinf;
q=100;
%clear the screen
clc
% load matrix with coefienct of the equations
% T_1,1 T_1,2 T_1,3 T_2,1 T_2,2 T_2,3 T_3,1 T_3,2 T_3,3
M= [g1      -1       0      -1       0      0       0       0       0;
    -0.5    g       -0.5    0       -1      0       0       0       0;
    0       -1      g1      0        0      -1      0       0       0;
    -0.5    0       0       g        -1     0       -0.5    0       0;
    0       1       0       1       -4      1       0       1       0;
    0       0       -0.5    0        -1     g       0       0      -.5;
    0       0       0       -1       0      0       g1      -1      0;
    0       0       0       0        2      0       1       -4      1;
    0       0       0       0        0     -1       0       -1      g1];

% load vector B with the right hand side of the equations
B=[2*g3;
    g3;
    2*g3;
    g3;
    0;
    g3;
    2*g3;
    (-2*q*dx)/k;
    2*g3;];
% multiply the inverse of M with B to solve
%for variables T1 though T4
%X=inv(M)*B;
X=inv(M)*B

% X =
%
%    21.4264
%    21.4276
%    21.4264
%    21.4287
%    21.4305
%    21.4287
%    21.4311
%    21.4371
%    21.4311
x(1,1:3)=X(1:3);
x(2,1:3)=X(4:6);
x(3,1:3)=X(7:9);
% standard interpolation is used for visualization only
B = imresize(x, 1000);
image(B,'CDataMapping','scaled')
% users may be ininterested in the following commands
% for further analisys
% C =  rand(4,2);
% image(,'CDataMapping','scaled')
%[x,y] = meshgrid(0:.2:2);
%Z = X.
% [X,Y] = meshgrid(0:.2:2);
% Z = X.*exp(-X.^2 - Y.^2);
% [DX,DY] = gradient(Z,.2,.2);
% figure
% contour(X,Y,Z)
% hold on
% quiver(X,Y,DX,DY)
% hold off
Matlab generated thermal image
Matlab generated thermal image
# maple code for algebraic manipulation</pre></pre>
restart;

NULL;
q[a] := (1/2)*h*dx*(T[inf]-T[1, 1])+(1/2)*h*dy*(T[inf]-T[1, 1]);
q[2*a] := (1/2)*k[1]*dy*(T[1, 2]-T[1, 1])/dx;
q[3*a] := (1/2)*k[1]*dx*(T[2, 1]-T[1, 1])/dy;
eq[1, 1] := q[a]+q[2*a]+q[3*a] = 0;
NULL;

q[b] := (1/2)*k[1]*dy*(T[1, 1]-T[1, 2])/dx; q[3*b] := (1/2)*k[1]*dy*(T[1, 3]-T[1, 2])/dx; q[2*b] := k[1]*dx*(T[2, 2]-T[1, 2])/dy; q[4*b] := h*dx*(T[inf]-T[1, 2]);
eq[1, 2] := q[b]+q[2*b]+q[3*b]+q[4*b] = 0;
NULL;

q := (1/2)*h*dx*(T[inf]-T[1, 3])+(1/2)*h*dy*(T[inf]-T[1, 3]); q[2*c] := (1/2)*k[1]*dy*(T[1, 2]-T[1, 3])/dx; q[3*c] := (1/2)*k[1]*dx*(T[1, 2]-T[1, 3])/dy;
eq[1, 3] := q+q[2*c]+q[3*c] = 0;
NULL;

q[d] := h*dy*(T[inf]-T[2, 1]); q[2*d] := (1/2)*k[2]*dx*(T[3, 1]-T[2, 1])/dy; q[3*d] := (1/2)*k[2]*dy*(T[2, 2]-T[2, 1])/dx+(1/2)*k[1]*dy*(T[2, 2]-T[2, 1])/dx; q[4*d] := (1/2)*k[1]*dx*(T[1, 1]-T[2, 1])/dy;
eq[2, 1] := q[d]+q[2*d]+q[3*d]+q[4*d] = 0;

NULL;

q[e] := (1/2)*k[2]*dy*(T[2, 1]-T[2, 2])/dx+(1/2)*k[1]*dy*(T[2, 1]-T[2, 2])/dx; q[3*e] := (1/2)*k[2]*dy*(T[2, 3]-T[2, 2])/dx+(1/2)*k[1]*dy*(T[2, 3]-T[2, 2])/dx; q[2*e] := k[2]*dx*(T[3, 2]-T[2, 2])/dy; q[4*e] := k[1]*dx*(T[1, 2]-T[2, 2])/dy;
eq[2, 2] := q[e]+q[2*e]+q[3*e]+q[4*e] = 0;

NULL;

q[f] := (1/2)*k[1]*dy*(T[2, 2]-T[2, 3])/dx+(1/2)*k[2]*dy*(T[2, 2]-T[2, 3])/dx;
q[2*f] := (1/2)*k[2]*dx*(T[3, 3]-T[2, 3])/dy; q[4*f] := (1/2)*k[1]*dx*(T[1, 3]-T[2, 3])/dy; q[3*f] := h*dy*(T[inf]-T[2, 3]);
eq[2, 3] := q[f]+q[2*f]+q[3*f]+q[4*f] = 0;
#
#
#
#
#
NULL;

q[g] := (1/2)*h*dy*(T[inf]-T[3, 1]); q[2*g] := (1/2)*h*dx*(T[inf]-T[3, 1]); q[3*g] := k[2]*dy*(T[3, 2]-T[3, 1])/dx; q[4*g] := k[2]*dx*(T[2, 1]-T[3, 1])/dy;
eq[3, 1] := q[g]+q[2*g]+q[3*g]+q[4*g] = 0;
NULL;
q[h] := k[2]*dy*(T[3, 2]-T[3, 3])/dx; q[2*h] := (1/2)*h*dx*(T[inf]-T[3, 3]); q[3*h] := (1/2)*h*dy*(T[inf]-T[3, 3]); q[4*h] := k[2]*dx*(T[2, 3]-T[3, 3])/dy;
eq[3, 3] := q[h]+q[2*h]+q[3*h]+q[4*h] = 0;
#
#
NULL;

#
eq[3, 2] := q[i]+q[2*i]+q[3*i]+q[4*i] = 0;

eq[1, 1];
eq1 := subs({T[1, 1] = T1, T[2, 1] = T2}, eq[1, 1]);
NULL;#
#
A := LinearAlgebra[GenerateMatrix]([eq[1, 1], eq[1, 2], eq[1, 3], eq[2, 1], eq[2, 2], eq[2, 3], eq[3, 1], eq[3, 2], eq[3, 3]], [T[1, 1], T[1, 2], T[1, 3], T[2, 1], T[2, 2], T[2, 3], T[3, 1], T[3, 2], T[3, 3]], augmented = true);
#
#
A2 := subs({h*dx = m2, h*dx = m7, h*dy = m3, h*dx*T[inf] = m9, h*dy*T[inf] = m8, k[1]*dx/dy = m1, k[1]*dy/dx = m4, k[2]*dx/dy = m6, k[2]*dy/dx = m5}, A);
with(CodeGeneration);
A3 := Matlab(A2);

#
<pre> 
Matlab generated thermal image
Matlab generated thermal image
Matlab code
% This is an example of two dimensional heat distribution
% on a square plate. For illustration,
% the number of nodes are limited 9.Reference [1] and [4] used.
% The figure below shows a square plate with the
% nodes at sides and top exposed to convective heat transfer.
% Heat enters the system at node T_3,2
%
%
%                     T_inf
%           o(T_1,1)----o(T_1,2)----o(T_1,3)
%           |           |           |
%           |           |           |
%           |           |           |
%           |           |           |
%           o(T_2,1)----o(T_2,2)----o(T_2,3)
%           |           |           |
%           |           |           |
%           |           |           |
%           |           |           |
%           o(T_3,1)----o(T_3,2)----o(T_3,3)
%                     q in
% Nodel equations for heat transfer of the system above
%
%
% clear all previous variables
% clear all
clc
clear all
%set the boundary conditions and calculate some of the coeficients used in
%the matrix. See maple code.
Tinf=20;
h=10;k(1)=200;dx=0.025;dy=0.025;k(2)=500;
m1=k(1)*(dx/dy);
m2=h*dx;
m3=h*dy;
m4=k(1)*(dy/dx);
m5=k(2)*(dy/dx);
m6=k(2)*(dx/dy);
m7=h*dx;
m8=h*dy*Tinf;
m9=h*dx*Tinf;
tr=1;
q(tr)=100;
% load matrix with coefienct of the equations
% T_1,1 T_1,2 T_1,3 T_2,1 T_2,2 T_2,3 T_3,1 T_3,2 T_3,3
A=[
    -m3 / 0.2e1 - m2 / 0.2e1 - m1 / 0.2e1 - m4 / 0.2e1 m4 / 0.2e1 0 m1 / 0.2e1 0 0 0 0 0 -m9 / 0.2e1 - m8 / 0.2e1;
    m4 / 0.2e1 -m1 - m4 m4 / 0.2e1 0 m1 0 0 0 0 0;
    0 m1 / 0.2e1 + m4 / 0.2e1 -m3 / 0.2e1 - m2 / 0.2e1 - m1 / 0.2e1 - m4 / 0.2e1 0 0 0 0 0 0 -m9 / 0.2e1 - m8 / 0.2e1;
    m1 / 0.2e1 0 0 -m6 / 0.2e1 - m3 - m4 / 0.2e1 - m1 / 0.2e1 - m5 / 0.2e1 m5 / 0.2e1 + m4 / 0.2e1 0 m6 / 0.2e1 0 0 -m8;
    0 m1 0 m5 / 0.2e1 + m4 / 0.2e1 -m5 - m4 - m1 - m6 m5 / 0.2e1 + m4 / 0.2e1 0 m6 0 0;
    0 0 m1 / 0.2e1 0 m5 / 0.2e1 + m4 / 0.2e1 -m6 / 0.2e1 - m3 - m4 / 0.2e1 - m1 / 0.2e1 - m5 / 0.2e1 0 0 m6 / 0.2e1 -m8;
    0 0 0 m6 0 0 -m3 / 0.2e1 - m5 - m6 - m2 / 0.2e1 m5 0 -m9 / 0.2e1 - m8 / 0.2e1;
    0 0 0 0 m6 0 m5 -0.2e1 * m5 - m6 m5 -q(tr) * dx / k(2);
    0 0 0 0 0 m6 0 m5 -m3 / 0.2e1 - m5 - m6 - m2 / 0.2e1 -m9 / 0.2e1 - m8 / 0.2e1;
    ];

% % split matrix with right hand side and left hand side equations
% left hand side of equations
A2(1:9,1:9)=A(1:9,1:9);
% rigth hand side
B=A(1:9,10);
% calculate unknown temperature
X=inv(A2)*B;
% creat image of calculated temperature
x(1,1:3)=X(1:3);
x(2,1:3)=X(4:6);
x(3,1:3)=X(7:9);
% resize image to improve visualization
B = imresize(x, 100);
C=image(B,'CDataMapping','scaled')
%imwrite(C,strcat('img\',num2str(j),'.png'))
%

References: [1] Heat transfer  Tenth edition J.P. Holman [2] Advanced Engineering Mathematics 10th edition Erwin Kreyszeg [3] Heat Transfer Essentials Third Edition Latif M. Jiji [4] Elementary Mathematical and Computational Tools for Electrical and Computer Engineers Using Matlab. Jamal T. Manassah