Image Source: Allied Market Research 

Numerical Analysis of Heat Transfer within a Cooling Fin

Objective

The goal of this analysis was to practice using numerical methods as solutions to differential equations, in order to prepare for situations where the analytical solutions are not known. 

Outcome

Finite-difference schemes and assumed boundary conditions were applied to convective temperature and heat transfer equations, then directly compared to analytical solutions. A table of errors was generated for interpreting accuracy. 


Findings

This analysis includes four figures and two tables. Below the figures and tables is the python code used to calculate the figures. 

It is clear that the temperature errors are much smaller than the heat transfer errors, because calculating the heat transfer from the temperature slope numerically will cause the errors to greatly expound. To increase accuracy, fourth-order finite difference schemes could be implemented. Increasing the number of nodes and modifying the step size would also cause results to converge better.

Figure 1: Shown is the temperature vs the position of a cooling fin modeled by an analytical solution and a second-order finite difference method. This model assumes that the ends of the fin are held at fixed temperatures.

Figure 2: The heat transfer vs the position of a cooling fin, also modeled by an analytical solution and a second-order finite difference method. The ends of the cooling fin are also held at fixed temperatures.

Figure 3: Shown is the temperature vs the position of a cooling fin modeled by an analytical solution and a second-order finite difference method. This time, the front end is held at a constant while the end of the fin is assumed to be adiabatic (zero heat transfer).

Figure 4: The Heat Transfer vs the position of a cooling fin modeled by an analytical solution and a second-order finite difference method. This model is additionally modeled with an adiabatic tip and a constant-temperature beginning

Python Code

'''

Heat Transfer of a Cooling Fins; Numerical vs

Analytical Models


Dallin Carter


The purpose of this code is to compare the difference

of temperature and heat transfer values across cooling

fins between numerical and analytical solutions.


'''


#Libraries

import numpy as np

import matplotlib.pyplot as plt

from scipy.linalg import solve


#Physical Conditions

length_fin = 0.15 #m

area_fin = 1.67E-4 #m^2

length_perimeter = 0.066 #m

k = 16 #W/m-K

h = 25 #W/m^2-K

temp_base = 500 #K

temp_freestream = 300 #K

temp_tip = 350 #K


#Numerical Conditions

dx = 0.03 #m

m = np.sqrt(h*length_perimeter/(k*area_fin))


'''

Analytical Functions:


The input is distance (in meters) across the fin, and returns the value

of Temperature (in Kelvin) or Heat Transfer (in Watts). Also included is

a linspace of xx, which has calculates the temperature/heat for many values

across the fin.

'''

def T1_ana(x):

    return ((temp_tip-temp_freestream)*np.sinh(m*x)+(temp_base-temp_freestream)*np.sinh(m*(length_fin-x)))/np.sinh(m*length_fin) + temp_freestream


def q1_ana(x):

    return (m*k*area_fin/np.sinh(m*length_fin))*((temp_base-temp_freestream)*np.cosh(m*(length_fin-x))-(temp_tip-temp_freestream)*np.cosh(m*x))


def T2_ana(x):

    return (np.cosh(m*(length_fin-x))/np.cosh(m*length_fin))*(temp_base-temp_freestream)+temp_freestream


def q2_ana(x):

    return m*k*area_fin*(np.sinh(m*(length_fin-x))/np.cosh(m*length_fin))*(temp_base-temp_freestream)


xx = np.linspace(0,length_fin,100)



'''First Boundary Conditions'''

n = 6

A = np.zeros([n,n])

b = np.zeros(n)

x = np.linspace(0,length_fin,n)


#Boundary Conditions for Temperature

A[0,0] = 1

A[-1,-1] = 1

b[0] = temp_base

b[-1] = temp_tip


# Fill in matrix

for i in range(1,n-1):

    A[i,i-1] = 1

    A[i,i] = -2-dx**2*m**2

    A[i,i+1] = 1

    b[i] = -dx**2*m**2*temp_freestream


# Solve for Temperature

T1 = solve(A,b)


#Plotting Temp vs Position

plt.plot(x,T1,marker='o', linestyle = ':', label = 'Numerical Solution')

plt.plot(xx,T1_ana(xx), label = 'Analytical Solution', color='k')

plt.xlabel('Position (m)')

plt.ylabel('Temperature (K)')

plt.legend()

plt.savefig('HeatTransfer1')

plt.show()


#Solve for Heat Transfer

q1 = np.zeros([n])

q1[0] = -k*area_fin*(-T1[2]+4*T1[1]-3*T1[0])/(2*dx)

q1[n-1] = -k*area_fin*(T1[-3]-4*T1[-2]+3*T1[-1])/(2*dx)

for i in range(1, n-1):

    q1[i] = -k*area_fin*(-T1[i-1]+T1[i+1])/(2*dx)


#Plotting Heat Transfer vs Position

plt.plot(x,q1,marker='o', linestyle = ':', label = 'Numerical Solution')

plt.plot(xx,q1_ana(xx), label = 'Analytical Solution', color='k')

plt.xlabel('Position (m)')

plt.ylabel('Heat Transfer (W)')

plt.legend()

plt.savefig('HeatTransfer2')

plt.show()



'''Second Boundary Conditions'''

#Boundary Conditions for Temperature

A[0,0] = 1

A[-1,-3] = 1

A[-1,-2] = -4

A[-1,-1] = 3

b[0] = temp_base

b[-2] = 0

b[-1] = 0


# Fill in matrix

for i in range(1,n-1):

    A[i,i-1] = 1

    A[i,i] = -2-dx**2*m**2

    A[i,i+1] = 1

    b[i] = -dx**2*m**2*temp_freestream


# Solve for Temperature

T2 = solve(A,b)


#Plot Temp vs Position

plt.plot(x,T2,marker='o', linestyle = ':', label = 'Numerical Solution')

plt.plot(xx,T2_ana(xx), label = 'Analytical Solution', color='k')

plt.xlabel('Position (m)')

plt.ylabel('Temperature (K)')

plt.legend()

plt.savefig('HeatTransfer3')

plt.show()


#Solve for Heat Transfer

q2 = np.zeros([n])

q2[0] = -k*area_fin*(-T2[2]+4*T2[1]-3*T2[0])/(2*dx)

q2[n-1] = -k*area_fin*(T2[-3]-4*T2[-2]+3*T2[-1])/(2*dx)

for i in range(1, n-1):

    q2[i] = -k*area_fin*(-T2[i-1]+T2[i+1])/(2*dx)


#Plot of Heat Transfer vs Position

plt.plot(x,q2,marker='o', linestyle = ':', label = 'Numerical Solution')

plt.plot(xx,q2_ana(xx), label = 'Analytical Solution',color='k')

plt.xlabel('Position (m)')

plt.ylabel('Heat Transfer (W)')

plt.legend()

plt.savefig('HeatTransfer4')

plt.show()


'''

Calculating Error and Printing the Solutions

'''

#Calculating the Errors

ErrorT1 = abs(T1 - T1_ana(x))/T1_ana(x)*100

ErrorT2 = abs(T2 - T2_ana(x))/T2_ana(x)*100

Errorq1 = abs(q1 - q1_ana(x))/q1_ana(x)*100

Errorq2 = abs(q2 - q2_ana(x))/q2_ana(x)*100


#Printing Solutions at all Nodes

#First LaTeX Table Print

for i in range(n):

    print(f"{x[i]:.3f} & {T1[i]:.2f} & {ErrorT1[i]:.2f} & {q1[i]:.2f} & {Errorq1[i]:.2f} \\\\")

print()


#Second LaTeX Table Print

for i in range(n):

    print(f"{x[i]:.3f} & {T2[i]:.2f} & {ErrorT2[i]:.2f} & {q2[i]:.2f} & {Errorq2[i]:.2f} \\\\")