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