1D_pfr_surfchem.ipynb (Source)
Data files used in this example: SiF4_NH3_mec.yaml1D 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 IDA solver to solve this system of differential algebraic equations (DAE).
Requires the scikits.odes.dae solver. The easiest way to install this package is from the conda-forge channel using conda.
In [1]:
import numpy as np
from scikits.odes import dae
import cantera as ct
%matplotlib inline
%config InlineBackend.figure_formats = ["svg"]
import matplotlib.pyplot as plt
plt.rcParams["figure.dpi"] = 120
print(f"Runnning Cantera version: {ct.__version__}")
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.yaml"
# 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])
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
gas_Si_N_interface.coverages = "HN_NH2(S): 1.0"
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:
\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 [4]:
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.get_net_production_rates(
"gas"
) # heterogeneous production rate of gas species
sdot_s = gas_Si_N_interface.get_net_production_rates("SI3N4")
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.
\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 [5]:
######## 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_g = gas_Si_N_interface.get_net_production_rates(
"gas"
) # heterogeneous production rate of gas species
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_g[:N] * W) / Ac
for i in range(gas.n_species):
b[1 + i] = (
wdot[i] * W[i] * Ac
+ sdot_g[i] * W[i] * perim
- gas.Y[i] * perim * np.sum(sdot_g[: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 [6]:
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)
print(solution)
Plot the results¶
In [7]:
# plot velocity of gas along the flow direction
f, ax = plt.subplots(3, 2, figsize=(9, 9))
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(r"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:
\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 [8]:
############################### initial conditions ##################################################################
# import the SiF4 + NH3 reaction mechanism
mech = "../data/SiF4_NH3_mec.yaml"
# 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
gas_Si_N_interface.coverages = "HN_NH2(S): 1.0"
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.get_net_production_rates(
"gas"
) # heterogeneous production rate of gas species
sdot_s = gas_Si_N_interface.get_net_production_rates("SI3N4")
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
sdot_Si = gas_Si_N_interface.get_net_production_rates(
"SiBulk"
) # bulk production rate
sdot_N = gas_Si_N_interface.get_net_production_rates("NBulk")
# 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 * (sdot_Si[0] * W_Si_b * h_Si + sdot_N[0] * W_N_b * h_N)
)
In [9]:
######## 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_g = gas_Si_N_interface.get_net_production_rates(
"gas"
) # heterogeneous production rate of gas species
sdot_s = gas_Si_N_interface.get_net_production_rates("SI3N4")
sdot_Si = gas_Si_N_interface.get_net_production_rates("SiBulk")
sdot_N = gas_Si_N_interface.get_net_production_rates("NBulk")
# 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_g * W) / Ac
for i in range(N):
b[1 + i] = (
wdot[i] * W[i] * Ac
+ sdot_g[i] * W[i] * perim
- gas.Y[i] * perim * np.sum(sdot_g * 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_g * W)
- perim * np.sum(sdot_Si[0] * W_Si * h_Si + sdot_N[0] * 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 [10]:
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 [11]:
f, ax = plt.subplots(4, 2, figsize=(9, 12))
# 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(r"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)