Note
Go to the end to download the full example code.
Simulating CH4 catalytic Combustion process#
Catalytic combustion was once an promising technology of ultra low NOx emission for gas power generation. The catalytic combustion process converts the gas-phase reactants to the gas-phase products on the catalyst surface while releasing heat, it consists of three major steps: (1) small fuel species and oxygen get adsorbed onto the catalyst surface, (2) the absorbed species react on the catalyst surface, and (3) the products such as carbon dioxide and water get released from the surface into the gas stream. Catalytic combustion can exhibit something similar to the gas-phase “ignition” phenomenon called “light-off” when the net heat release from the surface reactions suddenly spikes. This is where the catalytic combustion process transitions from being kinetic-limited to being transport-limited. Hence, this is why it is appropriate to use the PFR model in this project. By going through the “surface route”, the combustion process can avoid triggering the NOx formation reaction pathways, and the “flame” is anchored at the “light-off” location on the catalyst surface.
Catalytic combustion applications might utilize a two PFRs in series configuration for large fuel species. The first PFR is for pure gas-phase reactions called the pre-oxidizer, its main purpose is to partially break down the large hydrocarbon fuel species for the downstream catalytic combustor. The second PFR allows both gas-phase and surface reactions where the broken down fuel species are oxidized on the catalyst surface to generate heat.
This project attempts to obtain the solution profiles in order to catch the surface
light-off in the downstream catalytic combustor. Two Chemistry Set are created
in this project. The first mech_no_surface chemistry set is for the first PFR,
the pre_burner, and it is a modified GRI 3.0 combustion mechanism with the addition
of the platinum element (PT). For the downstream cat_combustor, the chemistry set
mech_catalytic includes the surface mechanism “surf_CH4_Cat_Combust_PT.inp” as well
as the same gas mechanism for the mech_no_surface chemistry set. The gas-phase
solution at the exit of the pre_burner can be applied to set up the inlet stream
of the cat_combustor as long as they use the same gas-phase mechanism.
Import PyChemkin packages and start the logger#
from pathlib import Path
import time
import ansys.chemkin.core as ck # Chemkin
from ansys.chemkin.core import Color
from ansys.chemkin.core.inlet import Mixture, Stream
from ansys.chemkin.core.logger import logger
# Chemkin PFR model (steady-state)
from ansys.chemkin.core.flowreactors.PFR import PFREnergyConservation as Pfr
import matplotlib.pyplot as plt # plotting
import numpy as np # number crunching
# check working directory
current_dir = str(Path.cwd())
logger.debug("working directory: " + current_dir)
# set verbose mode
ck.set_verbose(True)
# set interactive mode for plotting the results
# interactive = True: display plot
# interactive = False: save plot as a PNG file
global interactive
interactive = True
Create the chemistry sets#
This example uses the methane catalytic combustion mechanism on platinum that
comes with the standard Chemkin example,
reactor_network__two_stage_catalytic_combustor. The reaction mechanism
consists of two parts: the gas-phase mechanism, chem_GRImech3_PT.inp,
describing the gas-phase combustion of methane in air; and the surface mechanism,
surf_CH4_Cat_Combust_PT.inp, describing the reactions lead to the light-off
(surface ignition) of the adsorbed methane and the adsorbed oxygen (or surface C
and O atoms) on the platinum catalyst.
The mechanisms are provided in the examplesdata folder of the local
PyChemkin installation. You must provide the mechanism file names,
chem_GRImech3_PT.inp and surf_CH4_Cat_Combust_PT.inp with the correct path
example_mech_data to the Chemistry Set before pre-processing. The thermodynamic
data file from the default Ansys Chemkin installation will be used.
data_dir = Path(ck.ansys_dir) / "reaction" / "data"
# find the location of this example py file
try:
script_dir_obj = Path(__file__).parent.resolve()
except NameError:
script_dir_obj = Path(current_dir) / ".." / "data"
script_dir = str(script_dir_obj)
# use the relative path to locate the example mechanism folder
example_mech_data = script_dir_obj / ".." / "data"
print(f"example data = {str(example_mech_data)}")
Set up and run the gas phase only reactor#
Create a chemistry set based on the GRI 3.0 methane combustion mechanism with no surface mechanism input.
# gas-phase chemistry only
mech_no_surface = ck.Chemistry(label="GRI3")
# set mechanism input files
# including the full file path is recommended
mech_no_surface.chemfile = str(data_dir / "grimech30_chem.inp")
mech_no_surface.thermfile = str(data_dir / "grimech30_thermo.dat")
Preprocess the gas-phase only chemistry set#
preprocess the mechanism files.
_ = mech_no_surface.preprocess()
Set up the lean premixed fuel-air stream#
Instantiate a stream feed for the inlet gas mixture. The stream is a mixture with the addition of the inlet flow rate. You specify inlet gas properties in the same way as you set up a mixture.
Use the x method of the Stream object to specify the natural gas
composition with a recipe. In this project, the sccm method,
[cm3/min] @ the standard condition, is used set the inlet volumetric flow rate.
The air composition is copied from the pre-defined Air object in Pychemkin.
To create the premixed fuel-air inlet stream with the given equivalence ratio, the
x_by_equivalence_ratio method is applied.
Once the inlet stream lean_premix is created correctly, it is used to instantiate
the gas-only PFR pre-burner. Provide the necessary input parameters and run the
PFR. Use the get_last_solution_mixture() method to get the solution Stream
preburner_exhaust at the exit of the pre_burner. preburner_exhaust will
be used to define the inlet stream properties of the down stream catalytic combustor
cat_combustor.
# set the fuel composition and conditions
natural_gas = Mixture(mech_no_surface)
natural_gas.pressure = 2.5 * ck.P_ATM
natural_gas.temperature = 890.0
natural_gas.x = [
("CH4", 0.971),
("C2H6", 0.0092),
("C3H8", 0.0051),
("CO2", 0.0053),
("N2", 0.0094),
]
# set the air composition and conditions
air = Mixture(mech_no_surface)
air.pressure = natural_gas.pressure
air.temperature = natural_gas.temperature
air.x = ck.Air.x()
# create the lean premixed natural gas-air mixture for the pre-burner
lean_premix = Stream(mech_no_surface)
lean_premix.pressure = natural_gas.pressure
lean_premix.temperature = natural_gas.temperature
# set stream velocity [cm/sec]
lean_premix.velocity = 60.0
# set equivalence ratio
equiv = 0.3
# products from the complete combustion of the natural gas and the air
products = ["CO2", "H2O", "N2"]
# species mole fractions of added/inert mixture.
# You can also create an additives mixture here.
add_frac = np.zeros(mech_no_surface.kk, dtype=np.double) # no additives: all zeros
# use the equivalence ratio to set the stream composition
ierror = lean_premix.x_by_equivalence_ratio(
mech_no_surface, natural_gas.x, air.x, add_frac, products, equivalenceratio=equiv
)
# check fuel-oxidizer mixture creation status
if ierror != 0:
print("Error: Failed to create the fuel-oxidizer mixture.")
exit()
# instantiate the pre-burner object with the lean premixed natural gas-air stream
pre_burner = Pfr(lean_premix, label="Pre-Burner")
# pre-burner parameters
# reactor length [cm]
pre_burner.length = 100.0
# reactor diameter [cm]
pre_burner.diameter = 12.0
# solver parameters
pre_burner.force_nonnegative = True
# set the start wall time
start_time = time.time()
# run the pre-burner PFR model
runstatus = pre_burner.run()
# check run status
if runstatus != 0:
# Run failed.
print(Color.RED + ">>> Run failed. <<<", end=Color.END)
exit()
# Run succeeded.
print(Color.GREEN + ">>> Run completed. <<<", end=Color.END)
# post process the homogeneous stage solution
pre_burner.process_solution()
# get the number of solution time points
n_points = pre_burner.getnumbersolutionpoints()
# the exit/last solution stream from the pre-burner becomes the feed stream of the
# downstream catalytic combustor
preburner_exhaust = pre_burner.get_last_solution_mixture()
# list the composition of the unburned fuel-air mixture
print(f"preburner exhaust temperature = {preburner_exhaust.temperature} [K]")
print(f"preburner exhaust mass flow rate = {pre_burner.mass_flowrate} [g/sec]")
# preburner_exhaust.list_composition(mode="mole")
# save the mass flow rate [g/cm]
preburner_exhaust_mflrt = pre_burner.mass_flowrate
Set up and run the catalytic reactor#
Create a new chemistry set based on the GRI 3.0 methane combustion mechanism with the addition of the catalytic combustion surface mechanism.
mech_catalytic = ck.Chemistry(label="catalytic")
# set mechanism input files
# including the full file path is recommended
local_data_dir = "CH4_Cat_Combust"
mech_catalytic.chemfile = str(
example_mech_data / local_data_dir / "chem_GRImech3_PT.inp"
)
mech_catalytic.surffile = str(
example_mech_data / local_data_dir / "surf_CH4_Cat_Combust_PT.inp"
)
mech_catalytic.thermfile = str(data_dir / "grimech30_thermo.dat")
Preprocess the catalytic combustion chemistry set#
preprocess the mechanism files.
_ = mech_catalytic.preprocess()
Set up the inlet stream for the catalytic combustor#
Use the exhaust mixture from the pre-burner to define the inlet
stream cat_feed for the catalytic combustor cat_combustor.
Even though the Chemistry Sets for the two PFR have the same
gas phase mechanism, you should avoid applying or copying the
preburner_exhaust object to instantiate cat_combustor because
preburner_exhaust is obtained from the pre_burner which does
not have any surface chemistry data. You should assign individual properties
of the cat_feed from the preburner_exhaust.
cat_feed = Stream(mech_catalytic)
# use the exhaust mixture (the solution) from the gas-phase only pre-burner
# to set up the feed stream properties of the catalytic combustor
cat_feed.x = preburner_exhaust.x
cat_feed.pressure = preburner_exhaust.pressure
cat_feed.temperature = preburner_exhaust.temperature
# the mass flow rate of the feed stream to the catalytic burner
# equals the mass flow rate out of the gas phase only pre-burner
cat_feed.mass_flowrate = preburner_exhaust_mflrt
# list the composition of the unburned fuel-air mixture
print(f"cat burner feed temperature = {cat_feed.temperature} [K]")
print(f"cat burner feed mass flow rate = {cat_feed.mass_flowrate} [g/sec]")
# cat_feed.list_composition(mode="mole")
Create the PFR to predict the gas composition of the outlet stream#
Instantiate the PFR cat_combustor as a PFREnergyConservation
object.
cat_combustor = Pfr(cat_feed, label="cat_combustor")
Set up additional reactor model parameters#
Before you can run the simulation, you must provide reactor parameters, solver controls, and output instructions. For PFR, you must provide either the reactor length, the reactor diameter (or, separately, the flow area and the surface area per reactor length). You can also use profiles to assign certain reactor conditions, such as the reactor temperature, the reactor pressure, and the surface area.
# catalytic combustor parameters
# reactor length [cm]
cat_combustor.length = 10.0
# reactor cross-sectional flow area [cm2]
cat_combustor.flowarea = 5.0
# reactor chemically active surface area per reactor length [cm2/cm]
# cat_combustor.area = 3000.0
area_prof_x = [0.0, 0.5, 1.0, 25.0]
area_prof_a = [2.0e2, 1.0e3, 1.0e4, 1.0e4]
cat_combustor.set_surfacearea_profile(area_prof_x, area_prof_a)
# specify the preactor ressure profile to simulate the pressure drop across
# the reactor length
# x: [cm]
pres_prof_x = [0.0, 25.0]
# pressure: [dynes/cm2]
pres_prof_p = [2.5 * ck.P_ATM, 2.48 * ck.P_ATM]
cat_combustor.set_pressure_profile(pres_prof_x, pres_prof_p)
Set up surface chemistry parameters#
When the Chemistry Set of the Stream includes surface chemistry,
you need to provide some essential surface chemistry related model parameters.
For example, you need to specify the activate surface area as well as the
temperature for each material defined in the surface mechanism. You can also provide
the surface coverage, site fractions and bulk activities, of the materials at the
PFR inlet to improve convergence performance. Use the set_site_fraction() and
the set_bulk_activity() method to specify estimated site fractions and bulk
activities at the reactor entrance (x = 0). By default, all site species have the
same fraction and all bulk species activity is set to 1.
You can use surface chemistry methods such as get_material_names()` and
``get_site_species_names() to obtain the material names and the surface site
species symbols, respectively.
# get the number of surface material defined in the surface mechanism
n_material = cat_combustor.get_numb_material
# get all surface material names
cat_material_names = cat_combustor.get_material_names()
# set surface area fraction (by default, materials have the same area)
# site species symbols: list[list[str]]
site_names = cat_combustor.get_site_species_names()
# bulk species symbols: list[list[str]]
bulk_names = cat_combustor.get_bulk_species_names()
# list the material names and the site species and the bulk species symbols
# per material
print(f"number of surface materials = {n_material}")
for m in range(n_material):
print(f"material names: {cat_material_names[m]}")
print(f" list of site species:\n {site_names[m]}")
print(f" list of bulk species:\n {bulk_names[m]}")
# set up surface coverage of all surface materials
site_recipe = [[("O(S)", 1.0), ("PT(S)", 0.0)]]
bulk_recipe = [[("PT(B)", 1.0)]]
# set the surface condition of the surface material
for i, mname in enumerate(cat_material_names):
# set area fraction [-] (optional: by default, materials have the same area)
# set material temperature [K] if it is different from the gas temperature
# (optional: by default, it is the same as the mixture temperature)
# cat_combustor.set_material_temperature(mname, 1000.0)
# set reactor surface coverage of the material
# site fractions (optional: by default, all sites have the same fraction)
cat_combustor.set_site_fraction(mname, site_recipe[i])
# bulk activities (optional: by default, all bulk species activities = 1.0)
cat_combustor.set_bulk_activity(mname, bulk_recipe[i])
# scale surface reaction rates
cat_combustor.surface_ratemultiplier = 1.0
Set solver controls#
You can overwrite the default solver controls by using solver-related methods, such as those for tolerances.
# set solver the tolerances
cat_combustor.tolerances = (1.0e-8, 1.0e-6)
# solver non-negative mode
cat_combustor.force_nonnegative = False
# transient solver max time step size [cm]
# cat_combustor.set_solver_max_timestep_size(1.0e-2)
# set adaptive solution saving distance
cat_combustor.adaptive_solution_saving(mode=True, steps=20)
Run the catalytic combustor reactor model#
Once all the necessary input parameters are provided, run the cat_combustor with
the run() method. Use the process_solution() method to retrieve the raw
solution profiles. If you are interested in the exit solution, use the
get_last_solution_mixture() method to get the solution at the exit as a
Stream object combustor_exhaust.
Note
In Pychemkin, the surface material is referred by its name. Instead of
the surface material index, you loop over the surface material names.
Use get_material_names() method to get a list of surface material
names defined in the surface mechanism.
Note
For a site or a bulk species, use any Mixture or Stream created
with the Chemistry Set containing the surface mechanism to get its
species index. There are two type of index for surface species, the
global index includes all species (the gas plus the site and the bulk of
all materials); the local index includes the same type of surface species
(site or bulk) of all materials.
# get gas-species index
i_co = cat_feed.get_specindex("CO")
i_ch4 = cat_feed.get_specindex("CH4")
i_no = cat_feed.get_specindex("NO")
# get surface species index: <global index> <local index>
# global index includes all species: the gas, the site of all materials, and the bulk
# of all materials
# local index includes the same type of species (site or bulk) of all materials
i_h_s_global, i_h_s_local = cat_feed.get_surf_specindex("H(S)")
i_pt_s_global, i_pt_s_local = cat_feed.get_surf_specindex("PT(S)")
# run the catalytic combustor PFR model
runstatus = cat_combustor.run()
# check run status
if runstatus != 0:
# Run failed.
print(Color.RED + ">>> Run failed. <<<", end=Color.END)
exit()
# Run succeeded.
print(Color.GREEN + ">>> Run completed. <<<", end=Color.END)
# compute the total runtime
runtime = time.time() - start_time
print(f"Total simulation duration: {runtime} [sec].")
# postprocess the solution profiles
cat_combustor.process_solution()
# get the exit/last solution stream from the catalytic combustor
combustor_exhaust = cat_combustor.get_last_solution_mixture()
# display the gas composition at the reactor exit
combustor_exhaust.list_composition(mode="mole")
# get the number of data points in the solution profiles
n_points = cat_combustor.getnumbersolutionpoints()
print(f"number of solution time points = {n_points}")
# store solution profiles
# distance from the reactor entrance [cm]
distance = cat_combustor.get_solution_variable_profile("distance") #
# temperature solution profile along the reactor length [K]
temp_profile = cat_combustor.get_solution_variable_profile("temperature")
# total heat release rate from the gas-phase reactions [erg/sec]
gashrr_profile = cat_combustor.get_solution_variable_profile("gashrr")
gashrr_profile /= 1.0e3 * ck.ERGS_PER_JOULE
# total heat release rate from all surface reactions [erg/sec]
surfhrr_profile = cat_combustor.get_solution_variable_profile("surfhrr")
surfhrr_profile /= 1.0e3 * ck.ERGS_PER_JOULE
# gas phase CO solution profile [mass fraction]
co_profile = cat_combustor.get_solution_variable_profile("CO")
# gas phase CH4 solution profile [mass fraction]
ch4_profile = cat_combustor.get_solution_variable_profile("CH4")
# gas phase NO solution profile [mass fraction]
no_profile = cat_combustor.get_solution_variable_profile("NO")
# surface PT(S) site fraction solution profile [site fraction]
pt_s_profile = cat_combustor.get_solution_variable_profile("PT(S)")
# surface OH(S) site fraction solution profile [site fraction]
oh_s_profile = cat_combustor.get_solution_variable_profile("OH(S)")
# clean up
ck.done()
Plot the catalytic combustor solution profiles#
Plot the solution profiles along the cat_combustor length.
You could observe the surface “light-off” point around x = 1.0 by the spikes
in the solution profiles. The much bigger peak value of the surface heat release
rate in comparison to the peak gas-phase heat release rate indicates the
combustion process takes place at the catalyst surface.
plt.subplots(2, 2, sharex="col", figsize=(12, 6))
plt.suptitle("PFR Catalytic Combustor Solution", fontsize=16)
plt.subplot(221)
plt.plot(distance, temp_profile, "r-")
plt.ylabel("Temperature [K]")
plt.subplot(222)
plt.plot(distance, gashrr_profile, "b-")
plt.plot(distance, surfhrr_profile, "r-")
plt.legend(["Gas Phase", "Surface"], loc="upper right")
plt.ylabel("Total Heat Release Rate [kJ/sec]")
plt.subplot(223)
plt.plot(distance, co_profile, "g-")
plt.plot(distance, ch4_profile, "r--")
plt.plot(distance, no_profile, "b:")
plt.legend(["CO", "CH4", "NO"], loc="upper right")
plt.yscale("log")
plt.xlabel("Distance [cm]")
plt.ylabel("Gas Mass Fraction")
plt.subplot(224)
plt.plot(distance, pt_s_profile, "b-")
plt.plot(distance, oh_s_profile, "g:")
plt.legend(["PT(S)", "OH(S)"], loc="lower right")
plt.yscale("log")
plt.xlabel("Distance [cm]")
plt.ylabel("Surface Site Fraction")
# plot results
if interactive:
plt.show()
else:
plt.savefig("plot_PFR_catalytic_combustion.png", bbox_inches="tight")