I am trying to solve the following heat equation for a rod using an explicit Runge-Kutta method in time:
$$ \frac{\partial T_{i}}{\partial t}=-u \frac{\partial T_i}{\partial z}+k \frac{\partial^{2} T_i}{\partial z^{2}}+h\left(T_{w}-T_{i}\right) \frac{P}{A} $$
For discretising in space, I am using
$$ \frac{\partial T}{\partial z}=\frac{T_{i-1}+T_{i+1}}{2 \Delta z} $$
and
$$ \frac{\partial^{2} T}{\partial t^{2}}=\frac{\left(T_{i+1}-2 T_{i}+T_{i-1}\right)}{\Delta z^{2}} $$
I have written the following code.
import matplotlib.pyplot as plt
#Set space conditions
nz = 51 #number of space grid points
L = 1.0 #length of rod
dz = L/(nz-1)
z = np.linspace(-1/2*dz, L + 1/2*dz, nz)
z[0] = 0
z[-1] = L #make last gridpoint align with end of rod
#Initial time conditions
t0 = 0
tf = 100
nt = 1000
dt = (tf-t0)/nt
t1 = np.linspace(t0, tf, nt)
#Initial temperature conditions
T0 = 20 * np.ones(nz) #temperature at time t0
T = np.zeros(nz)
T[:] = T0[:] #initialise temperature
dT_dt = np.zeros(nz)
#Define problem
#Parameters
u = 1.0 #velocity, 1.0 for simplicity
rho = 1.0 #density, 1.0 for simplicity
h = 1.0 #heat transfer coefficient, 1.0 for simplicity
Tw = 1.0 #temperature of wall, 1.0 for simplicity
L = 1.0
A = 1.0
P = 1.0
k = 1.0
cp = 1.0
def heattransfer(t, T):
for i in range (1, nz-1):
dT_dt[i] = - u*(T[i+1]-T[i-1])/(2*dz) + (k/rho/cp)*(T[i+1]-2*T[i]+T[i-1])/dz**2+ (h*P/A/rho/cp)*(Tw - T[i])
def RK(f, T0, t0, tf, dt):
for i in range (nt - 1):
time.append(t) ; T_RK.append(T)
k1 = f(t[i], T[i])
k2 = f(t[i]+ dt/2, T[i] + k1/2)
k3 = f(t[i]+ dt/2, T[i] + k2/2)
k4 = f(t[i]+ dt, T[i] + k3)
dT = (k1 + 2*k2 + 2*k3 + k4)/6
T[i+1] = T[i] + dT*dt
return t, T
f = lambda t, T: heattransfer(t, T)
t, T = RK(f, T0, t0, tf, dt)
The code raises the following error: Error raised in code
What have I done wrong?