A Cantera Simulation with Reaction Sensitivity: Comparison with Native RMG Sensitivity Analysis and CHEMKIN Sensitivity Analysis

Note that this requires Cantera with SUNDIALS installed for sensitivity support. If you are using Anaconda, Cantera version >= 2.3.0 is required

Let's check what version of Cantera you have installed



In [ ]:

    
import cantera
print cantera.__version__



In [ ]:

    
from rmgpy.chemkin import *
from rmgpy.tools.canteraModel import *
from rmgpy.tools.plot import parseCSVData
from rmgpy.species import Species
from IPython.display import display, Image

Load the species and reaction from the RMG-generated chemkin file chem_annotated.inp and species_dictionary.txt file found in your chemkin folder after running a job.



In [ ]:

    
speciesList, reactionList = loadChemkinFile('data/minimal_model/chem_annotated.inp',
                                            'data/minimal_model/species_dictionary.txt',
                                           'data/minimal_model/tran.dat')

Set the reaction conditions



In [ ]:

    
# Find the species: ethane and methane
user_ethane = Species().fromSMILES('CC')
user_methane = Species().fromSMILES('C')
speciesDict = getRMGSpeciesFromUserSpecies([user_ethane, user_methane], speciesList)
ethane = speciesDict[user_ethane]
methane = speciesDict[user_methane]
sensitiveSpecies = [ethane, methane]

#reactorTypeList = ['IdealGasReactor']
reactorTypeList = ['IdealGasConstPressureTemperatureReactor']
molFracList=[{ethane: 1}]
Tlist = ([1300],'K')#,1500,2000],'K')
Plist = ([1],'atm')
reactionTimeList = ([0.5], 'ms')



In [ ]:

    
# Create cantera object, loading in the species and reactions
job = Cantera(speciesList=speciesList, reactionList=reactionList, outputDirectory='temp', sensitiveSpecies=sensitiveSpecies)
# The cantera file must be created from an associated chemkin file

# We can either load the Model from the initialized set of rmg species and reactions
job.loadModel()

# Or load it from a chemkin file by uncommenting the following line:
#job.loadChemkinModel('data/minimal_model/chem_annotated.inp',transportFile='data/minimal_model/tran.dat')

# Generate the conditions based on the settings we declared earlier
job.generateConditions(reactorTypeList, reactionTimeList, molFracList, Tlist, Plist)
# Simulate and plot
alldata = job.simulate()
job.plot(alldata)



In [ ]:

    
# Show the plots in the ipython notebook
for i, condition in enumerate(job.conditions):
    print 'Cantera Simulation: Condition {0} Mole Fractions'.format(i+1)
    display(Image(filename="temp/{0}_mole_fractions.png".format(i+1)))
    
    print 'Cantera Simulation: Condition {0} Ethane Reaction Sensitivity'.format(i+1)
    display(Image(filename="temp/{0}_ethane(1)_sensitivity.png".format(i+1)))



In [ ]:

    
# Let's compare against the same simulation in RMG
# Create an input file

input = '''
database(
    thermoLibraries = ['primaryThermoLibrary'],
    reactionLibraries = [],
    seedMechanisms = [],
    kineticsDepositories = 'default',
    kineticsFamilies = 'default',
    kineticsEstimator = 'rate rules',
)

species(
    label = "ethane",
    structure = SMILES("CC"))

species(
    label = "methane",
    structure = SMILES("C"))

simpleReactor(
    temperature = (1300,"K"), 
    pressure = (1,"atm"),
    initialMoleFractions={
        "ethane": 1
    },
    terminationTime = (0.5,"ms"),
    sensitivity=['ethane','methane']
)

model(
    toleranceMoveToCore = 0.04,
)

options(
    saveSimulationProfiles = True,
)

'''
f = open('temp/temp_input.py', 'wb')
f.write(input)
f.close()



In [ ]:

    
from rmgpy.tools.sensitivity import runSensitivity
runSensitivity('temp/temp_input.py', 'data/minimal_model/chem_annotated.inp', 'data/minimal_model/species_dictionary.txt')



In [ ]:

    
print 'RMG Native Simulation: Species Mole Fractions'
display(Image(filename="temp/solver/simulation_1_19.png"))

print 'RMG Native Simulation: Ethane Reaction Sensitivity'
display(Image(filename="temp/solver/sensitivity_1_SPC_1_reactions.png"))



In [ ]:

    
# Let's also compare against the same simulation and sensitivity analysis that was conducted in CHEMKIN
# and saved as a .csv file
time, dataList = parseCSVData('data/minimal_model/chemkin_mole_fractions.csv')
SimulationPlot(xVar=time,yVar=dataList).plot('temp/chemkin_mole_fractions.png')
print 'CHEMKIN Simulation: Mole Fractions'
display(Image(filename="temp/chemkin_mole_fractions.png"))



In [ ]:

    
time, dataList = parseCSVData('data/minimal_model/chemkin_sensitivity_ethane.csv')
ReactionSensitivityPlot(xVar=time,yVar=dataList).barplot('temp/chemkin_ethane_sensitivity.png')
print 'CHEMKIN Simulation: Ethane Reaction Sensitivity'
display(Image(filename="temp/chemkin_ethane_sensitivity.png"))