Custom rate expression in a 1-D plug flow reactor

[Jacob Miller]

Hello Cantera community,

I am currently trying to simulate reactions in a steady-state 1-D gas-phase plug flow reactor (only incorporating convection, not diffusion) using Cantera in Python. My kinetics are derived from empirical fitting of experimental data, so I would like to implement it as a custom rate expression. The actual chemistry I am modeling is happening on the surface of a heterogeneous catalyst, but I am perfectly fine modeling it as homogenous and incorporating terms for porosity. I have set up a script and a yaml file (both below), but I am having two problems:

1) For some reason, my conversion is increasing as the number of volume/distance increments I include in the for loop in the code increases. I am not sure why this is happening, but it may have something to do with the sim.step() command in that for loop.

2) My empirical rate expression (reactions[0]) is a function of species partial pressures. As of now, I don't see a way to easily incorporate species partial pressures into the rate expression. Currently, the code has a dummy rate expression dependent on only temperature because I'm trying to troubleshoot my first problem, but the rate expression I want to use is commented out right below it. Is it possible to compute an empirical rate expression that has contributions from species pressures?

Thanks in advance for the help!
Jake Miller

yaml file:

units:
  length: cm
  quantity: mol
  activation-energy: J/mol
  pressure: Pa

phases:
  - name: gas
    thermo: ideal-gas
    elements: [O, H, C, N]
    species: [BA, H2O, C7H14O, CO2, N2]
    kinetics: gas
    reactions: all
    #reactions: {reversible: no}
    #transport: mixture-averaged

species:
  - name: BA
    composition: {C: 4, H: 8, O: 2}
    thermo:
      model: Shomate
      temperature-ranges: [623.00, 673.00]
      data:
        - [187,0.19,0,0,0,-426900,471.0]

  - name: H2O
    composition: {H: 2, O: 1}
    thermo:
      model: Shomate
      temperature-ranges: [500.00, 1700.00, 6000.00]
      data:
        - [30.092, 6.832514, 6.793435, -2.53448, 0.082139, -250.881, 223.3967]
        - [41.96426, 8.622053, -1.49978, 0.098119, -11.15764, -272.1797, 219.7809]
     
        # - [30.092, 6.832514, 6.793435, -2.53448, 0.082139, -250.881, 223.3967, -241.8264]
        # - [41.96426, 8.622053, -1.49978, 0.098119, -11.15764, -272.1797, 219.7809, -241.8264]

  - name: C7H14O
    composition: {C: 7, H: 14, O: 1}
    thermo:
      model: Shomate
      temperature-ranges: [623.00, 673.00]
      data:
        - [289.2,0.266,0,0,0,-221720,616.7]
        # 623.00, 673.00,
        #   289.2, 0.266,     # Cp(T) = 289.2 + 0.266 * (T - 623) / J/K/mol
        #   -221720, 29.6,    # H(T) = -221720 + 29.6 * (T - 623) / J/mol
        #   616.7, 0.456]     # S(T) = 616.7 + 0.456 * (T - 623) / J/K/mol

  - name: CO2
    composition: {C: 1, O: 2}
    thermo:
      model: Shomate
      temperature-ranges: [100.00, 500.00, 2000.00]
      data:
        # - [24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, 228.2431, -393.5224]
        # - [41.96426, 8.622053, -1.49978, 0.098119, -11.15764, -272.1797, 219.7809, -241.8264]
        - [24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, 228.2431]
        - [41.96426, 8.622053, -1.49978, 0.098119, -11.15764, -272.1797, 219.7809]


  - name: N2
    composition: {N: 2}
    thermo:
      model: Shomate
      temperature-ranges: [298.00, 1200.00, 6000.00]
      data:
        - [28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 226.4168] #0 was H
        - [19.50583, 19.88705, -8.598535, 1.369784, 0.527601, -4.935202, 212.39]

reactions:
  - equation: 2 BA => H2O + C7H14O + CO2
    rate-constant: {A: 1.0, b: 0.0, Ea: 0.0}

Python script:

import csv
import os
import cantera as ct
import numpy as np

# Set CANTERA_DATA environment variable to the current directory
os.environ['CANTERA_DATA'] = os.getcwd()

# Unit conversion factors to SI

#######################################################################
# Input Parameters
#######################################################################

tc = 380  # Temperature in Celsius
length = 10000  # Catalyst bed length /m
area = 0.25 * np.pi  # Catalyst bed area /m^2, assuming a 1m diameter reactor
rhobulk = 5.68  # g/cm^3
ebed = 0.41568
theta0 = 0.59625
cat_mass_per_vol = 5.68 * ebed * theta0 * 1000  # Catalyst particle mass per unit volume, kg/m^3
m_ins = np.divide([1, 9, 0, 0, 0.1], 60)  # in kg/s
MWs = np.divide([88.11, 18.01, 44.01, 114.19, 28.02], 1000)  # in kg/mol
F_ins = np.divide(m_ins, MWs)  # in mol/s
ndot = np.sum(F_ins)  # in mol/s
P_in = 101325  # Pa
R = 8.314  # J/K/mol
velocity = ndot * R * (tc + 273) / (ebed * area)  # gas velocity, m/s
porosity = ebed  # Catalyst bed porosity
total_cat_mass = cat_mass_per_vol * area * length * (1 - porosity)  # total catalyst mass in kg

# Input file containing the surface reaction mechanism
yaml_file = 'but_ket_thermo.yaml'
output_filename = 'butyric_output.csv'

#####################################################################

t = tc + 273.15  # convert to Kelvin

# Import the model and set the initial conditions
gas = ct.Solution(yaml_file, 'gas')

reactions = gas.reactions()

# Define reaction
reactions[0] = ct.Reaction(
    equation= '2 BA => H2O + C7H14O + CO2',
    rate=lambda T: 9630000.0 * T**2.0 )
#     rate-constant: {A: 1.0, b: 0.0, Ea: 0.0}
#     rate-coefficients:
#       k: 2.60352e-4/60  # in mol/kgcat/s
#       K_BA: 8.682e-3  # in Pa^-1
#       K_W: 6.60e-4  # in Pa^-1
#     rate-expression: |
#       k * P_BA^2 / ( (1 + K_BA * P_BA)^2 * (1 + K_W * P_W) )
# )
   
gas = ct.Solution(thermo='ideal-gas', kinetics='gas',
                   species=gas.species(), reactions=reactions)

# Calculate mole fractions from molar flow rates
X = F_ins / ndot
species_names = ['BA', 'H2O', 'C7H14O', 'CO2', 'N2']
mole_fractions = dict(zip(species_names, X))

gas.TPX = t, P_in, mole_fractions
print(gas.net_production_rates)
mass_flow_rate = sum(m_ins)

# Create a new reactor
r = ct.FlowReactor(gas)
r.area = area
r.mass_flow_rate = mass_flow_rate
r.energy_enabled = True

sim = ct.ReactorNet([r])

output_data = []
n_steps = 1000
distance = np.linspace(0, length, n_steps)

print('    distance       X_BA        X_H2O        X_C7H14O   X_CO2')
for dist in distance:
    r_volume = area * dist  # Volume of the reactor up to the current distance
    catalyst_mass = cat_mass_per_vol * r_volume * (1 - porosity)  # Catalyst mass in the current reactor volume

    r.volume = r_volume  # Set the reactor volume
    sim.step()

    # Capture output data
    output_data.append([dist] + list(r.thermo['BA', 'H2O', 'C7H14O', 'CO2'].X))
   
    # Print the results at intervals
    if len(output_data) % 10 == 0:
        print(f'{dist * 1e3:10.3f}  {r.thermo["BA"].X[0]:10.5f}  {r.thermo["H2O"].X[0]:10.5f}  '
              f'{r.thermo["C7H14O"].X[0]:10.5f}  {r.thermo["CO2"].X[0]:10.5f}')

# Write the output to a CSV file
with open(output_filename, 'w', newline='') as outfile:
    writer = csv.writer(outfile)
    writer.writerow(['Distance (mm)', 'X_BA', 'X_H2O', 'X_C7H14O', 'X_CO2'])
    writer.writerows(output_data)

print(f"Results saved to '{output_filename}'")

[Ray Speth]

Hi Jake,

For the first issue, I think there’s some confusion in your integration loop. Following surf_pfr.py, when using sim.step, the current distance is an output of the integrator, not something you specify. So your integration loop should look more like

while sim.distance < length: sim.step() dist = sim.distance output_data.append([dist] + list(r.thermo['BA', 'H2O', 'C7H14O', 'CO2'].X))

Regarding your second question, this is exactly what the ExtensibleRate class was introduced for. Building on the example given in custom_reactions.py, you will want to modify the update method to store the partial pressures and any other thermodynamic properties needed to calculate the rate. This method should return True if anything affecting the rate calculation has changed, and False otherwise. This state information will then be available in the data object passed to the eval method.

Regards,
Ray

​