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

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

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


# 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.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)

1 Answer 1


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


  • $\begingroup$ 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. $\endgroup$
    – Curt F.
    Mar 16, 2015 at 14:50
  • $\begingroup$ 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. $\endgroup$
    – Curt F.
    Mar 16, 2015 at 14:53
  • $\begingroup$ 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. $\endgroup$
    – Curt F.
    Mar 16, 2015 at 15:24
  • $\begingroup$ 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. $\endgroup$
    – wigging
    Mar 16, 2015 at 15:38
  • $\begingroup$ 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). $\endgroup$
    – Curt F.
    Mar 16, 2015 at 15:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.