Runge-Kutta method for heat equation

422 Views Asked by At

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?