# # 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()