{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Chemical equilibrium with fixed volume and internal energy\n", "\n", "

Written by Allan Leal (ETH Zurich) on Jan 7th, 2022.
Last revised on Jan 21st, 2022

\n", "\n", "```{attention}\n", "Always make sure you are using the [latest version of Reaktoro](https://anaconda.org/conda-forge/reaktoro). Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these [update instructions](updating_reaktoro_via_conda) to get the latest version of Reaktoro!\n", "```\n", "\n", "Let's consider the combustion of CH{{_4}} again. However, this time, the combustion process happens in a rigid and adiabatic chamber of 10 cm3. Therefore, we will solve a chemical equilibrium problem in which we specify **volume** and **internal energy** as known properties at the equilibrium state (since these properties must be preserved within the chamber!).\n", "\n", "In this problem, both temperature and pressure are unknown and will be calculated together with the amounts of species that bring the system to a state of chemical equilibrium.\n", "\n", "Let's create our {{ChemicalSystem}} object (one with a single gaseous phase):" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from reaktoro import *\n", "\n", "db = NasaDatabase(\"nasa-cea\")\n", "\n", "gases = GaseousPhase(\"CH4 O2 CO2 CO H2O H2\")\n", "\n", "system = ChemicalSystem(db, gases)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we create an initial chemical state for this system in which CH{{_4}} and O{{_2}} exist in a 10 cm3 volume with mole fractions 0.75 and 0.25 respectively. You'll notice that we accomplish this by setting the amounts of CH{{_4}} and O{{_2}} to 0.75 and 0.25 moles, respectively, followed by **scaling the volume of the chemical state** to the desired 10 cm3." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INITIAL STATE\n", "+-----------------+-------------+------+\n", "| Property | Value | Unit |\n", "+-----------------+-------------+------+\n", "| Temperature | 298.15 | K |\n", "| Pressure | 100000 | Pa |\n", "| Charge: | 0 | mol |\n", "| Element Amount: | | |\n", "| :: H | 0.00121019 | mol |\n", "| :: C | 0.000302547 | mol |\n", "| :: O | 0.000201698 | mol |\n", "| Species Amount: | | |\n", "| :: CH4 | 0.000302547 | mol |\n", "| :: O2 | 0.000100849 | mol |\n", "| :: CO2 | 4.03395e-20 | mol |\n", "| :: CO | 4.03395e-20 | mol |\n", "| :: H2O | 4.03395e-20 | mol |\n", "| :: H2 | 4.03395e-20 | mol |\n", "+-----------------+-------------+------+\n" ] } ], "source": [ "state0 = ChemicalState(system)\n", "state0.temperature(25.0, \"celsius\")\n", "state0.pressure(1.0, \"bar\")\n", "state0.set(\"CH4\", 0.75, \"mol\")\n", "state0.set(\"O2\", 0.25, \"mol\")\n", "state0.scaleVolume(10.0, \"cm3\")\n", "\n", "print(\"INITIAL STATE\")\n", "print(state0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check the volume and internal energy in this initial chemical state. The next code block computes the chemical properties of the system in this initial state (using {{ChemicalProps}}). It also stores the value of volume and internal energy in the variables `V0` and `U0` respectively, which we will use later." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Initial volume of the gases is 1e-05 m3\n", "Initial internal energy of the gases is -23.5698 J\n" ] } ], "source": [ "props0 = ChemicalProps(state0)\n", "\n", "V0 = props0.volume() # the initial volume of the gases\n", "U0 = props0.internalEnergy() # the initial internal energy of the gases\n", "\n", "print(\"Initial volume of the gases is\", V0, \"m3\")\n", "print(\"Initial internal energy of the gases is\", U0, \"J\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our chemical equilibrium problem (to model this combustion process) needs to impose the following properties at equilibrium:\n", "\n", "* volume and\n", "* internal energy.\n", "\n", "The code below creates an {{EquilibriumSpecs}} object to reflect these specifications for the equilibrium constraints, which is then used to create our {{EquilibriumSolver}} object." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "specs = EquilibriumSpecs(system)\n", "specs.volume()\n", "specs.internalEnergy()\n", "\n", "solver = EquilibriumSolver(specs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step is to create an {{EquilibriumConditions}} object in which we specify the values of volume and internal energy at equilibrium (which should be those values at the initial state, and that's why we created variables `V0` and `U0` before!). We do this in the following code block, which also sets lower and upper bounds for temperature and pressure (to avoid negative or extremely large values for these properties during the equilibrium calculation!):" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "conditions = EquilibriumConditions(specs)\n", "conditions.volume(V0)\n", "conditions.internalEnergy(U0)\n", "\n", "conditions.setLowerBoundTemperature(25.0, \"celsius\")\n", "conditions.setUpperBoundTemperature(2000.0, \"celsius\")\n", "\n", "conditions.setLowerBoundPressure(1.0, \"bar\")\n", "conditions.setUpperBoundPressure(1000.0, \"bar\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{tip}\n", "It's good to set lower and upper bounds for temperature and pressure when solving chemical equilibrium problems in which these properties are unknown. In the course of the algorithm execution, these properties could become negative or extremely large, causing errors in the algorithm. At the moment, no default lower and upper bounds are set for temperature and pressure because sensible values do not exist for every possible application. For example, consider we had opted for 0.1 K and 0.1 Pa as default lower bound values for temperature and pressure. What would happen if you relied on a thermodynamic model (e.g., an equation of state) that only works for temperatures and pressures above 10 °C and 1 atm, respectively?! \n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have everything we need now to perform the equilibrium calculation, which is done next:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "FINAL STATE\n", "+-----------------+-------------+------+\n", "| Property | Value | Unit |\n", "+-----------------+-------------+------+\n", "| Temperature | 1080.39 | K |\n", "| Pressure | 583621 | Pa |\n", "| Charge: | 0 | mol |\n", "| Element Amount: | | |\n", "| :: H | 0.00121019 | mol |\n", "| :: C | 0.000302547 | mol |\n", "| :: O | 0.000201698 | mol |\n", "| Species Amount: | | |\n", "| :: CH4 | 0.000128967 | mol |\n", "| :: O2 | 1e-16 | mol |\n", "| :: CO2 | 9.68738e-06 | mol |\n", "| :: CO | 0.000163893 | mol |\n", "| :: H2O | 1.84304e-05 | mol |\n", "| :: H2 | 0.000328729 | mol |\n", "+-----------------+-------------+------+\n" ] } ], "source": [ "state = ChemicalState(state0) # let's create a copy of the initial state\n", "\n", "solver.solve(state, conditions)\n", "\n", "print(\"FINAL STATE\")\n", "print(state)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "tags": [ "remove-cell" ] }, "outputs": [ { "data": { "text/plain": [ "autodiff.real(1e-16, 0)" ] }, "metadata": { "scrapbook": { "mime_prefix": "", "name": "nO2" } }, "output_type": "display_data" } ], "source": [ "from myst_nb import glue\n", "glue(\"nO2\", float(state.speciesAmount(\"O2\")));\n", "options = EquilibriumOptions()\n", "glue(\"epsilon\", options.epsilon);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that O{{_2}} was entirely consumed in the process. Its amount ({glue:}`nO2` moles) is not zero because Reaktoro's equilibrium solver has a default lower bound of {glue:}`epsilon` moles for the species amounts. Thus, the equilibrium calculation converged with O{{_2}} *attached to its lower bound*. \n", "\n", "```{tip}\n", "You can use {{EquilibriumOptions}} to change this default lower bound for species amounts:\n", "\n", "~~~python\n", "options = EquilibriumOptions()\n", "options.epsilon = 1e-30\n", "\n", "solver = EquilibriumSolver(specs)\n", "solver.setOptions(options)\n", "~~~\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The previously printed state does not show volume and internal energy. To verify our equilibrium state stored in `state` has volume and internal energy equal to `V0` and `U0` respectively, we do the following:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Volume at initial state: 1e-05 m3\n", "Volume at final state: 1e-05 m3\n", "Internal energy at initial state: -23.5698 J\n", "Internal energy at final state: -23.5698 J\n" ] } ], "source": [ "V = state.props().volume()\n", "U = state.props().internalEnergy()\n", "\n", "print(\"Volume at initial state:\", V0, \"m3\")\n", "print(\"Volume at final state:\", V, \"m3\")\n", "\n", "print(\"Internal energy at initial state:\", U0, \"J\")\n", "print(\"Internal energy at final state:\", U, \"J\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The verification above demonstrates we found an equilibrium state in which volume and internal energy were preserved.\n", "\n", "Let's now check what is the final temperature and pressure when burning CH{{_4}} in that rigid and adiabatic chamber with our given initial conditions:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Temperature at final state: 807.237 °C\n", "Pressure at final state: 5.83621 bar\n" ] } ], "source": [ "T = units.convert(state.temperature(), \"K\", \"degC\") # convert from K to °C\n", "P = units.convert(state.pressure(), \"Pa\", \"bar\") # convert from Pa to bar\n", "\n", "print(\"Temperature at final state:\", T, \"°C\")\n", "print(\"Pressure at final state:\", P, \"bar\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This concludes this tutorial, which demonstrated how Reaktoro can be used to perform chemical equilibrium calculations with given volume and internal energy. \n", "\n", "Keep on reading to learn how to set your own constraints from scratch!" ] } ], "metadata": { "interpreter": { "hash": "e4e8b2f3ae27709963f14fd23a6560d362beea55eaec742263828e04d814e23c" }, "kernelspec": { "display_name": "Python 3.9.7 64-bit ('reaktoro-jupyter-book': conda)", "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.9.9" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }