This example is available as an ipynb (Jupyter Notebook) file in the main GitHub repository at https://github.com/pr-omethe-us/PyKED/blob/master/docs/shock-tube-example.ipynb
The ChemKED file used in this example can be found in the tests directory of the PyKED
repository at https://github.com/pr-omethe-us/PyKED/blob/master/pyked/tests/testfile_st_p5.yaml.
The data in this file comes from Stranic et al., describing
shock-tube ignition delays for tert-butanol.
The file starts with the meta information about the ChemKED file itself:
file-authors:
- name: Morgan Mayer
ORCID: 0000-0001-7137-5721
file-version: 0
chemked-version: 0.0.1
This is followed by the information about the reference for the data, and the information about the experimental apparatus:
reference:
doi: 10.1016/j.combustflame.2011.08.014
authors:
- name: Ivo Stranic
- name: Deanna P. Chase
- name: Joseph T. Harmon
- name: Sheng Yang
- name: David F. Davidson
- name: Ronald K. Hanson
journal: Combustion and Flame
year: 2012
volume: 159
pages: 516-527
experiment-type: ignition delay
apparatus:
kind: shock tube
institution: High Temperature Gasdynamics Laboratory, Stanford University
facility: stainless steel shock tube
This ChemKED file specifies multiple data points with some common
conditions, including a common mixture composition and common definition of
ignition delay. Therefore, a common-properties section is specified, followed by the
datapoints list:
common-properties:
composition: &comp
kind: mole fraction
species:
- species-name: t-butanol
InChI: 1S/C4H10O/c1-4(2,3)5/h5H,1-3H3
amount:
- 0.003333333
- species-name: O2
InChI: 1S/O2/c1-2
amount:
- 0.04
- species-name: Ar
InChI: 1S/Ar
amount:
- 0.956666667
ignition-type: &ign
target: OH*
type: 1/2 max
datapoints:
- temperature:
- 1459 kelvin
ignition-delay:
- 347 us
pressure:
- 1.60 atm
composition: *comp
ignition-type: *ign
equivalence-ratio: 0.5
- temperature:
- 1389 kelvin
ignition-delay:
- 756 us
pressure:
- 1.67 atm
composition: *comp
ignition-type: *ign
equivalence-ratio: 0.5
- temperature:
- 1497 kelvin
ignition-delay:
- 212 us
pressure:
- 1.55 atm
composition: *comp
ignition-type: *ign
equivalence-ratio: 0.5
- temperature:
- 1562 kelvin
ignition-delay:
- 105 us
pressure:
- 1.50 atm
composition: *comp
ignition-type: *ign
equivalence-ratio: 0.5
In this example, we will run constant-volume simulations at each
pressure and temperature condition in the datapoints list. Once again,
the ChemKED file specifies all the information required for the simulations except
for the chemical kinetic model, and Cantera can be used to simulate autoignition. For these simulations, we will be using the butanol mechanism developed by Sarathy et al. and available from the LLNL website. In Python, additional functionality can be imported into a script or session by the import
keyword. Cantera, NumPy, and PyKED must be imported into the session so that we can work with the
code. In the case of Cantera and NumPy, we will use many functions from these libraries, so we
assign them abbreviations (ct and np, respectively) for convenience. From PyKED, we
will only be using the ChemKED class, so this is all that is imported. In this example, we also import
Python's built-in multiprocessing library so that we can run the simulations on multiple cores. We
import the Pool class, which offers facilities to manage a jobs on a pool of processors:
In [1]:
import cantera as ct
import numpy as np
from multiprocessing import Pool
from pyked import ChemKED
Then, we define a function that will be mapped onto each job. This function takes the initial temperature, pressure, and mole fractions as input and returns the simulated ignition delay time, defined as the time the temperature increases by 400 K over the initial temperature, a simplified definition for this example. In general, the user could process the mole fraction of OH\* (provided that the kinetic model includes accurate chemistry for OH\*) to match the definition of ignition delay in the experiments.
In [2]:
# Suppress warnings from loading the mechanism file
ct.suppress_thermo_warnings()
def run_simulation(T, P, X):
gas = ct.Solution('LLNL_sarathy_butanol.cti')
gas.TPX = T, P, X
reac = ct.IdealGasReactor(gas)
netw = ct.ReactorNet([reac])
while reac.T < T + 400:
netw.step()
return netw.time
Then, we load the ChemKED file and generate a list of initial conditions that will be mapped onto
the run_simulation() function. We first define a convenience function to collect the input
from a single datapoint and return a Python tuple with the conditions (tuples are lightweight groups
of data in Python). In this function, we define a conversion between the species names in the ChemKED file and the species in the mechanism using a dictionary. Then, we use the built-in map function to apply the
collect_input function to each of the elements in the ck.datapoints list:
In [3]:
from urllib.request import urlopen
import yaml
st_link = 'https://raw.githubusercontent.com/pr-omethe-us/PyKED/master/pyked/tests/testfile_st_p5.yaml'
with urlopen(st_link) as response:
testfile_st = yaml.safe_load(response.read())
ck = ChemKED(dict_input=testfile_st)
def collect_input(dp):
T_initial = dp.temperature.to('K').magnitude
P_initial = dp.pressure.to('Pa').magnitude
species_conversion = {'t-butanol': 'tc4h9oh', 'O2': 'o2', 'Ar': 'ar'}
X_initial = dp.get_cantera_mole_fraction(species_conversion)
return (T_initial, P_initial, X_initial)
initial_conditions = list(map(collect_input, ck.datapoints))
Finally, we create the processor Pool (with 4 processes) and send the jobs out to run:
In [4]:
with Pool(processes=4) as pool:
ignition_delays = pool.starmap(run_simulation, initial_conditions)
for (T, P, X), tau in zip(initial_conditions, ignition_delays):
print(f'The ignition delay for T_initial={T} K, P_initial={P} Pa is: {tau} seconds')
The simulated ignition delay results are returned in the ignition_delays list. The results
are printed to the screen using a Python formatted string (f-string). In addition, the results can be plotted against the experimental data on an Arrhenius plot
In [5]:
%matplotlib notebook
import matplotlib.pyplot as plt
inv_temps = [1000/i[0] for i in initial_conditions]
exp_ignition_delays = [dp.ignition_delay.to('ms').magnitude for dp in ck.datapoints]
sim_ignition_delays = np.array(ignition_delays)*1.0E3
plt.figure()
plt.scatter(inv_temps, exp_ignition_delays, label='Experimental ignition delays')
plt.scatter(inv_temps, sim_ignition_delays, label='Simulated ignition delays', marker='s')
plt.legend(loc='best')
plt.yscale('log')
plt.ylabel('Ignition delay [ms]')
plt.xlabel('1000/T [1/K]');