Detailed Mech Simulation



In [ ]:

    
import os
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import display, Image
%matplotlib inline



In [ ]:

    
from rmgpy.chemkin import *
from rmgpy.species import Species



In [ ]:

    
from afm.canteraModel import Cantera

1. specify condition and simulate



In [ ]:

    
model = 'pdd_2014_pruning4_s4_a3ene_c11'

working_dir = os.path.join('../', 'data', 'pdd_chemistry', 'detailed', model)
chemkin_path = os.path.join(working_dir, 'chem_annotated.inp')
species_dict_path = os.path.join(working_dir, 'species_dictionary.txt')



In [ ]:

    
results_path = os.path.join(working_dir, 'results')
if not os.path.exists(results_path):
    os.mkdir(results_path)



In [ ]:

    
speciesList, reactionList = loadChemkinFile(chemkin_path,
                                            species_dict_path)



In [ ]:

    
speciesDict = {}
for spe in speciesList:
    speciesDict[spe.toChemkin()] = spe



In [ ]:

    
# Find the species: ArCCCCR, RCCCCR, RC
pdd = speciesDict['PDD(1)']

reactorTypeList = ['IdealGasConstPressureTemperatureReactor']
molFracList=[{pdd: 1}]
Tlist = ([673.15],'K')
Plist = ([350],'bar')
reactionTimeList = ([3600*14], 's')



In [ ]:

    
job = Cantera(speciesList=speciesList, reactionList=reactionList, outputDirectory='temp_detailed')
job.loadModel()
job.generateConditions(reactorTypeList, reactionTimeList, molFracList, Tlist, Plist)
alldata = job.simulate()

2. reactant conversion



In [ ]:

    
# prepare moles data
time, dataList, reactionSensitivityData = alldata[0]
TData = dataList[0]
PData = dataList[1]
VData = dataList[2]
total_moles = PData.data*VData.data/8.314/TData.data



In [ ]:

    
pdd_mf = dataList[7].data
print dataList[7].label
pdd_moles = pdd_mf*total_moles
pdd_conv = (pdd_moles[0]-pdd_moles)/pdd_moles[0]



In [ ]:

    
plt.plot(time.data, pdd_conv)
numpy.savetxt(os.path.join(results_path, 'reactant_conv.csv'), (time.data, pdd_conv))

3. molecular weight distribution



In [ ]:

    
molecular_weight_distri = []
for mf in dataList[3:]:
    spe = speciesDict[mf.label]
    molecular_weight_distri.append((spe.molecule[0].getMolecularWeight(), mf.data[-1]))



In [ ]:

    
mws = [tup[0]*1000 for tup in molecular_weight_distri]
molefracs = [tup[1] for tup in molecular_weight_distri]



In [ ]:

    
numpy.savetxt(os.path.join(results_path, 'mwd.csv'), (mws, molefracs))



In [ ]: