{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# RCM modeling with varying reactor volume\n", "\n", "This example is available as an ipynb (Jupyter Notebook) file in the main GitHub repository at https://github.com/pr-omethe-us/PyKED/blob/main/docs/rcm-example.ipynb\n", "\n", "The ChemKED file that will be used in this example can be found in the\n", "`tests` directory of the PyKED\n", "repository at https://github.com/pr-omethe-us/PyKED/blob/main/tests/testfile_rcm.yaml.\n", "Examining that file, we find the first section specifies the information about\n", "the ChemKED file itself:\n", "\n", "```yaml\n", "file-authors:\n", " - name: Kyle E Niemeyer\n", " ORCID: 0000-0003-4425-7097\n", "file-version: 0\n", "chemked-version: 0.4.0\n", "```\n", "\n", "Then, we find the information regarding the article in the literature from which\n", "this data was taken. In this case, the dataset comes from the work of\n", "[Mittal et al.][Mittal2006a]:\n", "\n", "```yaml\n", "reference:\n", " doi: 10.1002/kin.20180\n", " authors:\n", " - name: Gaurav Mittal\n", " - name: Chih-Jen Sung\n", " ORCID: 0000-0003-2046-8076\n", " - name: Richard A Yetter\n", " journal: International Journal of Chemical Kinetics\n", " year: 2006\n", " volume: 38\n", " pages: 516-529\n", " detail: Fig. 6, open circle\n", "experiment-type: ignition delay\n", "apparatus:\n", " kind: rapid compression machine\n", " institution: Case Western Reserve University\n", " facility: CWRU RCM\n", "```\n", "\n", "[Mittal2006a]: https://doi.org/10.1002/kin.20180\n", "\n", "Finally, this file contains just a single datapoint, which describes the experimental\n", "ignition delay, initial mixture composition, initial temperature, initial pressure,\n", "compression time, ignition type, and volume history that specifies\n", "how the volume of the reactor varies with time, for simulating the compression\n", "stroke and post-compression processes:\n", "\n", "```yaml\n", "datapoints:\n", "- temperature:\n", " - 297.4 kelvin\n", " ignition-delay:\n", " - 1.0 ms\n", " pressure:\n", " - 958.0 torr\n", " composition:\n", " kind: mole fraction\n", " species:\n", " - species-name: H2\n", " InChI: 1S/H2/h1H\n", " amount:\n", " - 0.12500\n", " - species-name: O2\n", " InChI: 1S/O2/c1-2\n", " amount:\n", " - 0.06250\n", " - species-name: N2\n", " InChI: 1S/N2/c1-2\n", " amount:\n", " - 0.18125\n", " - species-name: Ar\n", " InChI: 1S/Ar\n", " amount:\n", " - 0.63125\n", " ignition-type:\n", " target: pressure\n", " type: d/dt max\n", " rcm-data:\n", " compression-time:\n", " - 38.0 ms\n", " time-histories:\n", " - type: volume\n", " time:\n", " units: s\n", " column: 0\n", " volume:\n", " units: cm3\n", " column: 1\n", " values:\n", " - [0.00E+000, 5.47669375000E+002]\n", " - [1.00E-003, 5.46608789894E+002]\n", " - [2.00E-003, 5.43427034574E+002]\n", " ...\n", "```\n", "\n", "The values for the volume history in the `time-histories` key are truncated here to save space. One application of the\n", "data stored in this file is to perform a simulation using [Cantera][Cantera] to\n", "calculate the ignition delay, including the facility-dependent effects represented in the volume\n", "trace. All information required to perform this simulation is present in the ChemKED file, with the\n", "exception of a chemical kinetic model for H2/CO combustion.\n", "\n", "[Cantera]: https://cantera.org\n", "\n", "In Python, additional functionality can be imported into a script or session by the `import`\n", "keyword. Cantera, NumPy, and PyKED must be imported into the session so that we can work with the\n", "code. In the case of Cantera and NumPy, we will use many functions from these libraries, so we\n", "assign them abbreviations (`ct` and `np`, respectively) for convenience. From PyKED, we\n", "will only be using the `ChemKED` class, so this is all that is imported:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import cantera as ct\n", "import numpy as np\n", "\n", "from pyked import ChemKED" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we have to load the ChemKED file and retrieve the first element of the `datapoints`\n", "list. Although this file only encodes a single experiment, the `datapoints` attribute will\n", "always be a list (in this case, of length 1). The elements of the\n", "`datapoints` list are instances of the [`DataPoint`](https://pr-omethe-us.github.io/PyKED/chemked.html#pyked.chemked.DataPoint) class, which we store in the variable\n", "`dp`. To load the YAML file from the web, we also import and use the [PyYAML package](https://github.com/yaml/pyyaml), and the built-in `urllib` package, and use the `dict_input` argument to `ChemKED` to read the information." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from urllib.request import urlopen\n", "\n", "import yaml\n", "\n", "rcm_link = \"https://raw.githubusercontent.com/pr-omethe-us/PyKED/main/tests/testfile_rcm.yaml\"\n", "with urlopen(rcm_link) as response:\n", " testfile_rcm = yaml.safe_load(response.read())\n", "ck = ChemKED(dict_input=testfile_rcm)\n", "dp = ck.datapoints[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The initial temperature, pressure, and mixture composition can be read from the\n", "instance of the `DataPoint` class. PyKED uses instances of the Pint `Quantity` class to\n", "store values with units, while Cantera expects a floating-point value in SI\n", "units as input. Therefore, we use the built-in capabilities of Pint to convert\n", "the units from those specified in the ChemKED file to SI units, and we use the `magnitude`\n", "attribute of the `Quantity` class to take only the numerical part. We also retrieve the\n", "initial mixture mole fractions in a format Cantera will understand:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "T_initial = dp.temperature.to(\"K\").magnitude\n", "P_initial = dp.pressure.to(\"Pa\").magnitude\n", "X_initial = dp.get_cantera_mole_fraction()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With these properties defined, we have to create the objects in Cantera that represent the physical\n", "state of the system to be studied. In Cantera, the `Solution` class stores the thermodynamic,\n", "kinetic, and transport data from an input file in the CTI format. After the `Solution` object\n", "is created, we can set the initial temperature, pressure, and mole fractions using the `TPX`\n", "attribute of the `Solution` class. In this example, we will use the [GRI-3.0](https://combustion.berkeley.edu/gri-mech/version30/text30.html) as the chemical kinetic mechanism for H2/CO combustion. GRI-3.0 is built-in to Cantera, so no other input files are needed." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "gas = ct.Solution(\"gri30.xml\")\n", "gas.TPX = T_initial, P_initial, X_initial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the thermodynamic and kinetic data loaded and the initial conditions defined, we need to\n", "install the `Solution` instance into an `IdealGasReactor` which implements the equations\n", "for mass, energy, and species conservation. In addition, we create a `Reservoir` to represent\n", "the environment external to the reaction chamber. The input file used for the environment,\n", "`air.xml`, is also included with Cantera and represents an average composition of air." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "reac = ct.IdealGasReactor(gas)\n", "env = ct.Reservoir(ct.Solution(\"air.xml\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To apply the effect of the volume trace to the `IdealGasReactor`, a `Wall` must be\n", "installed between the reactor and environment and assigned a velocity. The `Wall` allows the\n", "environment to do work on the reactor (or vice versa) and change the reactor's thermodynamic state;\n", "we use a `Reservoir` for the environment because in Cantera, `Reservoir`s always have a\n", "constant thermodynamic state and composition. Using a `Reservoir` accelerates the solution\n", "compared to using two `IdealGasReactor`s, since the composition and state of the environment\n", "are typically not necessary for the solution of autoignition problems. Although we do not show the\n", "details here, a reference implementation of a class that computes the wall velocity given the volume\n", "history of the reactor is available in [CanSen](https://bryanwweber.github.io/CanSen/), in the\n", "[`cansen.profiles.VolumeProfile`](https://bryanwweber.github.io/CanSen/codedocs.html#cansen.profiles.VolumeProfile) class, which we import here:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from cansen.profiles import VolumeProfile\n", "\n", "exp_time = dp.volume_history.time.magnitude\n", "exp_volume = dp.volume_history.volume.magnitude\n", "keywords = {\"vproTime\": exp_time, \"vproVol\": exp_volume}\n", "ct.Wall(reac, env, velocity=VolumeProfile(keywords));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, the `IdealGasReactor` is installed in a `ReactorNet`. The `ReactorNet`\n", "implements the connection to the numerical solver ([CVODES](https://computation.llnl.gov/projects/sundials/cvodes) is\n", "used in Cantera) to solve the energy and species equations. For this example, it is best practice\n", "to set the maximum time step allowed in the solution to be the minimum time difference in the time array from the volume trace:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "netw = ct.ReactorNet([reac])\n", "netw.set_max_time_step(np.min(np.diff(exp_time)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To calculate the ignition delay, we will follow the definition specified in the ChemKED file for\n", "this experiment, where the experimentalists used the maximum of the time derivative of the pressure\n", "to define the ignition delay. To calculate this derivative, we need to store the state variables and the composition on each time step, so we initialize several Python lists to act as storage:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "time = []\n", "temperature = []\n", "pressure = []\n", "volume = []\n", "mass_fractions = []" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the problem is integrated using the `step` method of the `ReactorNet`. The\n", "`step` method takes one timestep forward on each call, with step size determined by the CVODES\n", "solver (CVODES uses an adaptive time-stepping algorithm). On each step, we add the relevant variables\n", "to their respective lists. The problem is integrated until a user-specified end time, in this case\n", "50 ms, although in principle, the user could end the simulation on any condition\n", "they choose:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "while netw.time < 0.05:\n", " time.append(netw.time)\n", " temperature.append(reac.T)\n", " pressure.append(reac.thermo.P)\n", " volume.append(reac.volume)\n", " mass_fractions.append(reac.Y)\n", " netw.step()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, the user would post-process the information in the pressure list to calculate\n", "the derivative by whatever algorithm they choose. We will plot the pressure versus the time of the simulation using the Matplotlib library:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "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 = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "\n", "plt.plot(exp_time, exp_volume / exp_volume[0], label=\"Experimental volume\", linestyle=\"--\")\n", "plt.plot(time, volume, label=\"Simulated volume\")\n", "plt.legend(loc=\"best\")\n", "plt.ylabel(\"Volume [m^3]\")\n", "plt.xlabel(\"Time [s]\");" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }