Chemical equilibrium with given element and charge amounts

Written by Allan Leal (ETH Zurich) on Jan 7th, 2022

So far, you most likely have only seen Reaktoro performing chemical equilibrium calculations with a given initial chemical state in disequilibrium such as the one below:

from reaktoro import *

db = NasaDatabase("nasa-cea")

# Create a chemical system with only a gas phase
system = ChemicalSystem(db, GaseousPhase("CH4 O2 CO2 CO H2O H2"))

# Create an initial state for this chemical system
state = ChemicalState(system)
state.temperature(1000, "celsius")
state.pressure(100, "bar")
state.set("CH4", 1.0, "mol")
state.set("O2",  0.5, "mol")

# Create an equilibrium solver and equilibrates `state`
solver = EquilibriumSolver(system)
solver.solve(state)

print(state)
+-----------------+-----------+------+
| Property        |     Value | Unit |
+-----------------+-----------+------+
| Temperature     |   1273.15 |    K |
| Pressure        |     1e+07 |   Pa |
| Charge:         |         0 |  mol |
| Element Amount: |           |      |
| :: H            |         4 |  mol |
| :: C            |         1 |  mol |
| :: O            |         1 |  mol |
| Species Amount: |           |      |
| :: CH4          |  0.376062 |  mol |
| :: O2           |     1e-16 |  mol |
| :: CO2          | 0.0936296 |  mol |
| :: CO           |  0.530308 |  mol |
| :: H2O          |  0.282433 |  mol |
| :: H2           |  0.965442 |  mol |
+-----------------+-----------+------+

In the background, when you called solver.solve(state) above, Reaktoro used the amounts of species in the initial state to compute the amounts of elements and charge. These are needed in an underlying mass conservation equation solved along with the Gibbs energy minimization problem.

There are applications in which you may have the amounts of elements and charge of the system already and you want to use this as input to the chemical equilibrium calculation.

Note

We use the term components to denote both chemical elements and charge (or any other entities used as building blocks for chemical species).

We demonstrate below how this can be achieved so that we produce exact same results as shown in the table above.

import numpy

Ne = system.elements().size()      # the number of elements in the system

iC = system.elements().index("C")  # the index of component C
iH = system.elements().index("H")  # the index of component H
iO = system.elements().index("O")  # the index of component O
iZ = Ne                            # the index of component charge

b = numpy.zeros(Ne + 1)  # the array with explicit amounts of elements and charge

b[iC] = 1.0  # 1 mol of C
b[iH] = 4.0  # 4 moles of H
b[iO] = 1.0  # 1 mol of O
b[iZ] = 0.0  # 0 mol of charge

state = ChemicalState(system)
state.temperature(1000, "celsius")
state.pressure(100, "bar")

solver.solve(state, b)

print(state)
+-----------------+-----------+------+
| Property        |     Value | Unit |
+-----------------+-----------+------+
| Temperature     |   1273.15 |    K |
| Pressure        |     1e+07 |   Pa |
| Charge:         |         0 |  mol |
| Element Amount: |           |      |
| :: H            |         4 |  mol |
| :: C            |         1 |  mol |
| :: O            |         1 |  mol |
| Species Amount: |           |      |
| :: CH4          |  0.376062 |  mol |
| :: O2           |     1e-16 |  mol |
| :: CO2          | 0.0936296 |  mol |
| :: CO           |  0.530308 |  mol |
| :: H2O          |  0.282433 |  mol |
| :: H2           |  0.965442 |  mol |
+-----------------+-----------+------+

That’s it. The table above should show an identical equilibrium state as computed before.

Tip

If you are handling reactive transport simulations in which the mass conservation equations are formulated in terms of chemical elements and electrical charge, the above procedure should be what you are looking for. At each mesh cell/node (or degree of freedom), you would have a vector b representing the concentrations of elements and charges at that specific point in space. The equilibrium calculation would provide the amounts of species distributed among all phases (e.g., aqueous, liquid, gaseous, mineral, etc.). The associated chemical properties at equilibrium could be used to compute densities, enthalpies, and other properties of the phases.