Chemical equilibrium with constraints
Contents
Chemical equilibrium with constraints#
Written by Allan Leal (ETH Zurich) on Jan 6th, 2022
Attention
Always make sure you are using the latest version of 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 to get the latest version of Reaktoro!
As we mentioned in the previous section, Reaktoro can be used for a wide range of chemical equilibrium problems. In this section, we show:
the fundamentals of specifying different equilibrium constraints when calculating chemical equilibrium; and
how the equilibrium solver can be configured to allow the system to be open to one or more substances.
We learned that a chemical equilibrium solver can be created with the class EquilibriumSolver, as shown in the example below:
from reaktoro import *
db = NasaDatabase("nasa-cea")
gases = GaseousPhase("CH4 O2 CO2 CO H2O H2")
system = ChemicalSystem(db, gases)
solver = EquilibriumSolver(system)
EquilibriumSolver objects created this way (i.e., as in solver = EquilibriumSolver(system)), are limited to chemical equilibrium calculations where:
temperature and pressure are given;
the system is closed (no mass transfer in or out of the system).
Thus, solver = EquilibriumSolver(system) creates a Gibbs energy minimization solver.
This is not always the desired behavior, however. If we want to impose the temperature and volume of the system, we will need to construct a specialized EquilibriumSolver object that can handle these constraint specifications.
Note
Chemical equilibrium problems in which temperature and volume are prescribed belong to the class of Helmholtz energy minimization problems [Smith and Missen, 1982]. Reaktoro, however, does not implement a Helmholtz energy minimization solver. Instead, Reaktoro implements a single parametric Gibbs energy minimization solver that accepts general equilibrium constraints, which can be configured to efficiently solve all other classes of equilibrium problems (including Helmholtz energy minimization problems!).
Specifying constrained properties#
Let’s create an equilibrium solver that can solve problems with prescribed temperature and volume:
specs = EquilibriumSpecs(system)
specs.temperature()
specs.volume()
solver = EquilibriumSolver(specs)
That’s it! We use EquilibriumSpecs to provide the constraint specifications for our EquilibriumSolver.
Tip
Access the link EquilibriumSpecs to find out all currently existing constraints you can impose! It is even possible to define your own constraints, instead of predefined ones, which we’ll cover in a dedicated tutorial given this is a more advanced use case.
It’s now time to demonstrate how we use this specialized equilibrium solver. We want to model the combustion of 1 mol of CH4 and 0.5 mol of O2 in a chamber of 10 cm³ at 1000 °C. We create below an initial chemical state with this gas composition:
state = ChemicalState(system)
state.set("CH4", 1.0, "mol")
state.set("O2", 0.5, "mol")
print("INITIAL STATE")
print(state)
INITIAL STATE
+-----------------+------------+------+
| Property | Value | Unit |
+-----------------+------------+------+
| Temperature | 298.1500 | K |
| Pressure | 1.0000 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 4.0000e+00 | mol |
| :: C | 1.0000e+00 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 1.0000e+00 | mol |
| :: O2 | 5.0000e-01 | mol |
| :: CO2 | 1.0000e-16 | mol |
| :: CO | 1.0000e-16 | mol |
| :: H2O | 1.0000e-16 | mol |
| :: H2 | 1.0000e-16 | mol |
+-----------------+------------+------+
This state is in disequilibrium. We use EquilibriumConditions below to specify the desired conditions of temperature and volume at the equilibrium state, and then equilibrate the state object:
conditions = EquilibriumConditions(specs)
conditions.temperature(1000.0, "celsius")
conditions.volume(10.0, "cm3")
solver.solve(state, conditions)
print("FINAL STATE")
print(state)
FINAL STATE
+-----------------+------------+------+
| Property | Value | Unit |
+-----------------+------------+------+
| Temperature | 1273.1500 | K |
| Pressure | 16603.2449 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 4.0000e+00 | mol |
| :: C | 1.0000e+00 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 7.1576e-01 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: CO2 | 2.2509e-01 | mol |
| :: CO | 5.9148e-02 | mol |
| :: H2O | 4.9067e-01 | mol |
| :: H2 | 7.7814e-02 | mol |
+-----------------+------------+------+
This printed state does not show volume. Let’s check the volume in the ChemicalProps object attached by EquilibriumSolver to the ChemicalState object state at the end of the computation. Let’s also print the computed pressure:
V = state.props().volume() # volume in m³
P = state.pressure() # pressure in Pa
V = units.convert(V, "m3", "cm3") # convert volume from m³ to cm³
P = units.convert(P, "Pa", "GPa") # convert pressure from Pa to GPa
print("Volume at equilibrium state is", V, "cm³!")
print("Pressure at equilibrium state is", P, "GPa!")
Volume at equilibrium state is 10 cm³!
Pressure at equilibrium state is 1.66032 GPa!
Note that the volume at the computed equilibrium state is exactly what we enforced in the EquilibriumConditions object: 10 cm³.
Attention
Don’t forget to set the expected conditions in the EquilibriumConditions object. For every condition marked to be known in the EquilibriumSpecs object (e.g., specs.temperature() and specs.volume()), make sure you give a value for each of them in the associated EquilibriumConditions object (e.g., conditions.temperature(1000.0, "celsius") and conditions.volume(10.0, "cm3"))! Not doing this can cause unexpected behavior or a runtime error.
Did you know?
The Gibbs energy minimization computation performed with:
solver = EquilibriumSolver(system)
solver.solve(state)
is executed in the background with the following code:
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
conditions = EquilibriumConditions(specs)
conditions.temperature(state.temperature())
conditions.pressure(state.pressure())
solver = EquilibriumSolver(specs)
solver.solve(state, conditions)
so you can type a lot less if all you need is equilibrium calculations with specified temperatures and pressures!
Specifying open conditions#
Let’s say we would like to find out how many mols of CH4 must react with 0.5 mol of O2 in the same chamber (with volume 10 cm³ and temperature 1000 °C) to obtain pressure 1 GPa (which, let’s assume, is the maximum pressure supported by the chamber material).
In this new equilibrium calculation, the system is open to the mass transfer of CH4, and this creates a new degree of freedom in the problem: the amount added/removed of this substance. That’s why we can simultaneously impose a value for temperature (1000 °C), volume (10 cm³), and pressure (1 GPa), which would be thermodynamically difficult otherwise.
The code below starts by creating a chemical state in which there is only 0.5 mols of O2:
state = ChemicalState(system)
state.set("O2", 0.5, "mol")
Next, we create an EquilibriumSpecs object in which we specify that temperature, volume, and pressure are constrained properties. We also specify that the system is open to CH4. This is all followed by the creation of a specialized and optimized EquilibriumSolver object to handle this kind of equilibrium calculations:
specs = EquilibriumSpecs(system)
specs.temperature()
specs.volume()
specs.pressure()
specs.openTo("CH4")
solver = EquilibriumSolver(specs)
We can now create the EquilibriumConditions object in which we provide the values for temperature, volume, and pressure, and then perform the calculation:
conditions = EquilibriumConditions(specs)
conditions.temperature(1000.0, "celsius")
conditions.volume(10.0, "cm3")
conditions.pressure(1.0, "GPa")
solver.solve(state, conditions)
print(state.props())
+----------------------------------------+--------------+-----------+
| Property | Value | Unit |
+----------------------------------------+--------------+-----------+
| Temperature | 1273.1500 | K |
| Pressure | 10000.0000 | bar |
| Volume | 1.0000e-05 | m3 |
| Gibbs Energy | -488420.1304 | J |
| Enthalpy | -184792.5765 | J |
| Entropy | 238.4853 | J/K |
| Internal Energy | -194792.5765 | J |
| Helmholtz Energy | -498420.1304 | J |
| Charge | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 1.5995e+00 | mol |
| :: C | 3.9988e-01 | mol |
| :: O | 1.0000e+00 | mol |
| Species Amount: | | |
| :: CH4 | 1.2748e-01 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: CO2 | 2.3311e-01 | mol |
| :: CO | 3.9293e-02 | mol |
| :: H2O | 4.9449e-01 | mol |
| :: H2 | 5.0306e-02 | mol |
| Mole Fraction: | | |
| :: CH4 | 1.3495e-01 | mol/mol |
| :: O2 | 1.0586e-16 | mol/mol |
| :: CO2 | 2.4676e-01 | mol/mol |
| :: CO | 4.1594e-02 | mol/mol |
| :: H2O | 5.2345e-01 | mol/mol |
| :: H2 | 5.3252e-02 | mol/mol |
| Activity Coefficient: | | |
| :: CH4 | 1.0000 | - |
| :: O2 | 1.0000 | - |
| :: CO2 | 1.0000 | - |
| :: CO | 1.0000 | - |
| :: H2O | 1.0000 | - |
| :: H2 | 1.0000 | - |
| Activity: | | |
| :: CH4 | 1.3495e+03 | - |
| :: O2 | 1.0586e-12 | - |
| :: CO2 | 2.4676e+03 | - |
| :: CO | 4.1594e+02 | - |
| :: H2O | 5.2345e+03 | - |
| :: H2 | 5.3252e+02 | - |
| lg(Activity): | | |
| :: CH4 | 3.1302 | - |
| :: O2 | -11.9753 | - |
| :: CO2 | 3.3923 | - |
| :: CO | 2.6190 | - |
| :: H2O | 3.7189 | - |
| :: H2 | 2.7263 | - |
| ln(Activity): | | |
| :: CH4 | 7.2075 | - |
| :: O2 | -27.5741 | - |
| :: CO2 | 7.8110 | - |
| :: CO | 6.0305 | - |
| :: H2O | 8.5630 | - |
| :: H2 | 6.2776 | - |
| Chemical Potential: | | |
| :: CH4 | -2.7843e+05 | J/mol |
| :: O2 | -5.8051e+05 | J/mol |
| :: CO2 | -6.2217e+05 | J/mol |
| :: CO | -3.2477e+05 | J/mol |
| :: H2O | -4.2294e+05 | J/mol |
| :: H2 | -1.2553e+05 | J/mol |
| Standard Volume: | | |
| :: CH4 | 0.0000e+00 | m3/mol |
| :: O2 | 0.0000e+00 | m3/mol |
| :: CO2 | 0.0000e+00 | m3/mol |
| :: CO | 0.0000e+00 | m3/mol |
| :: H2O | 0.0000e+00 | m3/mol |
| :: H2 | 0.0000e+00 | m3/mol |
| Standard Gibbs Energy (formation): | | |
| :: CH4 | -354726.7026 | J/mol |
| :: O2 | -288624.1836 | J/mol |
| :: CO2 | -704857.1814 | J/mol |
| :: CO | -388605.1648 | J/mol |
| :: H2O | -513583.4128 | J/mol |
| :: H2 | -191986.0517 | J/mol |
| Standard Enthalpy (formation): | | |
| :: CH4 | -14294.0889 | J/mol |
| :: O2 | 32390.4240 | J/mol |
| :: CO2 | -344886.6903 | J/mol |
| :: CO | -79597.1200 | J/mol |
| :: H2O | -204066.9932 | J/mol |
| :: H2 | 29072.1299 | J/mol |
| Standard Entropy (formation): | | |
| :: CH4 | 267.3940 | J/(mol*K) |
| :: O2 | 252.1420 | J/(mol*K) |
| :: CO2 | 282.7400 | J/(mol*K) |
| :: CO | 242.7114 | J/(mol*K) |
| :: H2O | 243.1107 | J/(mol*K) |
| :: H2 | 173.6309 | J/(mol*K) |
| Standard Internal Energy (formation): | | |
| :: CH4 | -14294.0889 | J/mol |
| :: O2 | 32390.4240 | J/mol |
| :: CO2 | -344886.6903 | J/mol |
| :: CO | -79597.1200 | J/mol |
| :: H2O | -204066.9932 | J/mol |
| :: H2 | 29072.1299 | J/mol |
| Standard Helmholtz Energy (formation): | | |
| :: CH4 | -354726.7026 | J/mol |
| :: O2 | -288624.1836 | J/mol |
| :: CO2 | -704857.1814 | J/mol |
| :: CO | -388605.1648 | J/mol |
| :: H2O | -513583.4128 | J/mol |
| :: H2 | -191986.0517 | J/mol |
| Standard Heat Capacity (constant P): | | |
| :: CH4 | 84.1523 | J/(mol*K) |
| :: O2 | 35.9257 | J/(mol*K) |
| :: CO2 | 56.9320 | J/(mol*K) |
| :: CO | 34.4671 | J/(mol*K) |
| :: H2O | 44.7425 | J/(mol*K) |
| :: H2 | 31.3023 | J/(mol*K) |
| Standard Heat Capacity (constant V): | | |
| :: CH4 | 84.1523 | J/(mol*K) |
| :: O2 | 35.9257 | J/(mol*K) |
| :: CO2 | 56.9320 | J/(mol*K) |
| :: CO | 34.4671 | J/(mol*K) |
| :: H2O | 44.7425 | J/(mol*K) |
| :: H2 | 31.3023 | J/(mol*K) |
+----------------------------------------+--------------+-----------+
The above table shows that we obtained an equilibrium state with desired values for temperature, pressure, and volume (in SI units). We can also see below the amount of CH4 that entered the system so that all these prescribed conditions could be attained:
print(f"Amount of CH4 titrated in: {state.equilibrium().titrantAmount('CH4')} mol")
Amount of CH4 titrated in: 0.3998837695333316 mol
Recap#
Reaktoro can solve chemical equilibrium problems with a variety of specifications. This guide demonstrated how this can be done for some relatively simple cases using classes:
We will later learn how these classes can be used to model more complicated equilibrium problems.