Shock tube modeling with constant volume

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 ignition delay for T_initial=1459 K, P_initial=162120.0 Pa is: 0.0006703341602991564 seconds
The ignition delay for T_initial=1389 K, P_initial=169212.75 Pa is: 0.0012624985834179291 seconds
The ignition delay for T_initial=1497 K, P_initial=157053.75 Pa is: 0.0005355282509854273 seconds
The ignition delay for T_initial=1562 K, P_initial=151987.5 Pa is: 0.00043835671071621276 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]');