Fragment Mech Simulation with Reattachment



In [ ]:

    
import os
import matplotlib.pyplot as plt
%matplotlib inline



In [ ]:

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



In [ ]:

    
from afm.simulator import OdeSimulator
import afm.utils
import afm.simulator

1. specify condition and simulate



In [ ]:

    
temperature = 673.15 # K
pressure = 350*3 # bar

initial_mol_fraction = {
    "ArCCCCR":1.0,
    "LCCCCR":1.0,
    "LCCCC":1.0
}
hr = 14
termination_time = 3600*hr # hrs



In [ ]:

    
model = 'two-sided_newcut1'
working_dir = os.path.join('../', 'data', 'pdd_chemistry', model)
chemkin_path = os.path.join(working_dir, 'chem_annotated.inp')
species_dict_path = os.path.join(working_dir, 'species_dictionary.txt')
smiles_dict_path = os.path.join(working_dir, 'fragment_smiles.txt')



In [ ]:

    
ode_simulator = OdeSimulator(chemkin_path,
                             species_dict_path,
                             smiles_dict_path,
                             temperature,
                             pressure)



In [ ]:

    
alldata = ode_simulator.simulate(initial_mol_fraction, termination_time)



In [ ]:

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

2. reactant conversion



In [ ]:

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

moles_dict = {}
for data in dataList[3:]:
    spe_label = data.label
    moles_dict[spe_label] = max(data.data[-1]*total_moles[-1],0)



In [ ]:

    
ArCCCCR_mf = dataList[3].data
print dataList[3].label
ArCCCCR_moles = ArCCCCR_mf*total_moles
ArCCCCR_conv = (ArCCCCR_moles[0]-ArCCCCR_moles)/ArCCCCR_moles[0]



In [ ]:

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

3. molecular weight distribution



In [ ]:

    
from afm.simulator import categorize_fragments



In [ ]:

    
fragmental_weight_distri = ode_simulator.get_molecular_weight_distribution(alldata)



In [ ]:

    
mws = [tup[0]*1000 for tup in fragmental_weight_distri]
moles = [tup[1] for tup in fragmental_weight_distri]

molefracs = moles/sum(moles)



In [ ]:

    
numpy.savetxt(os.path.join(results_path, 'mwd_{0}hr.csv'.format(hr)), (mws, molefracs))



In [ ]: