Use biomass pyrolysis kinetics in a Cantera reactor model

[Gavin]

I'm trying to use Cantera with a kinetic scheme for biomass pyrolysis to look at concentration changes over time in a batch reactor. An overview of the kinetics is shown below along with a reference to the paper. Notice that species concentrations are on a mass basis such as kg/m^3.

• wood = biomass which is typically pine
• gas = lumped species which contains light non-condensable gases
• tar = lumped species of condensable pyrolysis vapors
• char = fully pyrolyzed wood, basically carbon

Reference: Colomba Di Blasi. Analysis of Convection and Secondary Reaction Effects Within Porous Solid Fuels Undergoing Pyrolysis. Combustion Science and Technology, vol. 90, pp. 315-340, 1993.

Assuming an initial wood concentration of 1.0, I can solve the system of reaction rate equations with Python and plot the conversion with time as shown below.

Unfortunately, my attempt to use the kinetic scheme with Cantera is giving errors about incompatible phase types. My blasi.cti file contains the following:

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

    #  Phases data

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

    

    

    stoichiometric_solid(

        name = "wood",

        species = "wood",

        density = (700, "kg/m3")

    )

    

    ideal_gas(

        name = "gas",

        species = "gas"

    )

    

    ideal_gas(

        name = "tar",

        species = "tar"

    )

    

    stoichiometric_solid(

        name = "char",

        species = "char",

        density = (110, "kg/m3")

    )

    

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

    #  Species data

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

    

    species(

        name="wood"

    )

    

    species(

        name = "gas"

    )

    

    species(

        name = "tar"

    )

    

    species(

        name = "char"

    )

    

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

    #  Reaction data

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

    

    # Reaction 1

    reaction("wood => gas", [1.4345e4, 0, 88.6])

    

    # Reaction 2

    reaction("wood => tar", [4.125e6, 0, 112.7])

    

    # Reaction 3

    reaction("wood => char", [7.3766e5, 0, 106.5])

    

    # Reaction 4

    reaction("tar => gas", [4.28e6, 0, 108])

    

    # Reaction 5

    reaction("tar => char", [1.0e6, 0, 108])

and the Python file blasi_reactor.py that uses the above cti file is:

    import cantera as ct

    import matplotlib.pyplot as plt

    

    tk = 773.15     # temperature [K]

    p = 101325.0    # pressure [Pa]

    

    gas = ct.Solution('blasi.cti')

    gas.TP = tk, p

    r = ct.IdealGasConstPressureReactor(gas)

    

    sim = ct.ReactorNet([r])

    time = 0.0

    states = ct.SolutionArray(gas, extra=['t'])

    

    for n in range(50):

        time += 1.0

        sim.advance(time)

        states.append(r.thermo.state, t=time)

    

    plt.figure()

    plt.plot(states.t, states.X[:, gas.species_index('wood')])

    plt.plot(states.t, states.X[:, gas.species_index('gas')])

    plt.plot(states.t, states.X[:, gas.species_index('tar')])

    plt.plot(states.t, states.X[:, gas.species_index('char')])

    plt.xlabel('Time [s]')

    plt.ylabel('Concentration [kg/m^3]')

    plt.show()

The error message from Cantera is:

    Traceback (most recent call last):

      File "blasi_cantera.py", line 9, in <module>

        r = ct.IdealGasConstPressureReactor(gas)

      File "interfaces/cython/cantera/reactor.pyx", line 191, in cantera._cantera.Reactor.__init__

      File "interfaces/cython/cantera/reactor.pyx", line 28, in cantera._cantera.ReactorBase.__init__

      File "interfaces/cython/cantera/reactor.pyx", line 199, in cantera._cantera.Reactor.insert

      File "interfaces/cython/cantera/reactor.pyx", line 50, in cantera._cantera.ReactorBase.insert

    cantera._cantera.CanteraError:

    ***********************************************************************

    CanteraError thrown by IdealGasReactor::setThermoMgr:

    Incompatible phase type provided

    ***********************************************************************

How do I define lumped species such as wood, gas, tar, and char with Cantera? Is it even possible to use such a kinetic scheme in Cantera? I typically create my own pyrolysis models using Python but I would like to use the reactor features in Cantera. This would also allow me to compare results between Cantera and my personal Python models.

Note - I looked at the examples on the Cantera documentation website but everything is for well defined gas phase species where you know the elemental composition and NASA coefficients. I also posted this question on Stack Overflow.

[Rodolfo Rodrigues]

Hi Gavin,

You have to improve your CTI file in order to get results. Cantera verifies the element balance for chemical reactions so you have to set the atoms composition of all species as well as h, s, and cp values.

I attached a CTI file sample that uses wood gasification kinetics based on Di Blasi (2004) and Tinaut et al. (2008).

Regards,

Rodolfo

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To post to this group, send email to canter...@googlegroups.com.
Visit this group at https://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.

--

Rodolfo Rodrigues, Dr Eng

  Professor de Engenharia Química

  Universidade Federal do Pampa (UNIPAMPA)

  Bagé, Rio Grande do Sul, Brasil

[Steven DeCaluwe]

Many thanks, Rodolfo, for the detailed and helpful reply!

Gavin - Rodolfo is correct that Cantera insists upon reactions which conserve elements, which in turn requires that species provide elemental information.

If, for some reason, you wanted to take a less detailed approach than Rodolfo’s suggestion, you could provide some ‘dummy’ compositions for each of your lumped species (say, "atoms = ‘C:1’ “ for all of them).  Then your simple kinetic mechanism from below should run, I believe.  The mechanism which Rodolfo shared certainly seems more informative :), but there can also be a benefit to running a quick, simpler model.

Best,

Steven

On Dec 11, 2018, at 5:04 AM, 'Rodolfo Rodrigues' via Cantera Users' Group <canter...@googlegroups.com> wrote:

Hi Gavin,

You have to improve your CTI file in order to get results. Cantera verifies the element balance for chemical reactions so you have to set the atoms composition of all species as well as h, s, and cp values.

I attached a CTI file sample that uses wood gasification kinetics based on Di Blasi (2004) and Tinaut et al. (2008).

Regards,

Rodolfo

Em seg, 10 de dez de 2018 às 01:03, Gavin <wig...@gmail.com> escreveu:

I'm trying to use Cantera with a kinetic scheme for biomass pyrolysis to look at concentration changes over time in a batch reactor. An overview of the kinetics is shown below along with a reference to the paper. Notice that species concentrations are on a mass basis such as kg/m^3.

<kinetics-1.png>

<kinetics-2.png>

• wood = biomass which is typically pine
• gas = lumped species which contains light non-condensable gases
• tar = lumped species of condensable pyrolysis vapors
• char = fully pyrolyzed wood, basically carbon

Reference: Colomba Di Blasi. Analysis of Convection and Secondary Reaction Effects Within Porous Solid Fuels Undergoing Pyrolysis. Combustion Science and Technology, vol. 90, pp. 315-340, 1993.

Assuming an initial wood concentration of 1.0, I can solve the system of reaction rate equations with Python and plot the conversion with time as shown below.

<wood.cti>

[Gavin]

Steven,

I like your suggestion of running a quick simple model. However, setting the atoms to C:1 for each species is giving me the error shown below. I also attached my Python and CTI files.

Thanks for the help,

Gavin

cantera (master *) $ python blasi_reactor.py

Traceback (most recent call last):

  File "blasi_reactor.py", line 9, in <module>

    r = ct.IdealGasConstPressureReactor(gas)
  File "interfaces/cython/cantera/reactor.pyx", line 191, in cantera._cantera.Reactor.__init__
  File "interfaces/cython/cantera/reactor.pyx", line 28, in cantera._cantera.ReactorBase.__init__
  File "interfaces/cython/cantera/reactor.pyx", line 199, in cantera._cantera.Reactor.insert
  File "interfaces/cython/cantera/reactor.pyx", line 50, in cantera._cantera.ReactorBase.insert
cantera._cantera.CanteraError:
***********************************************************************
CanteraError thrown by IdealGasReactor::setThermoMgr:
Incompatible phase type provided
***********************************************************************

[Nick Curtis]

Hi Gavin,

I'm able to get rid of the 'Incompatible phase type provided' error by simply switching the IdealGasConstPressureReactor to a ConstPressureReactor, i.e.:

r = ct.ConstPressureReactor(gas)

However, this change then results in the following exception:

***********************************************************************
CanteraError thrown by CVodesIntegrator::integrate:
CVodes error encountered. Error code: -9
The right-hand side routine failed at the first call.

Exceptions caught during RHS evaluation:
no convergence. dt = 100
Components with largest weighted error estimates:
1: 6.74155e-295
0: 9.63079e-304
2: 1.58101e-313
***********************************************************************

I think the main issue is you need to figure out how to specify all the various species in a single phase (instead of the four separate phases you have now, two stoich_solid's and ideal_gas's).

When you load the model:

gas = ct.Solution('blasi.cti')

This defaults to the first phase listed in the cti file, i.e.:

stoichiometric_solid(

    name = "wood",

    elements = "C",

    species = "wood",

    density = (700, "kg/m3")

)

Which means it has no idea about any of your other species (and is why you're getting incorrect phase errors!).

For instance, if you do:

gas = ct.Solution('blasi.cti')
gas.TP = tk,

p
gas.X = 'wood:1, gas:1'

you'll run into the error:

***********************************************************************

(home) (xenial)nick@localhost:~/Documents/paper_repos/defense$ python test.py       

Traceback (most recent call last):

 File "test.py", line 9, in <module>
   gas.X = 'wood:1, gas:1'
 File "interfaces/cython/cantera/thermo.pyx", line 553, in cantera._cantera.ThermoPhase.X.__set__ (interfaces/cython/cantera/_cantera.cpp:18455)
cantera._cantera.CanteraError:  
***********************************************************************
CanteraError thrown by Phase::setMoleFractionsByName:
Unknown species 'gas'
***********************************************************************

One way to model this system might be using surface chemisty (see for example the catalytic combustion example and corresponding cti file.

Alternatively, you could convert your cti file to contain all species as ideal gas, e.g.:

ideal_gas(

    name = "gas",

    elements = "C",

    species = "tar gas wood char"

)

but this again gives me errors during integration with:

gas = ct.Solution('blasi.cti')

gas.TP = tk, p

gas.X = 'wood:1'

r = ct.IdealGasConstPressureReactor(gas)

***********************************************************************
CanteraError thrown by CVodesIntegrator::integrate:
CVodes error encountered. Error code: -9
The right-hand side routine failed at the first call.

Exceptions caught during RHS evaluation:
ydot contains non-finite elements:

ydot[1] = nan

Components with largest weighted error estimates:
0: 0
1: 0
2: 0
3: 0
4: 0
5: 0
***********************************************************************

that will require debugging.

HTH,

Nick

[Bryan W. Weber]

Gavin,

In the result you showed from using SciPy, how do you model the energy equation? The error that Nick showed (ydot[1] = nan) indicates that the temperature is not being evaluated properly: https://cantera.org/documentation/docs-2.4/sphinx/html/cython/zerodim.html#cantera.Reactor.get_state That's probably because you didn't provide enough information for Cantera to evaluate the temperature, specifically, how the c_p, h, and s vary with temperature.

Best,

Bryan

[Gavin]

Attached is a Python script (blasi.py) that uses the kinetic equations to plot conversion over time from some initial concentration. The only parameters are temperature and the initial concentration of wood. There is no energy equation just the reaction rate equations.

[Gavin]

The attached slide is an example of the material balance for a batch reactor and CSTR reactor at steady-state. The Python example I posted demonstrates the conversion of wood in a batch reactor using the pyrolysis rate equations. Similarly, I'm trying to implement the same thing in Cantera using wood pyrolysis kinetics in a batch reactor.

[Gavin]

Based on the previous posts, I changed the blasi.cti file to the following:

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

#  Phases data

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

ideal_gas(

    elements="C",

    species="wood gas tar char",

    reactions="all"

)

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

#  Species data

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

species(

    name = "wood",

    atoms = "C:1"

)

species(

    name = "gas",

    atoms = "C:1"

)

species(

    name = "tar",

    atoms = "C:1"

)

species(

    name = "char",

    atoms = "C:1"

)

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

#  Reaction data

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

# Reaction 1

reaction("wood => gas", [1.4345e4, 0, 88.6])

# Reaction 2

reaction("wood => tar", [4.125e6, 0, 112.7])

# Reaction 3

reaction("wood => char", [7.3766e5, 0, 106.5])

# Reaction 4

reaction("tar => gas", [4.28e6, 0, 108])

# Reaction 5

reaction("tar => char", [1.0e6, 0, 108])

And in the Python code, I turned off the energy equation:

import cantera as ct

import matplotlib.pyplot as plt

tk = 773.15     # temperature [K]

p = 101325.0    # pressure [Pa]

gas = ct.Solution('blasi.cti')

gas.TP = tk, p

r = ct.IdealGasConstPressureReactor(gas, energy='off')

sim = ct.ReactorNet([r])

time = 0.0

states = ct.SolutionArray(gas, extra=['t'])

for n in range(30):

    time += 1.0

    sim.advance(time)

    states.append(r.thermo.state, t=time)

plt.figure()

plt.plot(states.t, states.X[:, gas.species_index('wood')], label='wood')

plt.plot(states.t, states.X[:, gas.species_index('gas')], label='gas')

plt.plot(states.t, states.X[:, gas.species_index('tar')], label='tar')

plt.plot(states.t, states.X[:, gas.species_index('char')], label='char')

plt.xlabel('Time [s]')

plt.ylabel('Concentration [kg/m^3]')

plt.legend(loc='best')

plt.show()

This approach does not produce any errors, but the generated plot does not show the conversion at each time step.

Why are the lines constant for the entire time span?

[Bryan W. Weber]

Gavin,

You're not setting the initial composition in this example, is that a typo? What is the equilibrium composition (if you do gas.equilibrate('HP') what do you get)?

Best,

Bryan

[Gavin]

I defined the initial concentration in blasi_reactor.py as shown below. But the generated plot is still the same as before (all constant lines). Running gas.equilibrate('HP') in iPython did not return anything.

import cantera as ct

import matplotlib.pyplot as plt

tk = 773.15                         # temperature [K]

p = 101325.0                        # pressure [Pa]

c0 = 'wood:1 gas:0 char:0 tar:0'    # initial concentrations

gas = ct.Solution('blasi.cti')

gas.TPX = tk, p, c0

r = ct.IdealGasConstPressureReactor(gas, energy='off')

sim = ct.ReactorNet([r])

time = 0.0

states = ct.SolutionArray(gas, extra=['t'])

for n in range(25):

    time += 1.0

    sim.advance(time)

    states.append(r.thermo.state, t=time)

plt.figure()

plt.plot(states.t, states.X[:, gas.species_index('wood')], label='wood')

plt.plot(states.t, states.X[:, gas.species_index('gas')], label='gas')

plt.plot(states.t, states.X[:, gas.species_index('tar')], label='tar')

plt.plot(states.t, states.X[:, gas.species_index('char')], label='char')

plt.xlabel('Time [s]')

plt.ylabel('Concentration [kg/m$^3$]')

plt.legend(loc='best')

plt.show()

[Bryan W. Weber]

Hi Gavin,

First, the exact same plot? That indicates that the reactions did happen, they probably just happened much faster than your smallest time step, since clearly the mole fraction of wood is no longer 1.0, indeed, most of the mixture is the gas species. You may want to use the step function to record the reactor state after every time step and stop at 1 s.

Second, the equilibrate function doesn't return anything, but it solves for and sets the state of the Solution to the equilibrium case. I'm curious what the mole fractions are under those conditions. You'll need to print them manually.

Best,

Bryan

[Nick Curtis]

Hi Gavin,

Running the following code:

import cantera as ct

import matplotlib.pyplot as plt

tk = 773.15                         # temperature [K]

p = 101325.0                        # pressure [Pa]

c0 = 'wood:1 gas:0 char:0 tar:0'    # initial concentrations

gas = ct.Solution('blasi.cti')

gas.TPX = tk, p, c0

r = ct.IdealGasConstPressureReactor(gas, energy='off')

sim = ct.ReactorNet([r])

time = 0.0

states = ct.SolutionArray(gas, extra=['t'])

while sim.time <= 1e-5:

    sim.step()

    states.append(r.thermo.state, t=sim.time)

plt.figure()

plt.plot(states.t, states.X[:, gas.species_index('wood')], label='wood')

plt.plot(states.t, states.X[:, gas.species_index('gas')], label='gas')

plt.plot(states.t, states.X[:, gas.species_index('tar')], label='tar')

plt.plot(states.t, states.X[:, gas.species_index('char')], label='char')

plt.xlabel('Time [s]')

plt.ylabel('Concentration [kg/m$^3$]')

plt.legend(loc='best')

plt.show()

with your updated .cti file, yields the attached graph.

So at the very least, the simulation is doing something now.

My advice is to double and triple check the units, both in your CTI file (see here for default units, and this page for setting units) and in your simulation (IIRC, the model you use specifies that the "concentrations" are in kg/m^3, but the way you set them currently using  the mole-fractions likely won't preserve the correct density -- not sure if that matters)

Best,

Nick

[Gavin]

I specified the units in the CTI file which fixed some of the issues I was having.

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

# Phases data

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

units(

    length="m",

    mass="kg",

    quantity="mol",

    time="s",

    energy="J",

    act_energy="kJ/mol",

    pressure="Pa"

)

The Python code for running the above CTI file in a Cantera reactor model is shown below.

import cantera as ct

import matplotlib.pyplot as plt

import numpy as np

import timeit

# Cantera reactor example with Blasi biomass pyrolysis kinetics

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

ti = timeit.default_timer()

tk = 773.15                             # temperature [K]

p = 101325.0                            # pressure [Pa]

c0 = 'wood:1, gas:0, char:0, tar:0'     # initial concentrations [kg/m^3]

gas = ct.Solution('blasi_kinetics.cti')

gas.TPX = tk, p, c0

r = ct.IdealGasConstPressureReactor(gas, energy='off')

sim = ct.ReactorNet([r])

states = ct.SolutionArray(gas, extra=['t'])

# time vector to evaluate reaction rates [s]

time = np.linspace(0, 25, num=1000)

for tm in time:

    sim.advance(tm)

    states.append(r.thermo.state, t=tm)

print(f'time cantera: {timeit.default_timer() - ti:.4g} s')

# Plot

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

def style():

    plt.box(on=False)

    plt.grid(color='0.9')

    plt.tick_params(color='0.9')

    plt.tight_layout()

plt.figure()

plt.plot(states.t, states.X[:, gas.species_index('wood')], label='wood')

plt.plot(states.t, states.X[:, gas.species_index('gas')], label='gas')

plt.plot(states.t, states.X[:, gas.species_index('tar')], label='tar')

plt.plot(states.t, states.X[:, gas.species_index('char')], label='char')

plt.legend(loc='best')

plt.xlabel('Time [s]')

plt.ylabel('Concentration [kg/m$^3$]')

plt.title('Cantera example')

style()

plt.show()

The only remaining problem I'm having is adjusting the initial concentration in the Python script. For example, if I change 

c0 = 'wood:1, gas:0, char:0, tar:0'

to a different concentration of wood such as

c0 = 'wood:2, gas:0, char:0, tar:0'

the resulting plot is the same as if the initial wood concentration was 1. My Cantera example seems to ignore what is defined in c0. So how is the initial concentration being used by Cantera? Do I need to change the CTI file to adjust the initial concentration?

[Abhinav Verma]

These are mole fractions and cantera will normalize. So after normalization both become the same and no wonder you get the same results. 

Hth

Abhi

[Gavin]

I didn't realize Cantera will automatically normalize the mole fractions. I assume it does the same for mass fractions? Since the biomass pyrolysis kinetics are mass based, I changed the example as follows:

y0 = 'wood:1, gas:0, char:0, tar:0'     # initial mass fractions [-]

gas = ct.Solution('blasi_kinetics.cti')

gas.TPY = tk, p, y0

r = ct.IdealGasConstPressureReactor(gas, energy='off')

and changed the y-label to mass fraction instead of concentration:

plt.ylabel('Mass fractions [-]')

Unless there are more suggestions, all of my questions regarding this example have been answered. Thank you everyone for your help.

[Abhinav Verma]

Hi Gavin, 

you are right that cantera will automatically normazle the mass fractions too. 

best regards

Abhi