# 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:

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$ ?

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')

#------------------------------------------------------------------------------

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['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.

• Gavin, I'm not following. Each unit of G1 produced contains 5 units of water, according to your diagram. You could say there are 5 moles of water per mole of G1, or 5 molecules of water per "molecule" of G1, or whatever. You could also think of G1 being a molecule that is converted to 5 H2O + 6 char at infinite speed relative to the other system kinetics. Either way, unless I am misunderstanding what you mean by your diagram, G1 = 5 H2O + 6 char. – Curt F. Mar 16 '15 at 14:50
• The only other possibility I can think of is that the system units are currently in mass, but you want to find moles of water? In that case you would have to adjust by the formula weight of water. – Curt F. Mar 16 '15 at 14:53
• You would still calculate the mass balance the same way. Or, if you wanted, you could calculate it as total = cell + cella + char + water + HAA + .... + lvg, that is, with all the components expanded. But you can't count G1, G2, or any of their components more than once. So if you include G1, don't include water and char, etc. – Curt F. Mar 16 '15 at 15:24
• Keep in mind that the products returned from my cell function are in $kg/m^3$. So cell, cella, g1, g2, and lvg returned from the function are in units of $kg/m^3$ not moles. The total for the system should always equal 0.5*700 kg/m^3 (350 kg/m^3). I'm confused on how 5*H2O + 6*Char = G1 if G1 should be some amount as kg/m^3. If H2O is 5*g1 then it would be greater than g1. – wigging Mar 16 '15 at 15:38
• I agree there is some inconsistency then. What is the meaning of the coefficiencts "5", "6", and "0.61", etc. in your diagram? You could interpret them as mass coefficients, in which case you would have to calculate water at 5/11*g1 + 0.83*g2/sum(g2coeffs) and char as 6/11*g1 + 0.61*g2/sum(g2coeffs), and so on. Or you could interpret the coefficients as molar coefficients, in which case to find mass you would have to multiply by FWs. For example, 18*5/11*g1/(18+12) + 0.83*g2/sum(g2coeffs)*18/sum(g2FWs). – Curt F. Mar 16 '15 at 15:42