Chemical equilibrium with fixed volume and internal energy
Chemical equilibrium with fixed volume and internal energy#
Written by Allan Leal (ETH Zurich) on Jan 7th, 2022.
Last revised on Jan 21st, 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!
Let’s consider the combustion of CH4 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!).
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.
Let’s create our ChemicalSystem object (one with a single gaseous phase):
from reaktoro import *
db = NasaDatabase("nasa-cea")
gases = GaseousPhase("CH4 O2 CO2 CO H2O H2")
system = ChemicalSystem(db, gases)
Next, we create an initial chemical state for this system in which CH4 and O2 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 CH4 and O2 to 0.75 and 0.25 moles, respectively, followed by scaling the volume of the chemical state to the desired 10 cm3.
state0 = ChemicalState(system)
state0.temperature(25.0, "celsius")
state0.pressure(1.0, "bar")
state0.set("CH4", 0.75, "mol")
state0.set("O2", 0.25, "mol")
state0.scaleVolume(10.0, "cm3")
print("INITIAL STATE")
print(state0)
INITIAL STATE
+-----------------+------------+------+
| Property | Value | Unit |
+-----------------+------------+------+
| Temperature | 298.1500 | K |
| Pressure | 1.0000 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 1.2102e-03 | mol |
| :: C | 3.0255e-04 | mol |
| :: O | 2.0170e-04 | mol |
| Species Amount: | | |
| :: CH4 | 3.0255e-04 | mol |
| :: O2 | 1.0085e-04 | mol |
| :: CO2 | 4.0340e-20 | mol |
| :: CO | 4.0340e-20 | mol |
| :: H2O | 4.0340e-20 | mol |
| :: H2 | 4.0340e-20 | mol |
+-----------------+------------+------+
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.
props0 = ChemicalProps(state0)
V0 = props0.volume() # the initial volume of the gases
U0 = props0.internalEnergy() # the initial internal energy of the gases
print("Initial volume of the gases is", V0, "m3")
print("Initial internal energy of the gases is", U0, "J")
Initial volume of the gases is 1e-05 m3
Initial internal energy of the gases is -23.5698 J
Our chemical equilibrium problem (to model this combustion process) needs to impose the following properties at equilibrium:
volume and
internal energy.
The code below creates an EquilibriumSpecs object to reflect these specifications for the equilibrium constraints, which is then used to create our EquilibriumSolver object.
specs = EquilibriumSpecs(system)
specs.volume()
specs.internalEnergy()
solver = EquilibriumSolver(specs)
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!):
conditions = EquilibriumConditions(specs)
conditions.volume(V0)
conditions.internalEnergy(U0)
conditions.setLowerBoundTemperature(25.0, "celsius")
conditions.setUpperBoundTemperature(2000.0, "celsius")
conditions.setLowerBoundPressure(1.0, "bar")
conditions.setUpperBoundPressure(1000.0, "bar")
Tip
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?!
We have everything we need now to perform the equilibrium calculation, which is done next:
state = ChemicalState(state0) # let's create a copy of the initial state
solver.solve(state, conditions)
print("FINAL STATE")
print(state)
FINAL STATE
+-----------------+------------+------+
| Property | Value | Unit |
+-----------------+------------+------+
| Temperature | 1080.3867 | K |
| Pressure | 5.8362 | bar |
| Charge: | 0.0000e+00 | mol |
| Element Amount: | | |
| :: H | 1.2102e-03 | mol |
| :: C | 3.0255e-04 | mol |
| :: O | 2.0170e-04 | mol |
| Species Amount: | | |
| :: CH4 | 1.2897e-04 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: CO2 | 9.6874e-06 | mol |
| :: CO | 1.6389e-04 | mol |
| :: H2O | 1.8430e-05 | mol |
| :: H2 | 3.2873e-04 | mol |
+-----------------+------------+------+
Note that O2 was entirely consumed in the process. Its amount (1e-16 moles) is not zero because Reaktoro’s equilibrium solver has a default lower bound of 1e-16 moles for the species amounts. Thus, the equilibrium calculation converged with O2 attached to its lower bound.
Tip
You can use EquilibriumOptions to change this default lower bound for species amounts:
options = EquilibriumOptions()
options.epsilon = 1e-30
solver = EquilibriumSolver(specs)
solver.setOptions(options)
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:
V = state.props().volume()
U = state.props().internalEnergy()
print("Volume at initial state:", V0, "m3")
print("Volume at final state:", V, "m3")
print("Internal energy at initial state:", U0, "J")
print("Internal energy at final state:", U, "J")
Volume at initial state: 1e-05 m3
Volume at final state: 1e-05 m3
Internal energy at initial state: -23.5698 J
Internal energy at final state: -23.5698 J
The verification above demonstrates we found an equilibrium state in which volume and internal energy were preserved.
Let’s now check what is the final temperature and pressure when burning CH4 in that rigid and adiabatic chamber with our given initial conditions:
T = units.convert(state.temperature(), "K", "degC") # convert from K to °C
P = units.convert(state.pressure(), "Pa", "bar") # convert from Pa to bar
print("Temperature at final state:", T, "°C")
print("Pressure at final state:", P, "bar")
Temperature at final state: 807.237 °C
Pressure at final state: 5.83621 bar
This concludes this tutorial, which demonstrated how Reaktoro can be used to perform chemical equilibrium calculations with given volume and internal energy.
Keep on reading to learn how to set your own constraints from scratch!