This tutorial demonstrates how to use Reaktoro to perform a chemical equilibrium calculation with the help of
the EquilibriumSolver class. First, we import
everything from the reaktoro
package by
from reaktoro import *
We start from creating an object of Database class to use in the initialization of the chemical system.
database = Database("supcrt98.xml")
Next, we define which phases and species the chemical system should have. This is done using an instance of ChemicalEditor class.
editor = ChemicalEditor(database)
editor.addAqueousPhaseWithElementsOf("H2O NaCl CO2")
editor.addGaseousPhase(["H2O(g)", "CO2(g)"])
editor.addMineralPhase("Halite")
Here, AqueousPhase is created by specifying the list of compound or substance names, i.e., H2O, NaCl, and CO2, that might not necessarily represent names of species in the database. Function addAqueousPhaseWithElementsOf will brake the list of compounds into a list of element names and the database will be similarly searched for all species that could be formed out of those elements.
Then, GaseousPhase is composed of the names of the provided gaseous species H2O(g) and CO2(g), using function addGaseousPhase. These names must conform to those used in the database that was specified during the initialization of the ChemicalEditor object, otherwise, an exception will be thrown.
The MineralPhase object is created by specifying the names of the mineral species one by one. Analogously to the gaseous species, provided names must coincide with those used in the database (specified during the initialization of ChemicalEditor object) , otherwise, an exception will be thrown. In this case, the method addMineralPhase is used to create single pure mineral phases with halite.
To initialize the chemical system, we use class ChemicalSystem, which requires the instance of ChemicalEditor defined earlier.
system = ChemicalSystem(editor)
The equilibrium problem is described by the class EquilibriumProblem. Here, different properties, such as temperature, pressure, and amounts of compounds, can be provided.
problem = EquilibriumProblem(system)
problem.setTemperature(60, "celsius")
problem.setPressure(300, "bar")
problem.add("H2O", 1, "kg")
problem.add("CO2", 100, "g")
problem.add("NaCl", 0.1, "mol")
In particular, we set the temperature to 60 °C and pressure to 300 bar. To equilibrate the chemical problem, we also add 1kg of water, 100 g of carbon dioxide, and 0.1 mol of sodium chloride.
See tutorials [ChemicalEditor](cl.chemical-editor.ipynb) and [ChemicalSystem](cl.chemical-system.ipynb) for more detailed explanation of capabilities of these classes.
The temperature, pressure, and mole amounts of the elements can be obtained from the definition of equilibrium problem.
T = problem.temperature()
P = problem.pressure()
b = problem.elementAmounts()
Next, we create an object of EquilibriumSolver class that can be reused many times for equilibrium simulations.
solver = EquilibriumSolver(system)
For more detailed overview on the functionality of the class EquilibriumSolver, please check the tutorial EquilibriumSolver.
Similarly, an object of ChemicalState class, which stores the equilibrium state of the system, can be created.
state = ChemicalState(system)
Using method EquilibriumSolver::solve,
the equilibrium state with given (T, P, b) inputs is generated and stored in the object state
.
solver.solve(state, T, P, b)
Here, the object state
serves as the initial guess and the final state of the equilibrium calculation. If known,
the temperature must be provided in Kelvin, and the pressure is expected in Pascal. Vector b
provides the molar
amounts of the elements in the equilibrium partition. Alternatively, one can call this method with the given
equilibrium problem.
To save the calculated chemical equilibrium state into the text-file, we use the method ChemicalState::output.
state.output('state-T-60.txt')
In the saved file, one can note that the amount of halite is of the order $10^{-21}$, which indicates its dissolution in sodium chloride brine.
Calculate the new equilibrium state at temperature increased by 10 °C. For that, we use the previous equilibrium state as the initial guess for improved performance.
solver.solve(state, T + 10.0, P, b)
Save newly calculated chemical equilibrium state with T = 70 °C into the text-file.
state.output('state-T-70.txt')
In comparison to chemical speciation obtained for equilibrium calculation at 60 °C, the ionic strength, reduction potential, and alkalinity slightly decrease, whereas ph slightly increases.