1D_pfr_surfchem.ipynb (Source)

1D Plug Flow Reactor Model with Surface Chemistry

In this model, we will illustrate the derivation of the governing differential equations and algebraic constraints, calculation of the initial conditions of the variables and their spatial derivatives and use the scikits.odes.dae IDA solver to solve this system of differential algebraic equations (DAE).

In [1]:
from __future__ import print_function, division
import numpy as np
from scikits.odes import dae
import cantera as ct
import matplotlib.pyplot as plt
%matplotlib inline
print('Runnning Cantera version: ' + ct.__version__)
Runnning Cantera version: 2.4.0

Define gas species, bulk species, surface species and the interface

Here, we use a kinetic mechanism involving the chemical vapor deposition of silicon nitride (Si3N4) from SiF4 and NH3. 25 gas species, 6 surface species and 2 bulk species mechanism is applied by Richard S. Larson et al. 1996, SAND96-8211.

In [2]:
#import the SiF4 + NH3 reaction mechanism
mech = 'data/SiF4_NH3_mec.cti'
#import the models for gas and bulk
gas, bulk_Si, bulk_N = ct.import_phases(mech, ['gas', 'SiBulk', 'NBulk'])
#import the model for gas-Si-N interface
gas_Si_N_interface = ct.Interface(mech, 'SI3N4', [gas,bulk_Si,bulk_N])

Case 1: isothermal reactor

Define reactor conditions : temperature, pressure, fuel, and some important parameters

In [3]:
T0 = 1713  # Kelvin
p0 = 2 * ct.one_atm / 760.0  # Pa ~2 Torr
gas.TPX = T0, p0, "NH3:6, SiF4:1"
bulk_Si.TP = T0, p0
bulk_N.TP = T0, p0
gas_Si_N_interface.TP = T0, p0
D = 5.08e-2  # diameter of the tube [m]
Ac = np.pi * D**2/4  # cross section of the tube [m^2]
mu = 5.7e-5  # kg/(m-s) dynamic viscosity
perim = np.pi * D  # perimeter of the tube
# calculate the site fractions of surface species at the entrance of the tube at steady state
gas_Si_N_interface.advance_coverages(100.0)  # Here we assume after 100s, the system reaches the steady state
Zk_0 = gas_Si_N_interface.coverages
N = gas.n_species  # number of gas species
M = gas_Si_N_interface.n_species  # number of surface species

Define a residual function for IDA solver

For the isothermal tube with laminar flow, since the temperature of the flow and tube is constant, the energy conservation equation can be ignored. The governing equations include conservation of mass and species, momentum equation, equation of state, and the algebraic constraints that the net production rate of surface species by heterogeneous reactions are zero and that the sum of site fractions equals 1.

Here we define a residual function, an equation which should always evaluate to the zero vector, as the input of IDA solver, which listed as follows:

In [4]:
%%latex
\begin{align}
    R[0] &= u\frac{d\rho}{dz} + \rho\frac{du}{dz} - \frac{p'}{A_c}\sum^{K_g}\dot{s}_{k,g}W_{k,g} \\
    R[1] &= \rho u A_c\frac{dY_k}{dz} + Y_k p'\sum^{K_g}\dot{s}_{k,g}W_{k,g} - \dot{\omega_k}W_kA_c - \dot{s}_{k,g}W_{k,g} p' \\
    R[2] &= 2\rho u \frac{du}{dz} + u^2\frac{d\rho}{dz} + \frac{dP}{dz} + \frac{32u\mu}{D^2}\\
    R[3] &= P\bar{W} - \rho RT\\
    R[4] &= \dot{s}_{k,s} \\
    R[5] &= \sum_{phase}{Z_{k,s}} - 1
\end{align}
\begin{align} R[0] &= u\frac{d\rho}{dz} + \rho\frac{du}{dz} - \frac{p'}{A_c}\sum^{K_g}\dot{s}_{k,g}W_{k,g} \\ R[1] &= \rho u A_c\frac{dY_k}{dz} + Y_k p'\sum^{K_g}\dot{s}_{k,g}W_{k,g} - \dot{\omega_k}W_kA_c - \dot{s}_{k,g}W_{k,g} p' \\ R[2] &= 2\rho u \frac{du}{dz} + u^2\frac{d\rho}{dz} + \frac{dP}{dz} + \frac{32u\mu}{D^2}\\ R[3] &= P\bar{W} - \rho RT\\ R[4] &= \dot{s}_{k,s} \\ R[5] &= \sum_{phase}{Z_{k,s}} - 1 \end{align}

The detailed derivation of the DAE system can be found in my report.

In [5]:
def residual(z, vec, vecp, result):
    """ we create the residual equations for the problem
        vec = [u, rho, Yk, p, Zk]
        vecp = [dudz, drhodz, dYkdz, dpdz, dZkdz]
    """
    # temporary variables
    u = vec[0]  # velocity
    rho = vec[1]  # density
    Y = vec[2:2+N]  # vector of mass fractions of all gas species
    p = vec[2+N]  # pressure
    Z = vec[3+N:]  # vector of site fractions of all surface species   

    dudz = vecp[0]  # velocity spatial derivative
    drhodz = vecp[1]  # density spatial derivative
    dYdz = vecp[2:2+N]  # mass fraction spatial derivative
    dpdz = vecp[2+N]  # pressure spatial derivative

    # Use unnormalized mass fractions to avoid over-constraining the system
    gas.set_unnormalized_mass_fractions(Y)
    gas.TP = T0,p

    bulk_Si.TP = T0,p
    bulk_N.TP = T0,p

    # Use unnormalized site fractions (coverages) to avoid over-constraining the system
    gas_Si_N_interface.set_unnormalized_coverages(Z)
    gas_Si_N_interface.TP = T0,p

    # temporary variables (based on the given state)
    coverages = gas_Si_N_interface.coverages  # site fraction vector
    sdot_g = gas_Si_N_interface.net_production_rates[:N]  # heterogeneous production rate of gas species
    sdot_s = gas_Si_N_interface.net_production_rates[-M:]
    wdot_g = gas.net_production_rates  # homogeneous production rate of gas species
    W_g = gas.molecular_weights  # vector of molecular weight of gas species

    # mass continuity equation
    result[0] = u*drhodz + rho*dudz - perim*np.sum(sdot_g*W_g)/Ac

    # conservation of species
    for k in range(N):
        result[1+k] = (rho*u*Ac*dYdz[k] + Y[k]*perim*np.sum(sdot_g*W_g)
                       - wdot_g[k]*W_g[k]*Ac
                       - sdot_g[k]*W_g[k]*perim)
    # conservation of momentum
    result[1+N] = 2*rho*u*dudz + np.power(u,2)*drhodz + dpdz + 32*u*mu/D**2 

    # equation of state
    result[2+N] = gas.density - rho

    # algebraic constraints
    for j in range(M):
        result[3+N+j] = sdot_s[j]

    # replace the constraint with the condition sum(Zk) = 1 for the largest site fraction species
    index = np.argmax(coverages)
    result[3+N+index] = np.sum(coverages) - 1

Determine the initial values of the spatial derivatives of the unknowns which need to be used as the initial conditions for the IDA solver

The following linear equation system has been solved by np.linalg.solve, a linear solver, to calculate the initial values of the spatial derivatives of the unknowns.

In [6]:
%%latex
\begin{align}
    u_0\rho_0' + \rho_0 u_0' - \frac{p'}{A_c}\sum^{K_g}\dot{s}_{k,g}W_{k,g} &= 0\\
    \rho_0 u_0 A_c Y_{k,0}' + Y_{k,0} p'\sum^{K_g}\dot{s}_{k,g}W_{k,g} - \dot{\omega_k}W_kA_c - \dot{s}_{k,g}W_{k,g} p' &=0 \\
    2\rho_0 u_0 u_0' + u_0^2\rho_0' + P_0' + \frac{32u_0 \mu}{D^2} &=0\\
    -RT\rho_0' + \bar{W_0}P_0' - P_0\frac{\sum^{K_g}Y_{k}'/W_{k,g}}{(\sum^{K_g}Y_{k}/W_{k,g})^2} &= 0
\end{align}
\begin{align} u_0\rho_0' + \rho_0 u_0' - \frac{p'}{A_c}\sum^{K_g}\dot{s}_{k,g}W_{k,g} &= 0\\ \rho_0 u_0 A_c Y_{k,0}' + Y_{k,0} p'\sum^{K_g}\dot{s}_{k,g}W_{k,g} - \dot{\omega_k}W_kA_c - \dot{s}_{k,g}W_{k,g} p' &=0 \\ 2\rho_0 u_0 u_0' + u_0^2\rho_0' + P_0' + \frac{32u_0 \mu}{D^2} &=0\\ -RT\rho_0' + \bar{W_0}P_0' - P_0\frac{\sum^{K_g}Y_{k}'/W_{k,g}}{(\sum^{K_g}Y_{k}/W_{k,g})^2} &= 0 \end{align}

We assume the derivatives of the site fractions are equal to zero, although it is trivial for the IDA solver.

In [7]:
######## Solve linear system for the initial vecp ###########
"""
   a = coefficient of [u', rho', Yk', P']
   b = RHS constant of each conservation equations
"""
rho0 = gas.density  # initial density of the flow
u0 = 11.53  # m/s initial velocity of the flow
W = gas.molecular_weights
W_avg = gas.mean_molecular_weight
sdot = gas_Si_N_interface.net_production_rates  # heterogeneous molar production rate
wdot = gas.net_production_rates  # homogeneous molar production rate
################### a #########################
a = np.zeros((3+N,3+N))
a[0,:] = np.hstack((rho0, u0, np.zeros(1+N)))
for i in range(N):
    a[1+i,2+i] = rho0*u0*Ac
a[1+N,:] = np.hstack((2*rho0*u0, u0**2, np.zeros(N), 1))
coef = np.zeros(N)
for j in range(N):
    coef[j] = gas.P/W[j]/np.power(np.sum(gas.Y/W),2)
a[2+N,:] = np.hstack((0, ct.gas_constant*T0, coef, -W_avg))
################### b ###########################
b = np.zeros(3+gas.n_species)
b[0] = perim*np.sum(sdot[:N]*W)/Ac
for i in range(gas.n_species):
    b[1+i] = (wdot[i]*W[i]*Ac
              + sdot[i]*W[i]*perim
              - gas.Y[i]*perim*np.sum(sdot[:N]*W))
b[1+gas.n_species] = -32*u0*mu/D**2
b[2+gas.n_species] = 0
part_vecp0 = np.linalg.solve(a,b)

vecp0 = np.hstack((part_vecp0, np.zeros(M)))
vec0 = np.hstack((11.53, gas.density, gas.Y, gas.P, Zk_0))

Run the IDA solver to calculate the unknowns varying in the flow direction

In [8]:
solver = dae(
    'ida',
    residual, 
    first_step_size=1e-16,
    atol=1e-8,  # absolute tolerance for solution
    rtol=1e-8,  # relative tolerance for solution
    # If the given problem is of type DAE, some items of the residual vector
    # returned by the 'resfn' have to be treated as algebraic equations, and
    # algebraic variables must be defined. These algebraic variables are
    # denoted by the position (index) in the state vector y. All these
    # indexes have to be specified in the 'algebraic_vars_idx' array.
    algebraic_vars_idx=[np.arange(3+N,3+N+M,1)], 
    max_steps=5000,
    old_api=False  # Forces use of new api (namedtuple)
)

times = np.arange(0,0.7,0.01)
solution = solver.solve(times, vec0, vecp0)

Plot the results

In [9]:
# plot velocity of gas along the flow direction
f, ax = plt.subplots(3,2, figsize=(9,9), dpi=96)
ax[0,0].plot(times, solution.values.y[:,0], color='C0')
ax[0,0].set_xlabel('Distance (m)')
ax[0,0].set_ylabel('Velocity (m/s)')

# plot gas density along the flow direction
ax[0,1].plot(times, solution.values.y[:,1], color='C1')
ax[0,1].set_xlabel('Distance (m)')
ax[0,1].set_ylabel('Density ($\mathregular{kg/m^3}$)')
ax[0,1].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation

# plot major and minor gas species separately
minor_idx = []
major_idx = []
for i,name in enumerate(gas.species_names): 
    mean = np.mean(solution.values.y[:,2+i])
    if mean <= 0.01:
        minor_idx.append(i) 
    else:
        major_idx.append(i)

# plot minor species
for i in minor_idx:
    style = '-' if i < 10 else '--'    
    ax[1,0].plot(times, solution.values.y[:,2+i], label=gas.species_names[i], linestyle=style)
ax[1,0].legend(fontsize=7, loc='upper right')
ax[1,0].set_xlabel('Distance (m)')
ax[1,0].set_ylabel('Mass Fraction')
ax[1,0].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation

# plot major species
for j in major_idx:
    ax[1,1].plot(times,solution.values.y[:,2+j], label=gas.species_names[j])
ax[1,1].legend(loc='best')
ax[1,1].set_xlabel('Distance (m)')
ax[1,1].set_ylabel('Mass Fraction')

# plot the pressure of the gas along the flow direction
ax[2,0].plot(times, solution.values.y[:,2+N], color='C2')
ax[2,0].set_xlabel('Distance (m)')
ax[2,0].set_ylabel('Pressure (Pa)')

# plot the site fraction of the surface species along the flow direction 
for i,name in enumerate(gas_Si_N_interface.species_names):
    ax[2,1].plot(times, solution.values.y[:,3+N+i], label=name)
ax[2,1].legend()
ax[2,1].set_xlabel('Distance (m)')
ax[2,1].set_ylabel('Site Fraction')
f.tight_layout(pad=0.5)

Case 2: Adiabatic reactor

Since the application of isothermal reactor is not prevalent, to improve the model for real use, the adiabatic reator is considered. Here, the energy balance equation is also considered.

The heat flow rate into the system has two components. One is due to the heat flux $q_e$ from the surroundings to the outer tube wall (whose surface area per unit length is $a_e$) and accumulation of enthalpy in the bulk solid. The other is due to $q_i$, the heat flux to the gas from the inner tube wall, and accumulation of enthalpy in the surface species. The expression of energy balance equation for this problem is as follows:

In [10]:
%%latex
\begin{align}
    \rho u A_c c_p \frac{dT}{dz} +A_c \sum_{K_g}\dot{\omega}_k W_k h_k + p'\sum_{K_g}h_k\dot{s}_k W_k &= a_eq_e - p'\sum^{K_b}_{bulk}\dot{\omega}_kh_k\\&=p'q_i + p'\sum^{K_g}_{gas}\dot{s_k}W_kh_k
\end{align}
\begin{align} \rho u A_c c_p \frac{dT}{dz} +A_c \sum_{K_g}\dot{\omega}_k W_k h_k + p'\sum_{K_g}h_k\dot{s}_k W_k &= a_eq_e - p'\sum^{K_b}_{bulk}\dot{\omega}_kh_k\\&=p'q_i + p'\sum^{K_g}_{gas}\dot{s_k}W_kh_k \end{align}

Since the adiabatic reactor is considered, $q_e = 0$. Similar to the procedure for the isothermal reactor model, add the energy equation into the residual function and calculate the initial value of the spatial derivative of the temperature.

In [11]:
############################### initial conditions ##################################################################
# import the SiF4 + NH3 reaction mechanism
mech = 'data/SiF4_NH3_mec.cti'
# import the models for gas and bulk
gas, bulk_Si, bulk_N = ct.import_phases(mech,['gas','SiBulk','NBulk'])

# import the model for gas-Si-N interface
gas_Si_N_interface = ct.Interface(mech, 'SI3N4',[gas,bulk_Si,bulk_N])
T0 = 1713  # K
p0 = 2 * ct.one_atm / 760.0  # Pa ~2Torr
gas.TPX = T0, p0, "NH3:6, SiF4:1"
bulk_Si.TP = T0, p0
bulk_N.TP = T0, p0
gas_Si_N_interface.TP = T0, p0
D = 5.08e-2  # diameter of the tube [m]
Ac = np.pi * D**2/4  # cross section of the tube [m]
mu = 5.7e-5  # kg/(m-s) dynamic viscosity
perim = np.pi * D  # perimeter of the tube
# calculate the site fractions of surface species at the entrance of the tube at steady state
gas_Si_N_interface.advance_coverages(100.0)
Zk_0 = gas_Si_N_interface.coverages
######################################## IDA solver ###################################################################
def residual(z, vec, vecp, result):
    """ we create the residual equations for the problem
        vec = [u, rho, Yk, p, Zk, T]
        vecp = [dudz, drhodz, dYkdz, dpdz, dZkdz, dTdz]
    """
    # temporary variables
    u = vec[0]  # velocity
    rho = vec[1]  # density
    Y = vec[2:2+N]  # vector of mass fractions of all gas species
    p = vec[2+N]  # pressure
    Z = vec[3+N:-1]  # vector of site fractions of all surface species
    T = vec[-1]  # temperature
    
    dudz = vecp[0]  # velocity spatial derivative
    drhodz = vecp[1]  # density spatial derivative
    dYdz = vecp[2:2+N]  # mass fraction spatial derivative
    dpdz = vecp[2+N]  # pressure spatial derivative
    dTdz = vecp[-1]  # temperature spatial derivative
    
    h = gas.enthalpy_mass  # enthalpy of gas species per mass
    h_Si = bulk_Si.enthalpy_mass  # enthalpy of Si per mass
    h_N = bulk_N.enthalpy_mass  # enthalpy of N per mass
    
    # initial conditions
    gas.set_unnormalized_mass_fractions(Y)
    gas.TP = T, p
    
    bulk_Si.TP = T, p
    bulk_N.TP = T, p
    gas_Si_N_interface.set_unnormalized_coverages(Z)
    gas_Si_N_interface.TP = T, p
    
    # temporary variables (based on the current system state)
    coverages = gas_Si_N_interface.coverages  # site fraction vector
    sdot_g = gas_Si_N_interface.net_production_rates[:N]  # heterogeneous production rate of gas species
    sdot_s = gas_Si_N_interface.net_production_rates[-M:]  # molar production rate of surface speceis
    wdot_g = gas.net_production_rates  # homogeneous production rate of gas species
    W_g = gas.molecular_weights  # vector of molecular weight of gas species
    W_Si_b = bulk_Si.molecular_weights
    W_N_b = bulk_N.molecular_weights
    bdot = gas_Si_N_interface.net_production_rates[gas.n_species:gas.n_species+2]  # bulk production rate
   
    # mass continuity equation
    result[0] = u*drhodz+rho*dudz-perim*np.sum(sdot_g*W_g)/Ac
    # conservation of species
    for k in range(gas.n_species):
        result[1+k] = (rho*u*Ac*dYdz[k] + Y[k]*perim*np.sum(sdot_g*W_g)
                       - wdot_g[k]*W_g[k]*Ac
                       - sdot_g[k]*W_g[k]*perim)
    # conservation of momentum
    result[1+gas.n_species] = 2*rho*u*dudz + np.power(u,2)*drhodz + dpdz + 32*u*mu/D**2

    # equation of state
    result[2+gas.n_species] = gas.density - rho

    # algebraic constraints
    for j in range(M):
        result[3+N+j] = sdot_s[j]
    
    # replace the constraints with the condition sum(Zk) = 1 for the largest site fraction species
    index = np.argmax(coverages)
    result[3+N+index] = np.sum(coverages) - 1
    
    # energy equation
    result[3+N+M] = (rho*u*Ac*gas.cp*dTdz
        + Ac*np.sum(wdot_g*W_g*h)
        + perim*np.sum(h*sdot_g*W_g)
        + perim*(bdot[0]*W_Si_b*h_Si + bdot[1]*W_N_b*h_N))
In [12]:
######## Solve linear system for the initial values of vecp ###########
"""
   a = coefficient of [u', rho', Yk', P',T]
   b = RHS constant of each conservation equations
"""
rho0 = gas.density  # initial density of the flow
u0 = 11.53  # m/s initial velocity of the flow
W = gas.molecular_weights
W_avg = gas.mean_molecular_weight
sdot = gas_Si_N_interface.net_production_rates  # heterogeneous molar production rate
wdot = gas.net_production_rates  # homogeneours molar production rate
h = gas.enthalpy_mass
h_Si = bulk_Si.enthalpy_mass
h_N = bulk_N.enthalpy_mass
W_Si = bulk_Si.molecular_weights
W_N = bulk_N.molecular_weights
################### a #########################
a = np.zeros((4+N,4+N))
a[0,:] = np.hstack((rho0, u0, np.zeros(2+N)))
for i in range(N):
    a[1+i,2+i] = rho0*u0*Ac
a[1+N,:] = np.hstack((2*rho0*u0, u0**2, np.zeros(N), 1,0))
coef = np.zeros(N)
for j in range(N):
    coef[j] = gas.P/W[j]/np.power(np.sum(gas.Y/W),2)
a[2+N,:] = np.hstack((0, ct.gas_constant*T0, coef, -W_avg, 0))
a[3+N,:] = np.hstack((np.zeros(3+N), rho0*u0*Ac*gas.cp))
################### b ###########################
b = np.zeros(4+gas.n_species)
b[0] = perim*np.sum(sdot[:N]*W)/Ac
for i in range(N):
    b[1+i] = (wdot[i]*W[i]*Ac
              + sdot[i]*W[i]*perim
              - gas.Y[i]*perim*np.sum(sdot[:N]*W))
b[1+gas.n_species] = -32*u0*mu/D**2
b[2+gas.n_species] = 0
b[3+gas.n_species] = (- Ac*np.sum(wdot*W*h)
                      - perim*np.sum(h*sdot[:N]*W)
                      - perim*np.sum(sdot[N]*W_Si*h_Si + sdot[N+1]*W_N*h_N))
part_vecp0 = np.linalg.solve(a,b)

vecp0 = np.hstack((part_vecp0[:-1], np.zeros(M), part_vecp0[-1]))
vec0 = np.hstack((11.53, gas.density, gas.Y, gas.P, Zk_0, T0))
In [13]:
solver = dae(
    'ida',
    residual, 
    atol=1e-8,  # absolute tolerance for solution
    rtol=1e-8,  # relative tolerance for solution
    algebraic_vars_idx=[np.arange(3+N,3+N+M,1)], 
    max_steps=5000,
    one_step_compute=True,
    old_api=False
)

time = []
solution = []
state = solver.init_step(0.0, vec0, vecp0)
while state.values.t < 0.7:
    time.append(state.values.t)
    solution.append(state.values.y)
    state = solver.step(0.7)

time = np.array(time)
solution = np.array(solution)
In [14]:
f, ax = plt.subplots(4,2, figsize=(9,12), dpi=96)

# plot gas velocity along the flow direction
ax[0,0].plot(time, solution[:,0], color='C0')
ax[0,0].set_xlabel('Distance (m)')
ax[0,0].set_ylabel('Velocity (m/s)')

# plot gas density along the flow direction
ax[0,1].plot(time, solution[:,1], color='C1')
ax[0,1].set_xlabel('Distance (m)')
ax[0,1].set_ylabel('Density ($\mathregular{kg/m^3}$)')
ax[0,1].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation

# plot major and minor gas species separately
minor_idx = []
major_idx = []
for i,name in enumerate(gas.species_names): 
    mean = np.mean(solution[:,2+i])
    if mean <= 0.01:
        minor_idx.append(i) 
    else:
        major_idx.append(i)

# plot minor gas species along the flow direction
for i in minor_idx:
    style = '-' if i < 10 else '--'
    ax[1,0].plot(time, solution[:,2+i], label=gas.species_names[i], linestyle=style)
ax[1,0].legend(fontsize=7.5, loc='best')
ax[1,0].set_xlabel('Distance (m)')
ax[1,0].set_ylabel('Mass Fraction')
ax[1,0].ticklabel_format(axis='y', style='sci', scilimits=(-2,2))  # scientific notation

# plot major gas species along the flow direction
for j in major_idx:
    ax[1,1].plot(time, solution[:,2+j], label=gas.species_names[j])
ax[1,1].legend(fontsize=8, loc='best')
ax[1,1].set_xlabel('Distance (m)')
ax[1,1].set_ylabel('Mass Fraction')

# plot the pressure of the gas along the flow direction
ax[2,0].plot(time, solution[:,2+N], color='C2')
ax[2,0].set_xlabel('Distance (m)')
ax[2,0].set_ylabel('Pressure (Pa)')

# plot the site fraction of the surface species along the flow direction
for i,name in enumerate(gas_Si_N_interface.species_names):
    ax[2,1].plot(time, solution[:,3+N+i], label=name)
ax[2,1].legend(fontsize=8)
ax[2,1].set_xlabel('Distance (m)')
ax[2,1].set_ylabel('Site Fraction')

# plot the temperature profile along the flow direction
ax[3,0].plot(time, solution[:,-1], color='C3')
ax[3,0].set_xlabel('Distance (m)')
ax[3,0].set_ylabel('Temperature (K)')
f.tight_layout(pad=0.5)