Kinetic objects¶

The Kinetic class is responsible for calculating reaction rates and reaction enthalpies when solving the reactor equations.

[1]:

import numpy as np

import reactord as rd

We will use an example to illustrate its functionality. First, we define the kinetic constants. Reaction 1 is reversible and has two kinetic constant independents of temperature and reaction 2 is irreversible, so the pre-exponential factor and $\frac{E}{R}$ of the Arrhenius equations are given.

[2]:

# Reaction 1 (constant with temperature)
kd = 0.1 / 3600 / 0.001  # m⁶/s/mol²
ki = 0.05 / 3600 / 0.001  # 1/s

# Reaction 2 (variable with temperature)
a2 = 0.15 / 3600 / 0.001  # 1/s
e2 = 1000 # K

Reaction rate¶

Next, we set the reaction rates, which are defined by Python function with the format:

def reaction_rate(
    composition: reactord.substance.CompositionalArgument,
    temperature: float,
    kinetics_constants: dict
) -> float:
    # calculation of the reaction rate
    return evaluated_reaction_rate

The rate laws are defined as functions. for example, given the hypothetical reactions:

\[2A + B \leftrightarrow C\]

\[C \rightarrow D\]

Where the first one is a reversible reaction and the second is irreversible, and both are elemental and function of the concentrations of the substrates [mol/m³]. Also, we have the data that, for the first reaction, the rate of consumption of A is:

\[r_A = k_d C_A^2 C_B - k_i C_C\]

And the rate of consumption of C in the second reaction is:

\[r_2 = a_2 e^{\frac{-e_2}{T}} C_C\]

Since we have the reaction rate for the consumption of A in the first reaction. We have to fix the stoichiometric coefficient of A to 1:

\[A + \frac{1}{2} B \leftrightarrow \frac{1}{2} C\]

\[C \rightarrow D\]

Reaction rates are defined positively. ReactorD will automatically detect, for a substrate, if its stoichiometric coefficient is positive or negative. Also, in the case of substrate C, would take into account the production of C by reaction 1 and the consumption of C by reaction 2.

[3]:

def r_rate1(concentrations, temperature, cons):
    kd, ki = cons["kd"], cons["ki"]

    # IMPORTANT! "A", "B", "C" are the names of our substance objects
    ca = concentrations["A"]
    cb = concentrations["B"]
    cc = concentrations["C"]

    return kd * ca ** 2 * cb - ki * cc


def r_rate2(concentrations, temperature, cons):
    a2, e2 = cons["a2"], cons["e2"]

    cc = concentrations["C"]

    return a2 * np.exp(-e2 / temperature) * cc

Note that both functions follow the broadcasting rules of Numpy. This is mandatory.

Now we must define our substances. Imagine that all the substances have the same molar volume and it’s constant with temperature and pressure:

[4]:

def molar_volume(temperature, pressure):
    return np.full(np.size(temperature), 1e-5) # m3/mol

Substance objects are instantiated:

NOTE THAT THE NAME OF THE SUBSTANCES IS THE SAME USED IN THE REACTION RATES!!!!

[5]:

a = rd.Substance(name="A", volume_liquid=molar_volume)
b = rd.Substance(name="B", volume_liquid=molar_volume)
c = rd.Substance(name="C", volume_liquid=molar_volume)
d = rd.Substance(name="D", volume_liquid=molar_volume)

After that, the mixture of the substances is instantiated. In this case, we create an IdealSolution object:

[6]:

mix = rd.mix.IdealSolution([a, b, c, d])

Next, the kinetic object is created:

[7]:

kinetic = rd.Kinetic(
    mix=mix,
    reactions={
        "r1": {"eq": a + 1/2 * b > 1/2 * c, "rate": r_rate1},
        "r2": {"eq": c > d, "rate": r_rate2},
    },
    kinetic_constants={"kd": kd, "ki": ki, "a2": a2, "e2": e2},
    rates_argument="concentration",
)

What happened there?

The kinetic object receives: * A mix object: it uses the mix object to evaluate the properties of the mixture. In this case, molar volume, to calculate concentrations.

A reactions dictionary: In this parameter, we set the stoichiometry of our reactive system. The format of the dictionary is:
```
reaction_dict={"reaction_name": {"eq": algebraic_expression, "rate": rate_function}}
```
The main keys of the dictionary are the reaction names and are merely aesthetic.

Then, in the interior dictionaries we have the "eq" key where we define the stoichiometry of the reaction by algebraic operation with the substance objects previously defined.

The "rate" key is for the reaction rate function of that reaction.

Finally, there is another optional key "DH" to specify a constant reaction enthalpy for that reaction (constant with temperature and pressure). If this key is not specified, ReactorD will try to calculate the standard reaction enthalpy from the formation enthalpies of the substances, and as corresponds, correct the value to temperature and pressure as the mixture object specifies. In this case, our reaction is isothermic.

Since we have two reactions in our system, we specify two reaction_dict.
Kinetic constants: we specify the kinetic constant values with a dictionary. Note that the kinetic constants have the same name that they were used in the reaction rate functions previously defined.
Rate argument: This is the compositional argument of the reaction_rate. The available compositional arguments in ReactorD are:
- Concentration $[\frac {mol} {m^3}]$
- Partial pressure [Pa]
In this case, our reaction rates are a function of the concentration of the substrates, so we choose “concentration”.

Now we can visually represent the kinetic object:

[8]:

kinetic.irepr

$\displaystyle r1: A + 0.5 B \rightarrow 0.5 C$

$\displaystyle r2: C \rightarrow D$

Also, we can obtain the stoichiometric matrix of our system:

[9]:

kinetic.stoichiometry

[9]:

array([[-1. , -0.5,  0.5,  0. ],
       [ 0. ,  0. , -1. ,  1. ]])

This is why the reaction rates are defined positively. ReactorD will apply that stoichiometry matrix to our reaction rates.

Evaluating the kinetic¶

now we can evaluate both reaction rates at once, but we must specify a mole composition for our substances:

[10]:

moles = np.array([10, 10, 5, 5])

mole_fractions = mix.mole_fractions(moles)

temperature = 303.15
pressure = 101325

kinetic.evaluate(mole_fractions, temperature, pressure)

[10]:

array([4.34027777e+11, 3.84700593e+01])

Constant Reaction enthalpy¶

To define a constant reaction enthalpy for our reaction we can set that the first one is endothermic and has a value of 10,000 $\frac{J}{mol_A}$. The second one is exothermic and the value of the reaction enthalpy is -20,000 $\frac{J}{mol_C}$.

[11]:

kinetic = rd.Kinetic(
    mix=mix,
    reactions={
        "r1": {"eq": a + 1/2 * b > 1/2 * c, "rate": r_rate1, "DH": 10_000},
        "r2": {"eq": c > d, "rate": r_rate2, "DH": -20_000},
    },
    kinetic_constants={"kd": kd, "ki": ki, "a2": a2, "e2": e2},
    rates_argument="concentration",
)

We must initiate the reaction enthalpy of our kinetic. This will be automatically done for the reactors that need this value. This feature prevents the reactor search for values like heat capacities maybe not be defined by the user because, for example, the user wants to simulate an isothermic reactor and the heat capacities are not needed.

[12]:

kinetic.set_dh_function()

kinetic.dhs_evaluate(temperature, pressure)

[12]:

array([[ 10000],
       [-20000]])

With constant reaction enthalpies, doesn’t matter the temperature or pressure, the return value will remain constant.

[13]:

temperatures = np.array([300, 400, 500, 600])

pressures = np.array([100_000, 100_000, 100_000, 100_000])

kinetic.dhs_evaluate(temperatures, pressures)

[13]:

array([[ 10000,  10000,  10000,  10000],
       [-20000, -20000, -20000, -20000]])

Variable reaction enthalpies¶

Now a full example with variable reaction enthalpies.

Consider the combustion of methane and ethane, in gas phase, with nitrogen as inert.

\[CH_4 + 2 O_2 \rightarrow CO_2 + 2 H_2O\]

\[C_2H_6 + \frac{7}{2} O_2 \rightarrow 2 CO_2 + 3 H_2O\]

For example, are used toy reaction rates.

For the calculation of the reaction, enthalpies are necessary the standard ideal gas formation enthalpies of all the substrates, and the ideal gas heat capacities. For our luck, we can take them from the thermo library.

[18]:

# =============================================================================
# Substances
# =============================================================================
ch4 = rd.Substance.from_thermo_database("ch4", "methane")
o2 = rd.Substance.from_thermo_database("o2", "oxygen")
co2 = rd.Substance.from_thermo_database("co2", "carbon dioxide")
h2o = rd.Substance.from_thermo_database("h2o", "water")
c2h6 = rd.Substance.from_thermo_database("c2h6", "ethane")
# DON'T FORGET THE INERT ;)
n2 = rd.Substance.from_thermo_database("n2", "nitrogen")

# =============================================================================
# Mixture
# =============================================================================
mix = rd.mix.IdealGas([ch4, o2, co2, h2o, c2h6, n2])

# =============================================================================
# Reaction rates
# =============================================================================
def r1(c, t, cons):
    k1 = cons["k1"]
    c_ch4, c_o2 = c["ch4"], c["o2"]

    return k1 * c_ch4 * c_o2

def r2(c, t, cons):
    k2 = cons["k2"]
    c_c2h6, c_o2 = c["c2h6"], c["o2"]

    return k2 * c_c2h6 * c_o2

# =============================================================================
# Kinetic object
# =============================================================================
kinetic = rd.Kinetic(
    mix=mix,
    reactions={
        "r1": {"eq": ch4 + 2 * o2 > co2 + 2 * h2o, "rate": r1},
        "r2": {"eq": c2h6 + 7/2 * o2 > 2 * co2 + 3 * h2o, "rate": r2}
    },
    kinetic_constants={"k1": 1e-5, "k2": 1e-6},
    rates_argument="concentration"
)

# Start the reaction enthalpy function
kinetic.set_dh_function()

[15]:

kinetic.stoichiometry

[15]:

array([[-1. , -2. ,  1. ,  2. ,  0. ,  0. ],
       [ 0. , -3.5,  2. ,  3. , -1. ,  0. ]])

[19]:

moles = np.array(
    [
        [10, 15, 16], # ch4 moles
        [20, 10, 16], # o2 moles
        [50, 20, 30], # co2 moles
        [0, 20, 15],  # h2o moles
        [40, 20, 30], # c2h6 moles
        [10, 0, 30]   # n2 moles
    ]
)

mole_fractions = mix.mole_fractions(moles)

temperatures = np.array([500, 1000, 2000])
pressures = np.array([100_000, 200_000, 300_000])


kinetic.evaluate(mole_fractions, temperatures, pressures)

[19]:

array([[6.84754527e-05, 1.20128216e-04, 4.43928451e-05],
       [2.73901811e-05, 1.60170955e-05, 8.32365845e-06]])

[20]:

kinetic.dhs_evaluate(temperatures, pressures)

[20]:

array([[ -800832.99199485,  -801075.31914692,  -809803.25400112],
       [-1425727.62523699, -1427790.25945566, -1443262.71596654]])