This tutorial demonstrates how to use Reaktoro to perform a chemical equilibrium calculation with carbon species.
We start by importing the reaktoro
package:
from reaktoro import *
The default thermodynamic databases embedded into Reaktoro is SUPCRT92, so you do not have to initialize the
database db = Database('supcrt98.xml')
, unless you use an alternative one.
Class ChemicalEditor
provides convenient operations to initialize
ChemicalSystem and
ReactionSystem instances.
editor = ChemicalEditor()
Alternatively, the editor can initialize from the file:
# Initialize a thermodynamic database with supcrt98.xml
db = Database("supcrt98.xml")
# Define the editor of the chemical system
editor = ChemicalEditor(db)
To define a chemical system, aqueous, gaseous, and mineral phases must be added. For the aqueous phase, it can be done from a list of chemical element names. The database will be searched for all species that could be formed out of those elements:
editor.addAqueousPhaseWithElements("H O C Ca Cl Mg")
To add gaseous phase, we provide a specific list of gaseous species that must be used:
editor.addGaseousPhase(["H2O(g)", "CO2(g)", "H2(g)", "O2(g)", "CH4(g)"])
For the MineralPhase object, we add two pure mineral phases:
editor.addMineralPhase("Calcite")
editor.addMineralPhase("Dolomite")
To create an object of class ChemicalSystem using the chemical system definition details stored in the object editor, we use
system = ChemicalSystem(editor)
To output the details of the chemical system to the console, i.e., the list of elements, species, and different phases, we can execute
print(system)
For the specific information about the chemical system, one can use the following methods:
print("Number of chemical species: ", len(system.species()))
print("Number of phases: ", len(system.phases()))
print("Number of elements: ", len(system.elements()))
print("List of species in chemical system: \n")
for species in system.species():
print(species.name())
See tutorials [ChemicalEditor](cl.chemical-editor.ipynb) and [ChemicalSystem](cl.chemical-system.ipynb) for more detailed explanation of capabilities of these classes.
Next, we create an equilibrium problem with our prescribed equilibrium conditions for amounts of elements that are consistent with our intention of calculating reaction of calcite with injected 0.002 molal MaCl2 brine. Both temperature and pressure are assumed to be default values, i.e., 25 °C and 1 bar, respectively.
problem = EquilibriumProblem(system)
problem.add("H2O", 1, "kg")
problem.add("MgCl2", 0.002, "mol")
problem.add("CaCO3", 1, "mol")
In this step, we use the equilibrate function to calculate the chemical equilibrium state of the system with the given equilibrium conditions stored in the object problem.
state = equilibrate(problem)
Reaktoro uses an efficient Gibbs energy minimization computation to determine the species amounts that correspond
to a state of minimum Gibbs energy in the system while satisfying the prescribed amount conditions for the
temperature, pressure, and element amounts. The result is stored in the object state
of class
ChemicalState.
To output the result of equilibration to the console, we use
print(state)
We will obtain the table describing the chemical state of the system. For example, the molar amounts, molar fractions, activities, activity coefficients, and chemical potentials of the species. The molar volumes of the phases, the amounts of each element in the phases, and also the total phase molar amounts.
Alternatively, to save equilibrated state into a file, one can use method ChemicalState::output:
state.output("result.txt") # Output the equilibrium state into a file result.txt
To print the amounts of some specific aqueous species, one can use
# Print the amounts of some aqueous species
print("Amount of CO2(aq):", state.speciesAmount("CO2(aq)"))
print("Amount of HCO3-:", state.speciesAmount("HCO3-"))
print("Amount of CO3--:", state.speciesAmount("CO3--"))
Similarly, one can also print the amounts of the certain element, say carbon, in both aqueous and gaseous phases
# Print the amounts of element C in both aqueous and gaseous phases
print("Amount of C in aqueous phase:", state.elementAmountInPhase("C", "Aqueous"))
print("Amount of C in gaseous phase:", state.elementAmountInPhase("C", "Gaseous"))
print("Amount of C in calcite phase:", state.elementAmountInPhase("C", "Calcite"))