Lecture23-Heptane-Ignition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Installing Cantera\n",
    "For this notebook you will need [Cantera](http://www.cantera.org/), an open source suite of object-oriented software tools for problems involving chemical kinetics, thermodynamics, and/or transport processes.\n",
    "Fortunately a helpful chap named Bryan Weber has made Anaconda packages, so to install you can simply type \n",
    "```\n",
    "conda install -c bryanwweber cantera\n",
    "```\n",
    "at your terminal (if you can remember back to when you installed Anaconda!).\n",
    "If you are on Windows you will probably also need to install the Visual C++ Redistributable which you can get [from Microsoft here](https://www.microsoft.com/en-us/download/details.aspx?id=48145)\n",
    "There are other, more difficult, ways to install it. It is also already on the COE computer lab 274 Snell (though there you will have to `pip install jupyter` to get this notebook working)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First, import cantera, with the nickname `ct` to save us some typing later.\n",
    "import cantera as ct\n",
    "\n",
    "# Then the usual suspects:\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heptane combustion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the reduced n-heptane model from LLNL https://combustion.llnl.gov/archived-mechanisms/alkanes/heptane-reduced-mechanism. Save the files alongside this python notebook. These files are in \"CHEMKIN\" format. First, we have to convert them into a format that is usable by Cantera. \n",
    "This may take a while and issue some warnings, but then end by saying `Validating mechanism...PASSED`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from cantera import ck2cti\n",
    "ck2cti.main(['--input=heptanesymp159_mec.txt',\n",
    "             '--thermo=heptanesymp_therm.txt',\n",
    "             '--permissive',\n",
    "             '--output=heptanesymp159.cti'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, this mechanism is more detailed than any we have considered before!\n",
    "Now, let's create a 'Solution' phase in Cantera called `gas` from our mechanism."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas = ct.Solution('heptanesymp159.cti')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's examine some of the reactions and species in the mechanism.\n",
    "\n",
    "This will return the first 10 reactions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas.reaction_equations(np.arange(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this will return a list of the chemical species names:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print \" \".join(gas.species_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Knowing what all those species names mean is a [formidable challenge](http://www.northeastern.edu/comocheng/2014/04/nsf-grant-to-identify-discrepancies/) but we are [making headway](http://www.northeastern.edu/comocheng/2015/05/uscombustionmeeting/) (and more help is welcome).\n",
    "\n",
    "For now, lets loop through all the species looking for ones with 7 carbons and 16 hydrogen atoms, which should be all the isomers of heptane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for species in gas.species():\n",
    "    if species.composition == {'C':7, 'H':16}:\n",
    "        print species.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is only one!\n",
    "Based on the name beginning with 'n' let's assume it represents normal-heptane (all 7 carbons in a single chain with no branching), the fuel we want to simulate. Now we need to find the index number for this species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i_heptane = gas.species_names.index('nc7h16')\n",
    "print(\"heptane is species index {0}\".format(i_heptane))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To specify the state of a system we must supply two intensive variables (temperature, pressure, density, specific entropy, specific enthalpy, specific volume) and the composition (mass or mole fractions). We will set the temperature, pressure, and mole fractions. In cantera, mole fractions are `X` and mass fractions are `Y`. We can then print some properties of our gas system by typing  `gas()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas.TPX = 1000, 10e5, 'nc7h16:1.0'\n",
    "gas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find equilbrium you must specify which two intensive variables to hold constant. We'll find the equilbrium at constant Temperature and Pressure, then print the properties again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas.equilibrate('TP')\n",
    "gas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will recall from Thermodynamics II that a system going to equilibrium at constant T and P should minimize the specific Gibbs free energy of the system. Sure enough, it has gone down (compare the \"Gibbs function\" in the \"1 kg\" columns above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print gas.h - gas.T * gas.s\n",
    "print gas.g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets find the equilibrium composition at 1 bar pressure and a range of temperatures between 100 and 2000 K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "temperatures = np.arange(100,2000,20)\n",
    "# make a big array to store the results in\n",
    "equilibrium_mass_fractions = np.zeros((len(temperatures), gas.n_species))\n",
    "for i, T in enumerate(temperatures):\n",
    "    gas.TP = T, 1e5\n",
    "    gas.equilibrate('TP')\n",
    "    print T,\n",
    "    equilibrium_mass_fractions[i,:] = gas.Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now plot the equilibrium mass fractions as a function of temperature. With 160 lines, let's forgo the legend and instead label the biggest peaks directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(temperatures,equilibrium_mass_fractions)\n",
    "for i, name in enumerate(gas.species_names):\n",
    "    Y = equilibrium_mass_fractions[:,i]\n",
    "    if max(Y)> 0.08:\n",
    "        peakT = temperatures[Y.argmax()]\n",
    "        peakY = max(Y)\n",
    "        plt.text(peakT,peakY, name)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question (a) \n",
    "What do you notice about the species that peaks at 100K, and the ones that peak at 2000K? Can you explain or justify this?"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see some of the complexity hidden at low concentrations, let's plot the y axis on a logarithmic scale:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.semilogy(temperatures,equilibrium_mass_fractions)\n",
    "plt.ylim(1e-30,1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you think about how many reactions are equilibrated, it was remarkably quick!\n",
    "\n",
    "Now we'll add some air, which is mostly nitrogen and oxygen. First of all, find the names of anything with just 2 oxygen atoms or just 2 nitrogen atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for species in gas.species():\n",
    "    if species.composition in [{'O':2}, {'N':2}]:\n",
    "        print species.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look up and store the species indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i_oxygen = gas.species_names.index('o2')\n",
    "print(\"oxygen is species number {0}\".format(i_oxygen))\n",
    "i_nitrogen = gas.species_names.index('n2')\n",
    "print(\"nitrogen is species number {0}\".format(i_nitrogen))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question b)\n",
    "\n",
    "For a \"stoichiometric\" mixture of n-heptane and air (enough oxygen to reach complete combustion) how many moles of heptane and how many moles of nitrogen should you have for one mole of oxygen?  Assume air is 80% nitrogen and 20% oxygen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "oxygen_mole = 1. # moles oxygen\n",
    "nitrogen_mole = ####FIXME\n",
    "heptane_mole = ####FIXME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Some checks\n",
    "assert nitrogen_mole / oxygen_mole == 4, \"Assume air is 80% nitrogen and 20% oxygen\"\n",
    "assert oxygen_mole / heptane_mole == 3+1+3//5*3+8-5//3, \"C7H16 + ?? O2 => 8 H2O + 7 CO2\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use those to make a string for the '`X`' when we set `gas.TPX`. Although we call it a mole fraction, they don't need to add up to one: Cantera will normalize it, preserving the ratios.  Then print it, use it, and check it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_string = 'nc7h16:{0},o2:{1},n2:{2}'.format(heptane_mole, oxygen_mole, nitrogen_mole)\n",
    "print \"The 'X' will be set to {0!r}\".format(X_string)\n",
    "gas.TPX = 1000, 10e5, X_string\n",
    "gas()\n",
    "assert round(gas.concentrations[i_oxygen] / gas.concentrations[i_heptane], 2) == 11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question c)\n",
    "We can do an equilibrium analysis like before, but before you do,\n",
    "starting with a stoichiometric mixture of fuel and air\n",
    "what do you expect the equilibrium composition to mostly consist of?\n",
    "(Imagine all reactions are fast with no barriers)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "temperatures = np.arange(100,2000,20)\n",
    "# make a big array to store the results in\n",
    "equilibrium_mass_fractions = np.zeros((len(temperatures), gas.n_species))\n",
    "for i, T in enumerate(temperatures):\n",
    "    gas.TP = T, 1e5\n",
    "    gas.equilibrate('TP')\n",
    "    print T,\n",
    "    equilibrium_mass_fractions[i,:] = gas.Y\n",
    "plt.plot(temperatures,equilibrium_mass_fractions)\n",
    "for i, name in enumerate(gas.species_names):\n",
    "    Y = equilibrium_mass_fractions[:,i]\n",
    "    if max(Y)> 0.08:\n",
    "        peakT = temperatures[Y.argmax()]\n",
    "        peakY = max(Y)\n",
    "        plt.text(peakT,peakY, name)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Now we are done with equilbria, let's calculate some kinetics!\n",
    "\n",
    "Cantera can do complex networks of reactors with valves, flow controllers, etc.\n",
    "but we will make a simple \"reactor network\" with just one constant volume ideal gas batch reactor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gas.TPX = 800, 10e5, X_string\n",
    "\n",
    "reactor = ct.IdealGasReactor(gas)\n",
    "reactor_network = ct.ReactorNet([reactor])\n",
    "\n",
    "start_time = 0.0  #starting time\n",
    "end_time = 4e-3 # seconds\n",
    "n_steps = 251\n",
    "times = np.linspace(start_time, end_time, n_steps)\n",
    "concentrations = np.zeros((n_steps, gas.n_species))\n",
    "pressures = np.zeros(n_steps)\n",
    "temperatures = np.zeros(n_steps)\n",
    "\n",
    "print_data = True\n",
    "if print_data:\n",
    "    #this just gives headings\n",
    "    print('{0:>10s} {1:>10s} {2:>10s} {3:>14s}'.format(\n",
    "            't [s]','T [K]','P [Pa]','u [J/kg]')) \n",
    "\n",
    "for n, time in enumerate(times):\n",
    "    if time > 0:\n",
    "        reactor_network.advance(time)\n",
    "    temperatures[n] = reactor.T\n",
    "    pressures[n] = reactor.thermo.P\n",
    "    concentrations[n,:] = reactor.thermo.concentrations\n",
    "    if print_data:\n",
    "        print('{0:10.3e} {1:10.3f} {2:10.3f} {3:14.6e}'.format(\n",
    "                 reactor_network.time, reactor.T, reactor.thermo.P, reactor.thermo.u))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot some graphs to see how things look"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.plot(times*1e3, concentrations[:,i_heptane])\n",
    "plt.ylabel(\"Heptane concentration (kmol/m3)\")\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.ylim(0,)\n",
    "plt.show()\n",
    "plt.plot(times*1e3, pressures/1e5)\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.ylabel(\"Pressure (bar)\")\n",
    "plt.show()\n",
    "plt.plot(times*1e3, temperatures)\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.ylabel(\"Temperature (K)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although the timescale is milliseconds instead of hours, that looks remarkably like the thermal runaway reaction that caused the T2 laboratory explosion that we studied last lecture. This time, however, it's not a thermal runaway but a chemical runaway - it's the gradual accumulation of reactive radical species like `OH` that is auto-catalytic.\n",
    "\n",
    "Let's look at some of the other species:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# skip the zeroth species which is nitrogen\n",
    "plt.plot(times*1e3, concentrations[:,1:])\n",
    "plt.ylim(0,)\n",
    "plt.ylabel(\"Concentration\")\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "for i, name in enumerate(gas.species_names):\n",
    "    if i==0: continue\n",
    "    concentration = concentrations[:,i]\n",
    "    peak_concentration = max(concentration)\n",
    "    if peak_concentration > 0.001:\n",
    "        peak_time = times[concentration.argmax()]\n",
    "        plt.text(peak_time*1e3, peak_concentration, name)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.semilogy(times*1e3, concentrations)\n",
    "plt.ylim(1e-15,1)\n",
    "plt.ylabel(\"Concentration\")\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.semilogy(times*1e3, concentrations)\n",
    "plt.ylim(1e-4,1)\n",
    "\n",
    "# Add some labels\n",
    "for t in [1.5, 3]:\n",
    "    i = (times*1e3>t).nonzero()[0][0]\n",
    "    time = times[i]*1e3\n",
    "    for j, name in enumerate(gas.species_names):\n",
    "        concentration = concentrations[i,j]\n",
    "        if concentration > 1e-4:\n",
    "            plt.text(time, concentration, name)\n",
    "plt.ylabel(\"Concentration\")\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to identify when the ignition occurs, so that we could compare our simulation with an experiment.\n",
    "Some experiments monitor the concentration `OH` via laser absorption, which you can  pull from your simulation as `[OH]`, but other studies monitor the luminescence of excited `OH*` decaying to `OH`, which is proportional to the rate of formation (not concentration) of `OH*`, which is predominantly made by reaction of `CH` with `O2`, so it is pretty closely proportional to the product `[CH][O2]`, i.e. \"brightest flash of light\" is propontional to “peak `OH*` emission” which can be modeled as “peak in the product of `[CH]` and `[O2]`”.  Likewise photoemission from creation of excited `CH*` can be modeled reasonably as the product `[C2H][O]`. When modeling an experiment it's important to know precisely what the experimenter measurend and how they defined their derived parameters. For now we'll look for the peak in `OH*` emission:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i_ch = gas.species_index('ch')\n",
    "i_o2 = gas.species_index('o2')\n",
    "excited_oh_generation = concentrations[:,i_o2] * concentrations[:,i_ch]\n",
    "plt.plot(times*1e3, excited_oh_generation)\n",
    "plt.xlabel(\"Time (ms)\")\n",
    "plt.ylabel(\"Excited OH* emission (arbitrary units)\")\n",
    "plt.show()\n",
    "ignition_time = times[excited_oh_generation.argmax()]\n",
    "print \"Ignition delay time is {0} ms\".format(ignition_time * 1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now let's put it all together, into a function that takes temperature, pressure, and stoichiometry, and predicts ignition delay time for n-heptane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_ignition_delay(temperature, pressure = 10.,\n",
    "                       stoichiometry = 1.0, plot = False):\n",
    "    \"\"\"\n",
    "    Get the ignition delay time in miliseconds, at the specified\n",
    "    temperature (K), pressure (bar), and stoichiometry \n",
    "    (stoichiometric = 1.0, fuel-rich > 1.0, oxygen-rich < 1.0).\n",
    "    Default pressure is 10.0 bar, default stoichoimetry is 1.0.\n",
    "    If plot=True then it draws a plot (default is False).\n",
    "    \"\"\"\n",
    "    oxygen_mole = 1. \n",
    "    nitrogen_mole = 4*oxygen_mole\n",
    "    heptane_mole = stoichiometry/11\n",
    "    X_string = 'nc7h16:{0},o2:{1},n2:{2}'.format(heptane_mole, oxygen_mole, nitrogen_mole)\n",
    "\n",
    "    gas.TPX = temperature, pressure*1e5, X_string\n",
    "    reactor = ct.IdealGasReactor(gas)\n",
    "    reactor_network = ct.ReactorNet([reactor])\n",
    "\n",
    "    time = 0.0\n",
    "    end_time = 6e-3 \n",
    "    \n",
    "    # Slightly different this time. We'll use lists, so they can be any length\n",
    "    times = []\n",
    "    concentrations = []\n",
    "    pressures = []\n",
    "    temperatures = []\n",
    "    \n",
    "    print_data = True\n",
    "    while time < end_time:\n",
    "        time = reactor_network.time\n",
    "        times.append(time)\n",
    "        temperatures.append(reactor.T)\n",
    "        pressures.append(reactor.thermo.P)\n",
    "        concentrations.append(reactor.thermo.concentrations)\n",
    "        # take a timestep towards the end_time.\n",
    "        # the size of the step will be determined by the ODE solver\n",
    "        # depending on how quickly things are changing.\n",
    "        reactor_network.step(end_time)\n",
    "    \n",
    "    print \"Reached end time {0} in {1} steps\".format(times[-1], len(times))\n",
    "    # convert lists into arrays\n",
    "    concentrations = np.array(concentrations)\n",
    "    times = np.array(times)\n",
    "    pressures = np.array(pressures)\n",
    "    temperatures = np.array(temperatures)\n",
    "\n",
    "    if plot:\n",
    "        plt.plot(times*1e3, pressures/1e5)\n",
    "        plt.xlabel(\"Time (ms)\")\n",
    "        plt.ylabel(\"Pressure (bar)\")\n",
    "        plt.show()\n",
    "    i_ch = gas.species_index('ch')\n",
    "    i_o2 = gas.species_index('o2')\n",
    "    excited_oh_generation = concentrations[:,i_o2] * concentrations[:,i_ch]\n",
    "    if plot:\n",
    "        plt.plot(times*1e3, excited_oh_generation)\n",
    "        plt.xlabel(\"Time (ms)\")\n",
    "        plt.ylabel(\"Excited OH* emission (arbitrary units)\")\n",
    "        plt.show()\n",
    "    step_with_highest_oh_gen = excited_oh_generation.argmax()\n",
    "    if step_with_highest_oh_gen > 1:\n",
    "        ignition_time_ms = 1e3 * times[step_with_highest_oh_gen]\n",
    "        print \"At {0} K {1} bar, ignition delay time is {2} ms\".format(temperature, pressure, ignition_time_ms)\n",
    "        return ignition_time_ms\n",
    "    else:\n",
    "        print \"At {0} K {1} bar, no ignition detected\".format(temperature, pressure)\n",
    "        return np.infty\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_ignition_delay(1000, plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's repeat it at a range of temperatures and pressures, and plot some graphs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "temperatures = np.linspace(1000,1500.,25)\n",
    "ignition_delay_times = np.zeros_like(temperatures)\n",
    "for P in [10,50]:\n",
    "    for i,T in enumerate(temperatures):\n",
    "        ignition_delay_times[i] = get_ignition_delay(T, P)\n",
    "    plt.semilogy(1000./temperatures, ignition_delay_times, 'o-', label='{0} bar'.format(P))\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel(\"1000K / temperature\")\n",
    "plt.ylabel(\"Ignition delay time (ms)\")\n",
    "plt.ylim(1e-2,)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question c\n",
    "Explain why this look as you would expect from Arrhenius behaviour."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question d\n",
    "Repeat the analysis but going down to 600K (i.e. cover the range 600-1500K).\n",
    "Describe and try to explain what you find."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_ignition_delay(800, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "get_ignition_delay(690, 10, plot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}