Calculate species concentration from first-order kinetic reactions

I am trying to determine the concentration of products from the kinetic reactions of wood pyrolysis. The kinetic scheme separates the wood into cellulose, hemicellulose, and lignin. For example, the kinetic reactions for just the cellulose are shown below:


![kinetic reactions](https://i.stack.imgur.com/B8Zgz.png)


The $G1$ and $G2$ terms represent a group of chemical species, $LVG$ represents levoglucosan, and the $K$ terms represent the kinetic rate constants for each reaction.


I am using Python (see code below) to plot the concentration profiles and mass balance of the system. As shown in the plot below, the total mass remains constant (as it should). The plot represents each component as a percent of the original wood concentration which is given as $kg/m^3$. My next step is to determine the concentration ($kg/m^3$) of each chemical species. From the plot, I have the amount of Group 1 ($G1$) produced but I would like to determine the concentrations for $H\_2O$ and $Char$. For example, $5\, H\_2O + 6\, Char = G1$.


How can I determine the concentration (as $kg/m^3$) of each individual species in Groups $G1$ and $G2$ ?


![cellulose reactions](https://i.stack.imgur.com/dFiqp.png)


Python function for the cellulose reactions (see below).


Note that the activation energy units are in $\frac{kcal}{kmol}$ but the concentration of wood provided to the system is a mass basis represented as the density of the wood $\frac{kg}{m^3}$. However, when using the appropriate units for $R$ the rate constant $K$ is in units of $\frac{1}{s}$ where s is seconds. The amount of cellulose provided to the reactions is assumed to be 50% of the original wood as $0.5 \* 700 \frac{kg}{m^3}$. The concentrations returned from the cellulose function `cell` are in terms of $kg/m^3$. The plot (see above) displays the percent of the concentration relative to the original concentration of wood. For example, cellulose in the plot is calculated as $\frac{\rho\_{cell}}{\rho\_{wood}}\*100$ where $\rho$ is the concentration in terms of $kg/m^3$.



```
def cell(T, pw, dt, p):
    """
    T = temperature, K
    pw = vector of initial wood concentration, kg/m^3
    dt = time step, s
    p = total number of time steps
    """

    # array to store species concentrations as a density, kg/m^3
    # row = chemical species
    # column = concentration at time step
    spec = np.zeros([16, p])

    # initial cellulose concentration in wood, kg/m^3
    spec[0] = pw*0.5

    R = 1.987   # universal gas constant, kcal/kmol*K

    # A = pre-factor (1/s) and E = activation energy (kcal/kmol)
    A1 = 4e13;  E1 = 45000  # CELL -> CELLA
    A2 = 0.5e9; E2 = 29000  # CELLA -> products
    A3 = 1.8;   E3 = 10000  # CELLA -> LVG
    A4 = 4e7;   E4 = 31000  # CELL -> 5*H2O + 6*Char

    # reaction rate constant for each reaction, 1/s
    K1 = A1 * np.exp(-E1 / (R * T))         # CELL -> CELLA
    K2 = A2 * np.exp(-E2 / (R * T))         # CELLA -> G2
    K3 = A3 * T * np.exp(-E3 / (R * T))     # CELLA -> LVG
    K4 = A4 * np.exp(-E4 / (R * T))         # CELL -> G1

    # concentrations for each chemical species, kg/m^3
    for i in range(1, p):
        r1 = K1 * spec[0, i-1]
        r2 = K2 * spec[1, i-1]
        r3 = K3 * spec[1, i-1]
        r4 = K4 * spec[0, i-1]
        spec[0, i] = spec[0, i-1] - (r1+r4)*dt          # CELL
        spec[1, i] = spec[1, i-1] + r1*dt - (r2+r3)*dt  # CELLA
        spec[2, i] = spec[2, i-1] + r4*dt               # G1
        spec[3, i] = spec[3, i-1] + r2*dt               # G2
        spec[4, i] = spec[4, i-1] + r3*dt               # LVG

    # return species array concentrations, kg/m^3
    return spec

```

Python script to plot the reactions from the above cellulose function:



```
# Modules
#------------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as py
import funcRanzi as kn

py.close('all')

# Parameters from Papadikis 2010a
#------------------------------------------------------------------------------

rhow = 700  # density of wood, kg/m^3
Tinf = 773  # ambient temp, K

# Initial Calculations
#------------------------------------------------------------------------------

dt = 0.01   # time step, delta t
tmax = 4   # max time, s
t = np.linspace(0, tmax, num=tmax/dt)   # time vector
p = len(t)  # total number of time steps

# Calculate Concentrations of Chemical Species
#------------------------------------------------------------------------------

# vectors for wood, gas, tar, char concentrations as a density, kg/m^3
pw = np.zeros(len(t))   # wood 
pg = np.zeros(len(t))   # gas
pt = np.zeros(len(t))   # tar
pc = np.zeros(len(t))   # char

pw[:] = rhow    # initial wood concentration as density

# array of chemical species concentrations as a density, kg/m^3
spec = kn.cell(Tinf, pw, dt, p)

# concentration as percent relative to original wood, %
cell = spec[0]/rhow*100     # CELL
cella = spec[1]/rhow*100    # CELLA
g1 = spec[2]/rhow*100       # G1
g2 = spec[3]/rhow*100       # G2
lvg = spec[4]/rhow*100      # LVG

total = cell + cella + g1 + g2 + lvg

# Plot Results
#------------------------------------------------------------------------------

py.rcParams['xtick.major.pad'] = 8
py.rcParams['ytick.major.pad'] = 8
py.rcParams['lines.linewidth'] = 2
py.rcParams['axes.grid'] = True

py.figure(1)
py.plot(t, cell, label='cell')
py.plot(t, cella, label='cella')
py.plot(t, g1, label='g1')
py.plot(t, g2, label='g2')
py.plot(t, lvg, label='lvg')
py.plot(t, total, label='total')
py.title('Cellulose Reactions at T = {} K'.format(Tinf))
py.xlabel('Time (s)')
py.ylabel('Conversion (% dry basis)')
py.legend(loc='best', numpoints=1)
py.show()

```



According to the diagram you wrote, `g1` and `g2` are just stoichiometric groupings of certain components. Because there is no modeled back-reaction from those components to the rest of the system, you can find the component amounts by simple stoichiometry.


Asssuming the "stoichiometric" coefficients in the groupings are meant as *mass* coefficients, then the calculations would be:



```
py.plot(t, 6*g1/(6+5)+0.61g2/(sum([0.8, 0.2, 0.1, 0.25, 0.3, 0.21,...etc...]), label='char')
py.plot(t, 5*g1/(6+5)+0.83g2/(sum([0.8, 0.2, 0.1, 0.25, 0.3, 0.21,...etc...]), label='water')
py.plot(t, 0.88g2/(sum([0.8, 0.2, 0.1, 0.25, 0.3, 0.21,...etc...]), label='HAA')
py.plot(t, 0.2g2/sum([0.8, 0.2, 0.1, 0.25, 0.3, 0.21,...etc...]), label='glyox')

```

etc.


 3