Project 4: Heat Distribution to a Cooling Fin

Mark Schools, Sarah Polzer, Joshua Cho

Introduction

Heat sinks are used to move excess heat away from the point where it is generated. In this project, we modeled the steady state distribution along a rectangular fin of a heat sink. The heat energy entered the fin along part of one side. The goal was to design the fin to keep the temperature within safe tolerances.

The cooling fin is a rectangular slab with dimensions $Lx × Ly$ and width $\delta$ cm where $\delta$ is small. Temperature will be denoted by $u(x,j)$ and will be considered constant along the width dimension.

There are three different ways that heat moves: conduction, convection, and radiation. Conduction refers to the transfer of energy between neighboring molecules, convection refers to the molecules themselves moving, and radiation is the movement of energy through photons. We will not consider radiation.

Conduction through a material can be modeled with Fourier's first law; q is the heat energy per unit time in watts, A is the area of the material, ∇u is the temperature gradient, and K is the thermal conductivity of the material.

$q = -kA∇u$

Convection can be modeled with Newton's law of cooling; H is the convective heat transfer coefficient, $u_b$ is the ambient temperature of the surrounding fluid.

$q=-HA(u-u_b)$

The cooling fin is a rectangle $[0,L_x] × [0, L_y] × \delta$ . The energy entering the box per unit time equals the energy leaving. The heat flux into the box through the $\Delta y × \delta$ sides and $\Delta x × \delta$ sides is by conduction, and through the two $\Delta x × \Delta y$ sides is by convection, and through the two $\Delta x × \Delta y$ sides is by convection, which yields the steady state equation

$-K \Delta y \delta u_x(x,y)+K\Delta y \delta u_x (x + \Delta x ,y) - K\Delta x \delta u_y (x,y) +K\Delta x \delta u_y (x, y + \Delta y)- 2H \Delta x \Delta yu(x,y) = 0$

Based on the text, the convective boundary condition is being implied

$K u_{normal} = H u$

where $u_{normal}$ is the partial derivative in respect to the normal direction $\vec{n}$ . This boundary condition is known as the $\textbf{Robin}$ boundary condition, which involves both the function value and its derivative. We also assume that power enters the fin along one side according to Fouier's Law:

$u_{normal} = \frac{P}{L \delta K}$

Where P is total power and L is length of input.

the finite difference approximations with step sizes $\textit{h}$ and $\textit{k}$ are used to approximate the PDE as

$\frac{u_{i+1}-2u_{ij}+u_{i-1,j}}{h^2}+\frac{u_{i,j+1}-2u_{ij}+u_{i,j-1}}{k^2} = \frac{2H}{K \delta} u_{ij}$

Which is used for the interior points $(x_i,y_j)$ . The fin obeys the Robin conditions using the first derivative approximation

$f'(x) = \frac{-3f(x) +4f(x + h) - f(x + 2h)}{2h} + O(h^2)$

which to apply to the fin edges the normal direction translates to the bottom, top , left and right edges.

Assuming the power enters the left side of the fin, Fouier's law leads to the equation

$\frac{-3u_{1j}+4u{2j}-u_{3j}}{2h} = -\frac{P}{L \delta K}$

where 1 $\leq$ i $\leq$ m, and 1 $\leq$ j $\leq$ n

It is assumed that the fin is composed of aluminum, whose thermal conductivity is $K = 1.68 W/cm C$ . It is also assumed that the convective heat transfer coefficient is $H=0.0005 W/cm^2$ and that room temperature is $u_b = 20 C$

Question 1

Begin with a fin of dimensions 2 × 2 cm, with 1 mm thickness. Assume that 5W of power is input along the entire left edge, as if the fin were attached to dissipate power from a CPU chip with L = 2 cm side length. Solve the PDE (8.44) with M = N = 10 steps in the x and y directions. Use the mesh command to plot the resulting heat distribution over the xy-plane. What is the maximum temperature of the fin, in ◦C ?

To complete this problem - we modified poisson.m from the textbook.


x
% Program 8.5 Finite difference solver for 2D Poisson equation
% with Dirichlet boundary conditions on a rectangle
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps
% Output: matrix w holding solution values
% Example usage: w=poisson(0,1,1,2,4,4)
function w=poisson(xl,xr,yb,yt,M,N,Delta,L,P)
f=@(x,y) 0; % define input function data
%g1=@(x) log(x.^2+1); % define boundary values
%g2=@(x) log(x.^2+4); % Example 8.8 is shown
%g3=@(y) 2*log(y);
%g4=@(y) log(y.^2+1);
m=M+1;n=N+1; mn=m*n;
H=0.005; K=1.68;
h=(xr-xl)/M;h2=h^2;k=(yt-yb)/N;k2=k^2;
x=xl+(0:M)*h; % set mesh values
y=yb+(0:N)*k;
A=zeros(mn,mn);b=zeros(mn,1);
for i=2:m-1 % interior points
for j=2:n-1
A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;
A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-(2.0*H/(K*Delta));
A(i+(j-1)*m,i+(j-2)*m)=1/k2;A(i+(j-1)*m,i+j*m)=1/k2;
b(i+(j-1)*m)=f(x(i),y(j));
end
end
for i=1:m % bottom and top boundary points
j=1;A(i+(j-1)*m,i)=-(3/(2*k))+H/K;A(i+(j-1)*m,i+m)=2/k;A(i+(j-1)*m,i+2*m)=-1/(2*k); 
j=n;A(i+(j-1)*m,i+(n-1)*m)=-3/(2*k)+H/K;A(i+(j-1)*m,i+(n-2)*m)=2/(k);A(i+(j-1)*m,i+(n-3)*m)=-1/(2*k);
end
for j=2:n-1 % left and right boundary points
    i=1;A(i+(j-1)*m,1+(j-1)*m)=-3/(2*h);A(i+(j-1)*m,2+(j-1)*m)=2/h;A(i+(j-1)*m,3+(j-1)*m)=-1/(2*h);b(i+(j-1)*m)=-P/(L*Delta*K);
    i=1;A(i+(j-1)*m,1+(j-1)*m)=-3/(2*h)+H/K;A(i+(j-1)*m,2+(j-1)*m)=2/h;A(i+(j-1)*m,3+(j-1)*m)=-1/(2*h);
end
i=m;A(i+(j-1)*m,m+(j-1)*m)=-3/(2*h)+H/K;A(i+(j-1)*m,m-1+(j-1)*m)=2/h;A(i+(j-1)*m,m-2+(j-1)*m)=-1/(2*h);
end
v=A\b; % solve for solution in v labeling
w=reshape(v(1:mn),m,n); %translate from v to w
mesh(x,y,w)

We called poisson.m with the code below


xxxxxxxxxx
poisson(0,2,0,2,10,10,0.1,2,5.0)

Results:

The maximum temperature that we obtained was 164.9357

Question 2

Increase the size of the fin to 4 × 4 cm. Input 5W of power along the interval [0,2] on the left side of the fin, as in the previous step. Plot the resulting distribution, and find the maximum temperature. Experiment with increased values of M and N. How much does the solution change?

Once again we used poisson.m to complete this problem and made the following modifications

We added an if-else to poisson.m


xxxxxxxxxx
for j=2:n-1 % left and right boundary points
if j<=n/2
    i=1;A(i+(j-1)*m,1+(j-1)*m)=-3/(2*h);A(i+(j-1)*m,2+(j-1)*m)=2/h;A(i+(j-1)*m,3+(j-1)*m)=-1/(2*h);b(i+(j-1)*m)=-P/(L*Delta*K);
else
    i=1;A(i+(j-1)*m,1+(j-1)*m)=-3/(2*h)+H/K;A(i+(j-1)*m,2+(j-1)*m)=2/h;A(i+(j-1)*m,3+(j-1)*m)=-1/(2*h);
end

and we changed the way that we called poisson.m


xxxxxxxxxx
poisson(0,4,0,4,10,10,0.1,2,5.0)

Results:

The maximum temperature that we obtained was 69.4704.

We changed the values of M and N to 2,2; 7,7; and 20,20 by replacing M and N when calling poisson.m

For M and N at 2,2

For M and N at 2 - the maximum temperature was 0

for M and N at 7,7

for M and N at 7 - the maximum temperature was 55.1246

for M and N at 20,20

for M and N at 20 - the maximum temperature was 41.0257

Question 3

Find the maximum power that can be dissipated by a 4 × 4 cm fin while keeping the maximum temperature less than 80◦C. Assume that the bulk temperature is 20◦C and the power input is along 2 cm, as in steps 1 and 2.

We completed this problem by using the code below while calling the poisson function.


xxxxxxxxxx
P = 1;
w = 1;
%poisson(xl,xr,yb,yt,M,N,Delta,L,P)
%P = 5.0
while(max(w) < 60)
w=poisson(0,4,0,4,10,10,0.1,2,P);
P = P +.001;
end
P = P -.001;
w=poisson(0,4,0,4,10,10,0.1,2,P)
P

The maximum temperature that we obtained was 80.0076 degrees Celsius

The maximum power that we obtained was 6.0650

Question 4

Replace the aluminum fin with a copper fin, with thermal conductivity K = 3.85W/cm ◦C. Find the maximum power that can be dissipated by a 4 × 4 cm fin with the 2 cm power input in the optimal placement, while keeping the maximum temperature below 80 ◦C.

To complete this problem we made the following change to the code in poisson.m


xxxxxxxxxx
K=3.85;

We called the function using the same code in question 3

The maximum temperature that we obtained was 80.0029 The maximum power that we obtained 7.8120

Question 5

Plot the maximum power that can be dissipated in step 4 (keeping maximum temperature below 80 degrees) as a function of thermal conductivity, for 1 ≤ K ≤ 5 W/cm◦C.

To create our plot we wrote the code below and slightly modified poisson.m to make K a parameter


xxxxxxxxxx
K_vals = zeros(1,40);
P_vals = zeros(1,40);
count = 1;
for K=1:0.1:5
    while(max(w) < 60)
    w=poisson(0,4,0,4,10,10,0.1,2,P,K);
    P = P +.001;
    end
        P = P -.001;
    w=poisson(0,4,0,4,10,10,0.1,2,P,K);
    K_vals(1,count) = K;
    P_vals(1,count) = P;
    count = count + 1;
end
plot(K_vals,P_vals)

Here is the modification of poisson.m


xxxxxxxxxx
function w=poisson(xl,xr,yb,yt,M,N,Delta,L,P,K)


xxxxxxxxxx
H=0.005; %K=3.85;

Question 6

Redo step 4 for a water-cooled fin. Assume that water has a convective heat transfer coefficient of H = 0.1 W/cm2 ◦C, and that the ambient water temperature is maintained at 20◦C

We made the following change to poisson.m to make the fin water-cooled.


xxxxxxxxxx
H=0.1

We also modified the script to maintain the temperature at 20 degrees C

The maximum temperature that we obtained was 80.0002

The maximum power that we obtained was 36.2940