The heat equation in a 1D Cartesian system (Dirichlet BC)

We compare the FVM solution by PyFVTool to the analytic solutions, for two different types of boundary conditions.

import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt
import pyfvtool as pf

# Spoiler alert!

def T_analytic_dirichlet(x,t, alfa, T0, Ts):
    return (T0-Ts)*erf(x/np.sqrt(4*alfa*t))+Ts

def T_analytic_neumann(x,t, alfa, T0, k, qs):
    return T0 + qs/k*np.sqrt(4*alfa*t/np.pi)*np.exp(-x**2/(4*alfa*t))\
              - qs/k*x*(1-erf(x/np.sqrt(4*alfa*t)))

# Parameters
L = 1.0 # [m] domain length
k = 20.0 # 0.6 for water, 0.025 for air W/m/K
rho = 8000.0 # kg/m^3
c = 500.0 # J/kg/K (4200 for water, 1000 for air)
alfa = k/(rho*c) # heat diffusion
T0 = 300.0 # [K]
Ts = 350.0 # [K]
qs = 1000 # [W/m^2]
t_sim = L**2/(20*alfa) # [s]
time_steps = 50
dt = t_sim/time_steps #
Nx = 50 # number of cells

m = pf.Grid1D(Nx, L)

# Cell variable with initial condition
T = pf.CellVariable(m, T0)

# Boundary condition
left_bc = "Dirichlet"

BC = pf.BoundaryConditions(m)
if left_bc == "Dirichlet":
    T.BCs.left.a[:] = 0.0
    T.BCs.left.b[:] = 1.0
    T.BCs.left.c[:] = Ts
    T_analytic = lambda x,t: T_analytic_dirichlet(x, t, alfa, T0, Ts)
else:
    T.BCs.left.a[:] = k
    T.BCs.left.b[:] = 0.0
    T.BCs.left.c[:] = -qs
    T_analytic = lambda x,t: T_analytic_neumann(x, t, alfa, T0, k, qs)

# physical parameters
alfa_cell = pf.CellVariable(m, alfa)
alfa_face = pf.harmonicMean(alfa_cell)

t=0
while t<t_sim:
    pf.solvePDE(T, [ pf.transientTerm(T, dt, 1.0),
                    -pf.diffusionTerm(alfa_face)])
    t +=dt

x = m.facecenters.x
T_face = pf.linearMean(T)
T_num = T_face.xvalue
T_an = T_analytic(x, t_sim)

plt.figure(1)
plt.clf()
plt.title('Left BC: '+left_bc)
plt.plot(x, T_an, x, T_num, 'o')
plt.legend({'Analytic', 'PyFVTool FVM'})
plt.xlabel('x [m]')
plt.ylabel('T [K]')
plt.show()

er = np.sum(np.abs(T_num-T_an)/T_an)/Nx
print(er)

0.00021832771380885904