Chemical equilibrium: the basics#

Written by Allan Leal (ETH Zurich) on Jan 4th, 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!

Reaktoro performs chemical equilibrium calculations by minimizing the Gibbs energy of the system [Leal et al., 2017]. In a Gibbs energy minimization problem, temperature and pressure are prescribed and the chemical system is closed. Thus, without mass transfer in or out of the system, the initial amounts of chemical elements and electric charge are conserved in the reactive process. However, there are many different classes of equilibrium problems where:

  • temperature and/or pressure may be unknown (e.g., combustion in a rigid and adiabatic chamber where volume and internal energy are specified and temperature and pressure are calculated);

  • the chemical system is open to certain substances (e.g., aqueous solution in equilibrium with a mixture of gases with known partial pressures or fugacities, such as the atmosphere).

Reaktoro’s chemical equilibrium algorithm supports all these cases because it implements an algorithm to solve parametric Gibbs energy minimization problems with general equilibrium constraints. Thus, Reaktoro is not restricted to closed systems and prescribed temperature and pressure conditions, as we will see in the following sections and tutorials.

We give a simple example below in which we demonstrate the most basic form of performing chemical equilibrium calculations in Reaktoro. In this example, we want to model the combustion of 1 mol of CH4 and 0.5 mol of O2 at 1000 ºC and 100 bar . First, we define a chemical system consisting of only a gaseous phase:

from reaktoro import *

db = NasaDatabase("nasa-cea")

gases = GaseousPhase("CH4 O2 CO2 CO H2O H2")

system = ChemicalSystem(db, gases)

We then create a chemical state for this system representing the initial temperature, pressure, and composition conditions stated above:

state = ChemicalState(system)
state.temperature(1000, "celsius")
state.pressure(100, "bar")
state.set("CH4", 1.0, "mol")
state.set("O2",  0.5, "mol")

print("=== INITIAL STATE ===")
print(state)

=== INITIAL STATE ===
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |  1273.1500 |    K |
| Pressure        |   100.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 |
+-----------------+------------+------+

And we finally create an EquilibriumSolver object to perform the desired chemical equilibrium calculation:

solver = EquilibriumSolver(system)
solver.solve(state)  # equilibrate the `state` object!

print("=== FINAL STATE ===")
print(state)

=== FINAL STATE ===
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |  1273.1500 |    K |
| Pressure        |   100.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          | 3.7606e-01 |  mol |
| :: O2           | 1.0000e-16 |  mol |
| :: CO2          | 9.3630e-02 |  mol |
| :: CO           | 5.3031e-01 |  mol |
| :: H2O          | 2.8243e-01 |  mol |
| :: H2           | 9.6544e-01 |  mol |
+-----------------+------------+------+

Note that the initial and final amounts of the elements H, C, and O are identical (the same applies to the electric charge of the system, even though there are no charged species in this specific problem).

Note

If you’ve read the previous guides, you’ve probably already encountered the use of equilibrate(state) to bring state, an object of class ChemicalState, to a state of chemical equilibrium. This method exists for convenience and performs a Gibbs energy minimization on the entire chemical system. It is similar to using EquilibriumSolver as follows:

solver = EquilibriumSolver(system)
solver.solve(state)

Tip

If you require many equilibrium calculations, particularly in the context of reactive transport simulations, it’s worth considering the use of a persistent EquilibriumSolver object instead of making repeated equilibrate calls, each creating a new EquilibriumSolver object, which may lead to unnecessary computation overhead.

Checking if the calculation was successful#

Performing chemical equilibrium calculations is not a trivial task. It involves highly complicated algorithms to solve different mathematical problems (e.g., non-linear equations, matrix equations). Therefore, it is possible that a calculation could fail. This can happen due to many factors. However, the most common are:

  • the formulation of the equilibrium problem is ill-formed (e.g. specifying a set of conflicting equilibrium constraints that cannot be achieved chemically and/or thermodynamically)

  • the use of poor initial guesses.

So you might want to check at the end of each equilibrium calculation if it was successful:

state = ChemicalState(system)
state.temperature(1000, "celsius")
state.pressure(100, "bar")
state.set("CH4", 1.0, "mol")
state.set("O2",  0.5, "mol")

solver = EquilibriumSolver(system)
result = solver.solve(state)

print("Succeesful?", result.succeeded())

Succeesful? True

Continue reading to learn more advanced ways to perform chemical equilibrium calculations with Reaktoro!