from numpy import matrix
from numpy import linalg
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.1 "
#Cross section of the element in m is given as
b = 0.1; #breadth in m
H = 0.01; #height in m
#Temperature of surrrounding oil in C is given as
Tinfinity = 80;
#Correspoding heat transfer coefficient in W/m2-K is given as:
h = 42.0;
#Heat generation rate is given in W/m3 as
qg = 10**6;
#Temperature below which element needed to maintain in C is
T = 200.0;
# Thermal conductivity of iron in W/m-K is taken as
k = 64.0;
#Because of symmetry we need to consider only half of the thickness of the heating element
L = H/2.0; #Length in m
#We are defining five nodes at a distance of (i-1)*dx, where i=1,2,3,4,5
N = 5.0; #Total number of grid points
dx = L/(N-1); #dx in m
#Since no heat flows across the top face, it corresponds to a zero-heat
#flux boundary condition.
#Applying Eq. (2.1) to a control volume extending from x=L-dx/2 to x=L
#We get TN=TN-1 +qg*dx*dx/(2*k)
#At the left face, , we have a surface convection boundary condition to which Eq. (3.7) can be applied
#Determining all the matrix coefficients in Eq. (3.11)
a1 = 1; #Matrix coefficient a1 in SI units
b1 = 1/(1+(h*dx)/k); #Matrix coefficient b1 in SI units
c1 = 0; #Matrix coefficient c1 in SI units
d1 = (dx/k)*((h*Tinfinity+(qg*dx)/2)/(1+(h*dx)/k)); #Matrix coefficient d1 in SI units
a2 = 2;a3 = a2;a4 = a3;#Matrix coefficient a2 in SI units
b2 = 1;b3 = b2;b4 = b3;#Matrix coefficient b2 in SI units
c2 = 1;c3 = c2;c4 = c3;#Matrix coefficient c2 in SI units
d2 = ((dx*dx)*qg)/k;d3 = d2;d4 = d2;#Matrix coefficient d2 in SI units
a5 = 1;b5 = 0;c5 = 1;d5 = ((dx*dx)*qg)/(2*k);#Matrix coefficient a5 in SI units
#Umath.sing the algorithm given in Appendix 3 for solving the tridiagonal system, we find the temperature distribution given as:
print "Final temperature distribution in C is the following"
#From equation 3.11
#Matrix A in the Appendix 3
A = [[a1,-b1,0,0,0],[-c2,a2,-b2,0,0],[0,-c3,a3,-b3,0],[0,0,-c4,a4,-b4],[0,0,0,-c5,a5]]
#Matrix D in the Appendix 3
D = [[d1],[d2],[d3],[d4],[d5]];
#Temperature matrix where temp are in degree C as given by appnedix 3
T = ((linalg.inv(A))*D)
z1=0
z2=0
z3=0
z4=0
z5=0
for i in range(0,5):
z1=z1+T[i][0]
z2=z2+T[i][1]
z3=z3+T[i][2]
z4=z4+T[i][3]
z5=z5+T[i][4]
print round(z1,4),"\n",round(z2,4),"\n",round(z3,4),"\n",round(z4,4),"\n",round(z5,4)
# the answer in the book is slightly different due to approximation
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.2 "
# we have to determine minimum depth xm at which a water main must be buried to avoid freezing
#Initial temperature of soil in C is given as:
Ts = 20.0;
# Under the worst conditions anticipated it would be subjected to a surface
# temperature of -15C for a period of 60 days
#Max temperature in degree C
Tmax = -15.0;
#Time period in days
dt = 60.0;
#We will use the following properties for soil (at 300 K)
rho = 2050;#density in kg/m3
k = 0.52;#thermal conductivity in W/m-K
c = 1840;#specific heat in J/kg-K
alpha = 0.138*(10**(-6));#diffusivity in m2/sec
#Fourier number is defined as:
#Fo=dt*alpha/(dx*dx);
#Let us select a maximum depth of 6 m
#First, let us choose , giving dx=1.2m
dx = 1.2; #dx in m
dt = (30*24)*3600;#Days converted in seconds
#Temperature array for the old temperature in degree C
Tnew = [-15,20,20,20,20,20];
#Temperature array for the new temperature in degree C
Told = [-15,20,20,20,20,20];
#Fourier number is defined as:
Fo = (dt*alpha)/(dx*dx);
#Umath.sing eq. 3.15
#Initialmath.sing timestep for looping
timestep = 0;
for timestep in range(0,100):
for N in range (2,4):
#New temp in degree C
Tnew[N] = Told[N]+Fo*(Told[N+1]-2*Told[N]+Told[N-1]);
#Incrementing timestep
timestep = timestep+1;
print "With dx=1.2m, we have the following distribution"
#New temp in degree C
Tnew
print "Depth in m at which temperature would be 0 degree C would be"
#Depth in m
xm = (0-Tnew[0]/(Tnew[1]-Tnew[0]))*dx
print xm
# the answer in the textbook is wrong
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.3 "
#initial temperature of the sheet in C is given as:
Tinitial = 500.0;
#thickness of the sheet in m is given as
th = 0.02;
#density in kg/m3 is given for steel as
rho = 8500.0;
#specific heat in J/kg-K is given as
c = 460.0;
#thermal conductivity in W/m-K is given as
k = 20.0;
#The heat transfer coefficient in W/m2-K to the air is given as
h = 80.0;
#the ambient air temperature in degree C is
Tinfinity = 20.0;
#Final temperature required to achieve in C is
Tfinal = 250.0;
#The transient cooling of stainless steel sheet can be modeled as a semi-infinite slab
#because the thickness of the sheet is much smaller than its width and length.
L = th/2.0; #Length in m
#Finding chart solution
#Biot number shall be
Bi = (h*L)/k;
#Since Bi<0.1 and hence the sheet can be treated as a lumped capacitance.
#To use fig. 2.42 on page 135, we need to calculate the following value:
value = (Tfinal-Tinfinity)/(Tinitial-Tinfinity); #value required
#So, now umath.sing fig. 2.42, we have alpha*dt/(L*L)=19
#BY the definition of thermal diffusivity,in SI units we have
alpha = k*1.0/(rho*c);
print "By chart solution, time required in seconds comes out to be"
#time required in seconds
t = ((19.0*L)*L)/alpha
print round(t,2)
#Proceeding to the numerical solution
#consider half the sheet thickness,with x=0 being the math.exposed left face and
#x=L being the sheet center-line
#Umath.sing 20 control volumes
N = 21.0; #Total number of grid points
dx = L/20.0; #dx in m
Told=numpy.zeros((1,21))
Tnew=numpy.zeros((1,21))
#Old temperature array
for N in range(0,20):
#Old temp in degree C
Told[0,N] = Tinitial;
#New temp in degree C
Tnew[0,N] = Tinitial;
#Increment of Time in sec
dt = 5.57;
#Condition of looping
while Told[0,20]>250:
#C1 of governing equation in SI units
C1 = (alpha*dt)/(dx*dx);
#C2 of governing equation in SI units
C2 = ((2*h)*dt)/((rho*c)*dx);
#C3 of governing equation in SI units
C3 = 2*C1;
#New temp in C as given by the equations of finite difference method
Tnew = (Told[0]+C2*(Tinfinity-Told[0])+C3*(Told[1]-Told[0]));
t=t+5.57 # increment
Tnew = Told[0,20]+C3*(Told[0,19]-Told[0,20]);
for N in range(2,20):
#New temp in C as given by the equations of finite difference method
Tnew = t+dt+(Tnew,N,Told(N)+C1*(Told(N+1)-2*Told(N)+Told(N-1)));
#Modified time for new loop
t = t+dt;
# L.67: No simple equivalent, so mtlb_fprintf() is called.
print "As per numerical solution time comes out to be ",round(t,2)," seconds\n"
print "This time is about 1.5% less than the chart solution"
# the solution in the book is slightly different due to approximation
import numpy
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.4 "
#Dimensions of the cross section in inches
l = 1;
b = 1;
#Dividing domain such that there are four nodes in x and y direction
dx = 1/3.0; #dx in inches
dy = 1/3.0; #dy in inches
#Assigning Temperature in C for top and bottom surface
T=numpy.zeros((4,4))
for i in range(0,3):
T[0,i] = 0;
T[3,i] = 0;
#Assigning Temperature in C for side surfaces
for j in range(0,3):
T[j,0] = 50;
T[j,3] = 100;
#Assigning Temperature in C for interior nodes
for i in range(0,2):
for j in range(0,2):
T[i,j] = 0;
#Defining looping parameter
step = 0;
for step in range (0,50):
#Umath.sing governing equations of finite difference
T[2,1] = 0.25*(50+0+T[1,2]+T[2,1]);
T[1,1] = 0.25*(50+0+T[2,1]+T[1,2]);
T[1,2] = 0.25*(100+0+T[2,1]+T[1,2]);
T[2,2] = 0.25*(100+0+T[1,1]+T[2,2]);
#print "At steady state, Final temperature of the cross section in C would be"
#New temp distribution in degree C
print'Temperature T(2,2) in degree C is ',round(T[1,1],2)
print'Temperature T(2,3) in degree C is ',round(T[2,1],2)
print'Temperature T(3,2) in degree C is ',round(T[1,2],2)
print'Temperature T(3,3) in degree C is ',round(T[2,2],2)
import numpy
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.5 "
#Thermal conductivity of alloy bus bar in W/m-K is given as
k = 20;
#Heat generation rate in W/m3 is given as
qg = 10**6;
#dimensions of the bar in m is given as
L = 0.1;#Length in m
b = 0.05;#Width in m
d = 0.01;#Thickness in m
#For top edge, heat transfer coefficient in W/m2K and ambient temperature
#in C are
h = 75;
Tinfinity = 0;
#We are taking a total of 11 nodes in the direction of length and 6 nodes
#in the direction of width
dx = 0.01; #dx in m
dy = 0.01; #dy in m
Told=numpy.zeros((6,12))
Tnew=numpy.zeros((6,12))
#Assigning a guess temperature of 25C to all nodes
for i in range(0,6):
for j in range(0,12):
#Old temp. in degree C
Told[i,j] = 25;
#Assigning temperature on the left and right hand side
for i in range(0,6):
#Old temp. in degree C
Told[i,0] = 40;
Told[i,10] = 10;
#New temp. in degree C
Tnew[i,0] = 40;
Tnew[i,10] = 10;
#Intitalisation of looping parameter
p = 0;
#Iteration to find temperature distribution
while p<500:
#Equation for all interior nodes
for i in range(1,5):
for j in range(1,10):
#New temp. in degree C
Tnew[i,j] = 0.25*(Told[i-1,j]+Told[i+1,j]+Told[i,j-1]+Told[i,j+1]+((qg*dx)*dx)/k);
#Equation for top wall
for j in range(1,10):
#New temp. in degree C
Tnew[0,j] = (h*Tinfinity+(qg*dx)/2+(k*(0.5*(Told[1,j-1]+Told[1,j+1])+Told[1,j]))/dx)/(h+(2*k)/dx);
#Equation for bottom wall
for j in range(1,10):
#New temp. in degree C
Tnew[5,j] = 0.25*(Told[5,j-1]+Told[5,j+1])+0.5*Told[4,j]+((qg*dx)*dx)/(4*k);
for i in range(0,6):
for j in range(0,11):
#Assigning Old Temp=New Temp
Told[i,j] = round(Tnew[i,j],2);
#New looping parameter incremented
p = p+1;
print "The temperature distribution in the bar in C is the following"
#Old temp. in degree C
for i in range(0,11):
print "Node",i+1,"= ",Told[0,i]
#Finding maximum temperature
Tmax = Told[0,0];
for i in range(0,6):
for j in range(0,11):
if Told[i,j]>Tmax:
Tmax = Told[i,j];
else:
Tmax = Tmax;
print "The maximum temperature in C in the alloy bus bar is"
#maximum temperature in C
print Tmax
#Finding heat transfer rate
dz = 0.01; #dz in m
#Defining areas
A=numpy.zeros((1,11))
for i in range(1,11):
A[0,i] = dx*dz; #Area in m2
q=numpy.zeros((1,11))
for i in range(0,11):
#heat transfer rate in W
q[0,i] = round((h*A[0,i])*(Tnew[0,i]-Tinfinity),3);
print "The heat transfer rate from the top edge in W is given by"
#heat transfer rate in W
for i in range(0,11):
print "node",i+1,"= ",q[0,i]
# the answer in the textbook is incorrect in the calculations
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.6 "
#Thermal diffusivity in m2/s
alpha = 0.000008;
#%Thermal conductivity of alloy bus bar in W/m-K is given as
k = 20;
#density*specific heat product in SI units
pc = k/alpha;
#dimensions of the bar in m is given as
L = 0.1;#Length in m
b = 0.05;#Width in m
d = 0.01;#Thickness in m
#Heat generation rate in W/m3 is given as
qg = 10**6;
#Assigning temperature on the left and right hand side
for i in range (0,6): #i is the looping parameter
#Old temp. in degree C
Told[i,0] = 40;
Told[i,10] = 10;
#New temp. in degree C
Tnew[i,0] = 40;
Tnew[i,10] = 10;
#Assigning a guess temperature of 20C to all nodes
for i in range(0,6):#i is the looping parameter
for j in range(0,11):#j is the looping parameter
#Guess temp. in degree C
Told[i,j] = 20;
Tnew[i,j] = 20;
#Initialimath.sing time
m = 0;
#For top edge, heat transfer coefficient in W/m2K and ambient temperature
#in C are
h = 75;
Tinfinity = 0;
#We are taking a total of 11 nodes in the direction of length and 6 nodes
#in the direction of width
dx = 0.01; #dx in m
dy = 0.01; #dy in m
#Largest permissible time step in sec is
tmax = 1/((2*alpha)*(1/(dx*dx)+1/(dy*dy)));
m=1140; # explicit time in secs
#Rounding it off to nearest integer
t = 3; #timestep in seconds
#Condition for convergence
while abs(Tnew[4,5]-Told[4,5])<0.0001:
#Equation for all interior nodes
for i in range(1,5):
for j in range (1,10):
#New temp. in degree C
Tnew[i,j] = (Told[i,j]+(alpha*t)*((Tnew[i+1,j]+Tnew[i-1,j])/(dx*dx)+(Tnew[i,j+1]+Tnew[i,j-1])/(dy*dy)+qg/k))/(1+((2*alpha)*t)*(1/(dx*dx)+1/(dy*dy)));
#Equation for top wall
for j in range (1,10):
#New temp. in degree C
Tnew[0,j] = (Told[0,j]+((2*t)/((dx*dx)*pc))*(k*((Tnew[0,j+1]+Tnew[0,j-1])/2+Tnew[1,j]))+((qg*dx)*dx)/2+(h*dx)*Tinfinity)/(1+((2*t)/((dx*dx)*pc))*(2*k+h*dx));
#Equation for bottom wall
for j in range (1,10):
#New temp. in degree C
Tnew[5,j] = (Told[5,j]+((2*t)/((dx*dx)*pc))*(k*((Tnew[5,j+1])+Tnew[5,j-1])/2+Tnew[4,j]))+(((qg*dx)*dx)/2)/(1+((2*t)/((dx*dx)*pc))*(2*k));
#New time in sec
m = m+t;
print "Time required to reach steady state using explicit method is 1140 seconds"
print "Time required to reach steady state using implicit method with deltaT=0.3 sec is "
print m,"seconds"
print "Principles of Heat Transfer, 7th Ed. Frank Kreith et. al Chapter - 3 Example # 3.7 "
# Heat Transfer coefficient is given in W/m2-K as:
h = 200;
# Radius of cylinder in m is given as:
R0 = 0.05;
# Thermal conductivity in W/m-K is given as:
k = 20;
# Thermal diffusivityt in m2/sec is given as:
alpha = 10**(-5);
# Therefore the biot number is given as:
Bi = (h*R0)/k;
# Ambient water bath temperature in C is given as:
Tinfinity = 0;
# Initial temperature of centre line is given as:
T0 = 500;
# Final Temperature of centre line is given as:
Tr = 100;
# Therefore the value of (Tr-Tinfinity)/(T0-Tinfinity) is:
value = (Tr-Tinfinity)/(T0-Tinfinity); #Required value
# Umath.sing above value and biot number, from Figure 2.43 (a) on page 137, we have
# alpha*t/(R0*R0)=1.8
print "Therefore from chart solution, time taken in seconds shall be"
#Time taken in seconds
t = ((1.8*R0)*R0)/alpha
print t
# Proceeding to the numerical solution
#Because of symmetry we need to consider only one quarter of the circular cross section
#The vertical and horizontal radii are then adiabatic surfaces.
#We will have a total of nine types of control volume
#Each of the control volume energy balance equations can be solved
#The coefficient on Tfor control volume type 7 is:
#(dx*dx/(alpha*dt)) -2 -2*h*dx/5
#and for it to be positive
# value of t we use in the numerical solution must be smaller than this
# maximum value. The calculation is continued until the temperature for the control vol-ume nearest the cylinder axis is less than 100°C
print "And using numerical solution the time in seconds comes out to be"
#Time taken in seconds
tfinal = 431
print tfinal
print "which is about 4% less than the chart solution of 450 s."