0
$\begingroup$

I have a system of chemical reactions where the rate equations represent a batch reactor model. The model is a system of ODEs which is solved with the SciPy solve_ivp function. My example (see below) currently contains 5 reactions but I need to add many more reactions from the overall scheme. With my current approach it has been tedious to keep track of the chemical species, mass fractions, rate parameters, etc.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Parameters
# -----------------------------------------------------------------------------

species = (
    'CELL', 'CELLA', 'CH2OHCHO', 'CHOCHO', 'CH3CHO', 'C6H6O3', 'C2H5CHO',
    'CH3OH', 'CH2O', 'CO', 'GCO', 'CO2', 'H2', 'H2O', 'GCH2O', 'HCOOH',
    'CH2OHCH2CHO', 'CH4', 'GH2', 'CHAR', 'C6H10O5', 'GCOH2', 'GMSW', 'HCE1',
    'HCE2'
)

# time vector to evaluate reaction rates [s]
time = np.linspace(0, 2.0)

# Solve batch reactor system of equations
# -----------------------------------------------------------------------------

def dcdt_debiagi(t, y):
    """
    Rate equations for Debiagi 2018 biomass pyrolysis kinetics.
    """
    R = 1.9859          # universal gas constant [cal/(mol K)]
    T = 773.15          # temperature [K]

    cell = y[0]         # cellulose mass fraction [-]
    cella = y[1]        # active cellulose mass fraction [-]
    gmsw = y[22]        # softwood hemicellulose mass fraction [-]

    # Cellulose reactions and rate constants
    # 1) CELL -> CELLA
    # 2) CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR
    # 3) CELLA -> C6H10O5
    # 4) CELL -> 4.45 H2O + 5.45 CHAR + 0.12 GCOH2 + 0.18 GCH2O + 0.25 GCO + 0.125 GH2 + 0.125 H2
    k1 = 1.5e14 * np.exp(-47_000 / (R * T))
    k2 = 2.5e6 * np.exp(-19_100 / (R * T))
    k3 = 3.3 * T * np.exp(-10_000 / (R * T))
    k4 = 9e7 * np.exp(-31_000 / (R * T))

    # Hemicellulose reactions and rate constants
    # 5) GMSW -> 0.70 HCE1 + 0.30 HCE2
    k5 = 1e10 * np.exp(-31_000 / (R * T))

    # species reaction rate equations where r = dc/dt
    # mass fractions associated with each species are also given
    r_CELL = -k1 * cell
    r_CELLA = k1 * cell - k2 * cella - k3 * cella
    r_CH2OHCHO = k2 * cella * 0.1481
    r_CHOCHO = k2 * cella * 0.0107
    r_CH3CHO = k2 * cella * 0.0462
    r_C6H6O3 = k2 * cella * 0.1944
    r_C2H5CHO = k2 * cella * 0.1254
    r_CH3OH = k2 * cella * 0.0395
    r_CH2O = k2 * cella * 0.0278
    r_CO = k2 * cella * 0.0846
    r_GCO = k2 * cella * 0.008637 + k4 * cell * 0.04319
    r_CO2 = k2 * cella * 0.1167
    r_H2 = k2 * cella * 0.001616 + k4 * cell * 0.001554
    r_H2O = k2 * cella * 0.1033 + k4 * cell * 0.4944
    r_GCH2O = k2 * cella * 0.009259 + k4 * cell * 0.03333
    r_HCOOH = k2 * cella * 0.005677
    r_CH2OHCH2CHO = k2 * cella * 0.02284
    r_CH4 = k2 * cella * 0.004947
    r_GH2 = k2 * cella * 0.001243 + k4 * cell * 0.001554
    r_CHAR = k2 * cella * 0.04889 + k4 * cell * 0.4037
    r_C6H10O5 = k3 * cella
    r_GCOH2 = k4 * cell * 0.02222
    r_GMSW = -k5 * gmsw
    r_HCE1 = k5 * gmsw * 0.7
    r_HCE2 = k5 * gmsw * 0.3

    # system of ODEs
    dcdt = (
        r_CELL, r_CELLA, r_CH2OHCHO, r_CHOCHO, r_CH3CHO, r_C6H6O3, r_C2H5CHO,
        r_CH3OH, r_CH2O, r_CO, r_GCO, r_CO2, r_H2, r_H2O, r_GCH2O, r_HCOOH,
        r_CH2OHCH2CHO, r_CH4, r_GH2, r_CHAR, r_C6H10O5, r_GCOH2, r_GMSW, r_HCE1,
        r_HCE2
    )
    return dcdt

# initial mass fractions [-]
y0 = np.zeros(len(species))
y0[0] = 0.6
y0[22] = 0.4

# solution for system of equations for a batch reactor
sol = solve_ivp(dcdt_debiagi, (time[0], time[-1]), y0, t_eval=time)

# Print
# ----------------------------------------------------------------------------

print('--- Initial ---')
print(f'CELL        {sol.y[0][0]}')
print(f'GMSW        {sol.y[22][0]}')

print('--- Final ---')
for i, sp in enumerate(species):
    print(f'{sp:11} {sol.y[i][-1]:.4f}')

Ideally I would like to define all the reactions in an array, parse that array, and output a matrix that can be used by SciPy solve_ivp to determine the concentration of species at each time step. I'm thinking something like what is shown below:

R = 1.9859      # universal gas constant [cal/(mol K)]
T = 773.15      # temperature [K]

reactions = [
    'CELL -> CELLA',
    'CELLA -> 0.40 CH2OHCHO + 0.03 CHOCHO + 0.17 CH3CHO + 0.25 C6H6O3 + 0.35 C2H5CHO + 0.20 CH3OH + 0.15 CH2O + 0.49 CO + 0.05 GCO + 0.43 CO2 + 0.13 H2 + 0.93 H2O + 0.05 GCH2O + 0.02 HCOOH + 0.05 CH2OHCH2CHO + 0.05 CH4 + 0.1 GH2 + 0.66 CHAR',
    'CELLA -> C6H10O5'
]

names = {
    'CELL': 'C6H10O5',
    'CELLA': 'C6H10O5',
    'GCO': 'CO',
    'GCH2O': 'CH2O',
    'GH2': 'H2',
    'CHAR': 'C'
}

arrhenius = [
    (1.5e14, 0, 47_000),
    (2.5e6, 0, 19_100),
    (3.3, 1, 10_000)
]

How can I define a system of chemical equations which can be implemented with solve_ivp without having to write every rate equation by hand for each chemical species?

$\endgroup$
2
  • $\begingroup$ Is there a reason you don't want to use Cantera to do this (i.e. this is a homework assignment, or... something else)? $\endgroup$ – tpg2114 Oct 30 '19 at 15:40
  • $\begingroup$ I am aware of Cantera but I don't enjoy using its API. The Python interface is not very "pythonic" maybe since it was originally developed in C++. My example is part of a larger project which requires reactor models not available in Cantera. $\endgroup$ – wigging Oct 30 '19 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.