{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Non-Ideal Shock Tube Example\n", "## Ignition delay time computations in a high-pressure reflected shock tube reactor\n", " \n", "In this example we will illustrate how to setup and use a constant volume, adiabatic reactor to simulate reflected shock tube experiments. This reactor will then be used to compute the ignition delay of a gas at any temperature and pressure. The example very explicitly follows the form set in batch_reactor_ignition_delay_NTC.pynb, which does very similar calculations, but with an IdealGasReactor. All credit is due to the developer of that example. This example generalizes that work to use a Reactor with no pre-assumed EoS. One can also run ideal gas phases through this simulation, simply by specifying a cti file with that thermodynamic EoS.\n", "\n", "Other than the typical Cantera dependencies, plotting functions require that you have matplotlib installed, and data storing and analysis requires pandas. See https://matplotlib.org/ and http://pandas.pydata.org/index.html, respectively, for additional info.\n", " \n", "The example here demonstrates the calculations carried out by G. Kogekar, et al., \"Impact of non-ideal behavior on ignition delay and chemical kinetics in high-pressure shock tube reactors,\" Combust. Flame., 2017.\n", "\n", "The reflected shock tube reactor is modeled as a constant-volume, adiabatic reactor. The heat transfer and the work rates are therefore both zero. With no mass inlets or exits, the 1st law energy balance reduces to:\n", "\n", "\\begin{equation*}\n", "\\frac{dU}{dt} = \\dot{Q} - \\dot{W} = 0.\n", "\\end{equation*}\n", " \n", "Because of the constant-mass and constant-volume assumptions, the density is also therefore constant:\n", "\n", "\\begin{equation*}\n", "\\frac{d\\rho}{dt} = 0.\n", "\\end{equation*}\n", "\n", "Along with the evolving gas composition, then, the thermodynamic state of the gas is defined by the initial total internal energy $U = mu = m\\sum_k\\left(Y_ku_k\\right)$, where $u_k$ and $Y_k$ are the specific internal energy (kJ/kg) and mass fraction of species $k$, respectively. \n", "\n", "The species mass fractions evolve according to the nety chemical production rates due to homogeneous gas-phase reactions:\n", "\n", "\\begin{equation*}\n", "\\frac{dY_k}{dt} = \\frac{W_k}{\\rho}\\dot{\\omega}_k,\n", "\\end{equation*}\n", "\n", "where $W_k$ is the molecular weight of species $k$ $\\left({\\rm kg}\\,{\\rm kmol}^{-3}\\right)$, $\\rho$ is the (constant) gas-phase density $\\left({\\rm kg}\\,{\\rm m^{-3}}\\right)$, and $\\dot{\\omega}_k$ is the net production rate of species $k$ $\\left({\\rm kmol}\\,{\\rm m^{-3}}\\,{\\rm s^{-1}}\\right)$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Runnning Cantera version: 2.4.0a1\n" ] } ], "source": [ "from __future__ import division\n", "from __future__ import print_function\n", "\n", "import pandas as pd\n", "import numpy as np\n", "\n", "import time\n", "\n", "import cantera as ct\n", "print('Runnning Cantera version: ' + ct.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the gas\n", "\n", "In this example we will choose a stoichiometric mixture of n-dodecane and air as the gas. For a representative kinetic model, we use that developed by Wang, Ra, Jia, and Reitz (https://www.erc.wisc.edu/chem_mech/nC12-PAH_mech.zip) by [H.Wang, Y.Ra, M.Jia, R.Reitz, Development of a reduced n-dodecane-PAH mechanism and its application for n-dodecane soot predictions, $Fuel$ 136 (2014) 25–36].\n", "\n", "To fun a different model or use a different EoS, simply replace this cti file with a different mechanism file." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "**** WARNING ****\n", "For species c5h11, discontinuity in cp/R detected at Tmid = 1000\n", "\tValue computed using low-temperature polynomial: 32.0653\n", "\tValue computed using high-temperature polynomial: 32.1675\n", "\n", "\n", "**** WARNING ****\n", "For species c4h4, discontinuity in cp/R detected at Tmid = 1000\n", "\tValue computed using low-temperature polynomial: 16.6543\n", "\tValue computed using high-temperature polynomial: 16.6734\n", "\n", "\n", "**** WARNING ****\n", "For species A3-, discontinuity in cp/R detected at Tmid = 1000\n", "\tValue computed using low-temperature polynomial: 52.2095\n", "\tValue computed using high-temperature polynomial: 52.364\n" ] } ], "source": [ "gas = ct.Solution('data/WangMechanismRK.cti')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define reactor conditions : temperature, pressure, fuel, stoichiometry" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define the reactor temperature and pressure:\n", "reactorTemperature = 1000 #Kelvin\n", "reactorPressure = 40.0*101325.0 #Pascals\n", "\n", "# Set the state of the gas object:\n", "gas.TP = reactorTemperature, reactorPressure\n", "\n", "# Define the fuel, oxidizer and set the stoichiometry:\n", "gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', oxidizer={'o2':1.0, 'n2':3.76})\n", "\n", "# Create a reactor object and add it to a reactor network\n", "# In this example, this will be the only reactor in the network\n", "r = ct.Reactor(contents=gas)\n", "reactorNetwork = ct.ReactorNet([r])\n", "\n", "# now compile a list of all variables for which we will store data\n", "stateVariableNames = [r.component_name(item) for item in range(r.n_vars)]\n", "\n", "# Use the above list to create a DataFrame\n", "timeHistory = pd.DataFrame(columns=stateVariableNames)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define useful functions" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def ignitionDelay(df, species):\n", " \"\"\"\n", " This function computes the ignition delay from the occurence of the\n", " peak in species' concentration.\n", " \"\"\"\n", " return df[species].argmax()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computed Ignition Delay: 4.093e-04 seconds. Took 3.12s to compute\n" ] } ], "source": [ "#Tic\n", "t0 = time.time()\n", "\n", "# This is a starting estimate. If you do not get an ignition within this time, increase it\n", "estimatedIgnitionDelayTime = 0.005\n", "t = 0\n", "\n", "counter = 1;\n", "while(t < estimatedIgnitionDelayTime):\n", " t = reactorNetwork.step()\n", " if (counter%20 == 0):\n", " # We will save only every 20th value. Otherwise, this takes too long\n", " # Note that the species concentrations are mass fractions\n", " timeHistory.loc[t] = reactorNetwork.get_state()\n", " counter+=1\n", "\n", "# We will use the 'oh' species to compute the ignition delay\n", "tau = ignitionDelay(timeHistory, 'oh')\n", "\n", "#Toc\n", "t1 = time.time()\n", "\n", "print('Computed Ignition Delay: {:.3e} seconds. Took {:3.2f}s to compute'.format(tau, t1-t0))\n", "\n", "# If you want to save all the data - molefractions, temperature, pressure, etc\n", "# uncomment the next line\n", "# timeHistory.to_csv(\"time_history.csv\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot the result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import modules and set plotting defaults" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "plt.rcParams['axes.labelsize'] = 16\n", "plt.rcParams['xtick.labelsize'] = 12\n", "plt.rcParams['ytick.labelsize'] = 12\n", "plt.rcParams['figure.autolayout'] = True" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure illustrating the definition of ignition delay" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "plt.figure()\n", "plt.plot(timeHistory.index, timeHistory['oh'],'-o',color='b',markersize=4)\n", "plt.xlabel('Time (s)',fontname='Times New Roman')\n", "plt.ylabel('$\\mathdefault{OH\\, mass\\, fraction,}\\, Y_{OH}}$',fontname='Times New Roman')\n", "\n", "# Figure formatting:\n", "plt.xlim([0,0.00075])\n", "ax = plt.gca()\n", "font = plt.matplotlib.font_manager.FontProperties(family='Times New Roman',size=14)\n", "ax.annotate(\"\",xy=(tau,0.005), xytext=(0,0.005),arrowprops=dict(arrowstyle=\"<|-|>\",color='r',linewidth=2.0),fontsize=14,)\n", "plt.annotate('Ignition Delay Time (IDT)', xy=(0,0), xytext=(0.00004, 0.00525), family='Times New Roman',fontsize=16);\n", "\n", "for tick in ax.xaxis.get_major_ticks():\n", " tick.label1.set_fontsize(12)\n", " tick.label1.set_fontname('Times New Roman')\n", "for tick in ax.yaxis.get_major_ticks():\n", " tick.label1.set_fontsize(12)\n", " tick.label1.set_fontname('Times New Roman')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Illustration : NTC behavior\n", "In the paper by Kogekar, et al., the reactor model is used to demonstrate the impacts of non-ideal behavior on IDTs in the **N**egative **T**emperature **C**oefficient region, where observed IDTs, counter to intuition, increase with increasing temperature." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define the temperatures for which we will run the simulations" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Make a list of all the temperatures we would like to run simulations at\n", "T = [1800, 1600, 1400, 1200, 1100, 1075, 1050, 1025, 1000, 975, 950, 925, 900, 850, 825, 800,\n", " 750, 700]\n", "\n", "estimatedIgnitionDelayTimes = np.ones(len(T))\n", "\n", "# Set the initial guesses to a common value. We could probably speed up simulations \n", "# by tuning this guess, but as seen in the figure above, the 'extra' time after igntion \n", "# does not add many data points or simulation steps. The time savings would be small.\n", "estimatedIgnitionDelayTimes[:] = 0.005\n", "\n", "# Now create a dataFrame out of these\n", "ignitionDelays = pd.DataFrame(data={'T':T})\n", "ignitionDelays['ignDelay'] = np.nan" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the code above for each temperature, and save the IDT for each." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computed Ignition Delay: 6.354e-07 seconds for T=1800K. Took 2.02s to compute\n", "Computed Ignition Delay: 1.586e-06 seconds for T=1600K. Took 2.11s to compute\n", "Computed Ignition Delay: 5.786e-06 seconds for T=1400K. Took 2.34s to compute\n", "Computed Ignition Delay: 3.911e-05 seconds for T=1200K. Took 2.67s to compute\n", "Computed Ignition Delay: 1.326e-04 seconds for T=1100K. Took 2.94s to compute\n", "Computed Ignition Delay: 1.839e-04 seconds for T=1075K. Took 2.77s to compute\n", "Computed Ignition Delay: 2.533e-04 seconds for T=1050K. Took 2.90s to compute\n", "Computed Ignition Delay: 3.365e-04 seconds for T=1025K. Took 3.04s to compute\n", "Computed Ignition Delay: 4.093e-04 seconds for T=1000K. Took 3.13s to compute\n", "Computed Ignition Delay: 4.289e-04 seconds for T=975K. Took 3.06s to compute\n", "Computed Ignition Delay: 3.911e-04 seconds for T=950K. Took 3.34s to compute\n", "Computed Ignition Delay: 3.407e-04 seconds for T=925K. Took 3.21s to compute\n", "Computed Ignition Delay: 3.145e-04 seconds for T=900K. Took 3.30s to compute\n", "Computed Ignition Delay: 3.233e-04 seconds for T=850K. Took 3.51s to compute\n", "Computed Ignition Delay: 3.439e-04 seconds for T=825K. Took 3.48s to compute\n", "Computed Ignition Delay: 3.852e-04 seconds for T=800K. Took 3.63s to compute\n", "Computed Ignition Delay: 6.824e-04 seconds for T=750K. Took 3.92s to compute\n", "Computed Ignition Delay: 2.056e-03 seconds for T=700K. Took 4.04s to compute\n" ] } ], "source": [ "for i, temperature in enumerate(T):\n", " # Setup the gas and reactor\n", " reactorTemperature = temperature\n", " reactorPressure = 40.0*101325.0\n", " gas.TP = reactorTemperature, reactorPressure\n", " gas.set_equivalence_ratio(phi=1.0, fuel='c12h26', oxidizer={'o2':1.0, 'n2':3.76})\n", " r = ct.Reactor(contents=gas)\n", " reactorNetwork = ct.ReactorNet([r])\n", "\n", " # Create and empty data frame\n", " timeHistory = pd.DataFrame(columns=timeHistory.columns)\n", "\n", " t0 = time.time()\n", "\n", " t = 0\n", " counter = 0\n", " while t < estimatedIgnitionDelayTimes[i]:\n", " t = reactorNetwork.step()\n", " if not counter % 20:\n", " timeHistory.loc[t] = r.get_state()\n", " counter += 1\n", "\n", " tau = ignitionDelay(timeHistory, 'oh')\n", " t1 = time.time()\n", "\n", " print('Computed Ignition Delay: {:.3e} seconds for T={}K. Took {:3.2f}s to compute'.format(tau, temperature, t1-t0))\n", "\n", " ignitionDelays.set_value(index=i, col='ignDelay', value=tau)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure: ignition delay ($\\tau$) vs. the inverse of temperature ($\\frac{1000}{T}$). " ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " this.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('