# # Author: Larry Engelhardt (July, 2016)
# Material properties for a few different materials are listed below
#
# 1. Thermal conductivity: `k_t`
# 2. Specific heat capacity: `c`
# 3. Density: `rho`
#
# Stainless Steel
k_t = 16 # W / (m C)
c = 500 # J / (kg C)
rho = 8000 # kg / m^3
# Bakelite
k_t = 0.2 # W / (m C)
c = 920 # J / (kg C)
rho = 1300 # kg / m^3
# Oak
k_t = 0.17 # W / (m C)
c = 2000 # J / (kg C)
rho = 700 # kg / m^3
# Fiberglass
k_t = 0.04 # W / (m C)
c = 700 # J / (kg C)
rho = 2000 # kg / m^3
# BEGINNING THE PROGRAM:
from pylab import *
# Discretizing space
L = 0.15 # Length of handle in meters
Nx = 76 # Number of positions in space
x = linspace(0, L, Nx) # Array of positions
dx = L / (Nx - 1) # Values of delta x (meters)
# Discretizing time
tMax = 600 # Maximum time in seconds
Nt = 1801 # Number of time steps
t = linspace(0, tMax, Nt) # Array of time values
dt = tMax / (Nt - 1) # Size of time step (seconds)
T0 = 72 # Initial temperature (Farenheit)
dT_stove = 300 # Temperature increase of the stove (F)
# Stainless Steel's material properties
k_t = 16 # W / (m C)
c = 500 # J / (kg C)
rho = 8000 # kg / m^3
r = k_t * dt / (c * rho * dx**2) ## THE PARAMETER! (dimensionless) ##
print('r:', r) # Check to make sure that r < 0.5
# Preparing 2D array of temperatures
T = zeros((Nt, Nx)) # 2D array for temperature. NOTE THE ORDER: [time, space]
T[0, :] = T0 # Initial values (time t = 0)
tau = 60 # Used in tanh function, units of seconds
T[:, 0] = T0 + dT_stove * tanh(t/tau) # Edge of rod (at x=0) in deg. F
T = (T - 32) / 1.8 # Convert all temperatures to Celsius
# Verifying that $T(t)$ for the pan ($x=0$) looks correct:
Tpan = T0 + dT_stove * tanh(t/tau)
plot(t, Tpan, lw=2)
grid(True)
xlabel('Time (sec)')
ylabel('Pan Temperature (F)')
grid(True); show()
# Computing the temperature of the handle for all `x` and `t`:
for j in range(1, Nt): # Loop thru TIME
for i in range(1, Nx-1): # Loop thru SPACE
T[j,i] = r*T[j-1, i+1] + (1 - 2*r)*T[j-1,i] + r*T[j-1, i-1]
T[j, -1] = (1-r)*T[j-1,i] + r*T[j-1,i-1] # End of rod (x=L)
T = 1.8*T + 32 # Convert back to F from C
x = x * 100 # Convert from meters to cm
# Plotting $T(x)$ for a few moments in time
jValues = [30, 180, 540, -1] # A few values of the time index (-1 for last element)
for j in jValues:
plot(x, T[j,:], linewidth=2, label='Time (seconds): '+str(t[j]))
legend()
title('Stainless Steel')
xlim([0, 15])
xlabel('Position (cm)')
ylabel('Temperature ($^{\circ}$F)')
grid(True);
savefig('StainlessSteel_T_vs_x.png')
show()
# Creating a contour plot of $T(x, t)$
colorvals = linspace(70, 160, 91) # colors for contour plot (min, max, num)
contourf(x, t, T, colorvals) # (x-axis, y-axis, z(x,y))
title('Stainless Steel')
xlabel('x (cm)')
ylabel('Time (sec)')
colorbar()
grid(True)
savefig('StainlessSteel_contour.png')
show()
# Repeating using Bakelite
# Discretizing space
L = 0.15 # Length of handle in meters
Nx = 76 # Number of positions in space
x = linspace(0, L, Nx) # Array of positions
dx = L / (Nx - 1) # Values of delta x (meters)
# Discretizing time
tMax = 600 # Maximum time in seconds
Nt = 1801 # Number of time steps
t = linspace(0, tMax, Nt) # Array of time values
dt = tMax / (Nt - 1) # Size of time step (seconds)
T0 = 72 # Initial temperature (Farenheit)
dT_stove = 300 # Temperature increase of the stove (F)
# Bakelite's material properties
k_t = 0.2 # W / (m C)
c = 920 # J / (kg C)
rho = 1300 # kg / m^3
r = k_t * dt / (c * rho * dx**2) ## THE PARAMETER! (dimensionless) ##
print('r:', r) # Check to make sure that r < 0.5
# Preparing 2D array of temperatures
T = zeros((Nt, Nx)) # 2D array for temperature. NOTE THE ORDER: [time, space]
T[0, :] = T0 # Initial values (time t = 0)
tau = 60 # Used in tanh function, units of seconds
T[:, 0] = T0 + dT_stove * tanh(t/tau) # Edge of rod (at x=0) in deg. F
T = (T - 32) / 1.8 # Convert all temperatures to Celsius
for j in range(1, Nt): # Loop thru TIME
for i in range(1, Nx-1): # Loop thru SPACE
T[j,i] = r*T[j-1, i+1] + (1 - 2*r)*T[j-1,i] + r*T[j-1, i-1]
T[j, -1] = (1-r)*T[j-1,i] + r*T[j-1,i-1] # End of rod (x=L)
T = 1.8*T + 32 # Convert back to F from C
x = x * 100 # Convert from meters to cm
jValues = [30, 180, 540, -1] # A few values of the time index (-1 for last element)
for j in jValues:
plot(x, T[j,:], linewidth=2, label='Time (seconds): '+str(t[j]))
legend()
title('Bakelite')
xlim([0, 15])
xlabel('Position (cm)')
ylabel('Temperature ($^{\circ}$F)')
grid(True);
savefig('StainlessSteel_T_vs_x.png')
show()
colorvals = linspace(70, 160, 91) # colors for contour plot (min, max, num)
contourf(x, t, T, colorvals) # (x-axis, y-axis, z(x,y))
title('Bakelite')
xlabel('x (cm)')
ylabel('Time (sec)')
colorbar()
grid(True)
savefig('Bakelite_contour.png')
show()