Chemical equilibrium with custom constraints
Contents
Chemical equilibrium with custom constraints#
Written by Allan Leal (ETH Zurich) on Jan 21th, 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 already implements many different types of chemical equilibrium constraints. The EquilibriumSpecs class allows you to constrain the following properties in an equilibrium state:
the temperature of the system
the pressure of the system
the volume of the system
the internal energy of the system
the enthalpy of the system
the Gibbs energy of the system
the Helmholtz energy of the system
the entropy of the system
the electric charge of the system
the amount of an element in the system
the amount of an element in a phase
the mass of an element in the system
the mass of an element in a phase
the amount of a phase
the mass of a phase
the volume of a phase
the electric charge of a phase
the chemical potential of a species
the activity of a species
the fugacity of a gas
the pH of an aqueous phase
the pMg of an aqueous phase
the pE of an aqueous phase
the Eh of an aqueous phase
Tip
This list of properties that can be constrained may be more extensive now! For a more up-to-date list, check the API of class EquilibriumSpecs.
These options may not suffice to you. This tutorial will demonstrate how your own constraints can be defined and used by the EquilibriumSolver class to honor the conditions you want to see attained at an equilibrium state.
Before we proceed, we need to understand how equilibrium constraints are formulated in Reaktoro and used by its equilibrium solver. Let’s consider you want to impose a value for the enthalpy of the system, \(H_{\mathrm{desired}}\). A chemical state that honors all imposed equilibrium constraints and mass/charge conservation conditions is sought during the chemical equilibrium calculation. This happens iteratively, and in a given iteration of the algorithm, the enthalpy of the system, \(H_{\mathrm{current}}\), may still be different than that you want, \(H_{\mathrm{desired}}\). At some point, the algorithm will converge to a chemical state in which the equilibrium constraint residual satisfies:
where \(\epsilon_{\mathrm{tolerance}}\) is a small positive tolerance value (e.g., 10-6).
We will learn in the next sections that equilibrium constraints in Reaktoro can be defined by writing residual expressions like the one above.
Tip
Reaktoro is not restricted to impose only the value of a property at equilibrium, such as shown in the previous example in which enthalpy is constrained. Because equilibrium constraints are formulated in terms of residual expressions, you can write whatever you want there, even expressions involving multiple properties!
Creating your own volume and internal energy constraints#
In a previous example, we modeled the combustion of CH4 in a rigid and adiabatic chamber. We’ll repeat this example here. However, this time, we’ll define the volume and internal energy constraints from scratch. This example should be helpful for you to understand the fundamentals of defining functional constraints and attaching them to the equilibrium solver.
So, let’s start as usual with the definition of our chemical system and the creation of an initial chemical state for it (this will be identical to what we did here!).
from reaktoro import *
db = NasaDatabase("nasa-cea")
gases = GaseousPhase("CH4 O2 CO2 CO H2O H2")
system = ChemicalSystem(db, gases)
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")
Let’s calculate the chemical properties of this system in this initial state and save the initial volume and internal energy in the variables V0 and U0 for later use:
props0 = ChemicalProps(state0)
V0 = props0.volume()
U0 = props0.internalEnergy()
Let’s now finally define our own constraints (from scratch!). Below is the complete sequence of steps with just a few comments. More details will be given later.
specs = EquilibriumSpecs(system)
idxV = specs.addInput("V") # add "V" as the symbol for a new input condition to the equilibrium problem
idxU = specs.addInput("U") # add "U" as the symbol for a new input condition to the equilibrium problem
volumeConstraint = EquationConstraint()
volumeConstraint.id = "VolumeConstraint" # give some identification name to the constraint (it's up to you how you call this)
volumeConstraint.fn = lambda props, w: props.volume() - w[idxV] # the residual function defining V(current) - V(desired)
internalEnergyConstraint = EquationConstraint()
internalEnergyConstraint.id = "InternalEnergyConstraint" # give some identification name to the constraint (it's up to you how you call this)
internalEnergyConstraint.fn = lambda props, w: props.internalEnergy() - w[idxU] # the residual function defining U(current) - U(desired)
specs.addConstraint(volumeConstraint)
specs.addConstraint(internalEnergyConstraint)
solver = EquilibriumSolver(specs)
As usual, we start with the creation of an EquilibriumSpecs object. Besides the code comments above, it’s worth commenting further that:
idxVandidxUare the indices of the newly added inputs to the equilibrium problem (namelyVandU)lambda props, w: props.volume() - w[idxV]is a lambda function defining our residual expression for the volume constraint: \(V_\mathrm{current} - V_\mathrm{desired}\)lambda props, w: props.internalEnergy() - w[idxV]is a lambda function defining our residual expression for the internal energy constraintpropsis theChemicalPropsobject containing the current chemical properties of the system from where we get current volume and internal energy: \(U_\mathrm{current} - U_\mathrm{desired}\)wis the array with the values of the inputs introduced in the equilibrium problem from where we get the desired volume and internal energyw[idxV]andw[idxU]are the values of introduced inputsVandU(the desired values of volume and internal energy)
It’s now time to create our EquilibriumConditions object, where specify the values for the inputs V and U:
conditions = EquilibriumConditions(specs)
conditions.set("V", V0) # remember the symbol "V" introduced before? you're using it here!
conditions.set("U", U0) # remember the symbol "U" introduced before? you're using it here!
conditions.setLowerBoundTemperature(25.0, "celsius")
conditions.setUpperBoundTemperature(2000.0, "celsius")
conditions.setLowerBoundPressure(1.0, "bar")
conditions.setUpperBoundPressure(1000.0, "bar")
After this, we are ready to solve the chemical equilibrium problem that models the combustion of CH4 in a rigid and adiabatic chamber:
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 |
+-----------------+------------+------+
That’s it! You may want to check if this computed equilibrium state is identical to that found in this tutorial, where we used predefined volume and internal energy constraints to solve the same problem.
Creating your own pH constraint#
We’ve seen in this previous tutorial how to use the predefined equilibrium constraint for pH in class EquilibriumSpecs. We also saw that using that pH constraint causes the system to be open to H+. However, what if we want to fix the pH of an aqueous solution by titrating it with another substance?
We consider thus the following case. At 30 °C and 5 atm, how much CO2 must be titrated into an aqueous solution saturated with mineral calcite (CaCO3) to obtain pH 7? This problem can be solved by formulating a chemical equilibrium calculation in which temperature, pressure, and pH are constrained and that the system is open to the mass transfer of CO2.
To solve this problem, we start with the definition of a chemical system with the following phases:
an aqueous phase (with aqueous species that can be formed from elements H, O, Na, Cl, C, and Ca);
a gaseous phase (with only CO2(g) and H2O(g) species)
a solid phase (representing the mineral calcite)
We also set some activity models for the aqueous and gaseous phases (check this tutorial for more information about specifying activity models):
from reaktoro import *
db = PhreeqcDatabase("phreeqc.dat")
solution = AqueousPhase(speciate("H O Na Cl C Ca"))
solution.set(ActivityModelPhreeqc(db))
gases = GaseousPhase("CO2(g) H2O(g)")
gases.set(ActivityModelPengRobinsonPhreeqcOriginal())
calcite = MineralPhase("Calcite")
system = ChemicalSystem(db, solution, gases, calcite)
Let’s now define a helper function that computes pH for given chemical properties of a system:
def pH(props: ChemicalProps):
return -props.speciesActivityLg("H+") # this results in -log10(a[H+]), which is pH
We’ll use this convenience function next (in the lambda function defining the residual expression), while creating our EquilibriumSpecs object in which we provide the specifications for the equilibrium problem we want to formulate.
Tip
The method call props.speciesActivityLg("H+") will be performed several times. If you would like to avoid an index search for H+ in every such call, you could do the following instead:
idxH = system.species().index("H+")
def pH(props: ChemicalProps):
return -props.speciesActivityLg(idxH)
We want to perform a chemical equilibrium calculation in which the following properties are constrained:
temperature;
pressure; and
pH.
Moreover, our system is open to CO2. These specifications are all provided next in the EquilibriumSpecs object we create:
specs = EquilibriumSpecs(system)
specs.addInput("pH") # add "pH" as the symbol for a new input condition to the equilibrium problem
pHConstraint = EquationConstraint()
pHConstraint.id = "pHConstraint" # give some identification name to the constraint (it's up to you how you call this)
pHConstraint.fn = lambda props, w: pH(props) - w[0] # the residual function defining pH(current) - pH(desired)
specs.temperature()
specs.pressure()
specs.addConstraint(pHConstraint)
specs.openTo("CO2")
solver = EquilibriumSolver(specs)
We’ll now create an initial chemical state for the system representing a 1 molal NaCl aqueous solution mixed with a sufficient amount of calcite to saturate the fluid at 25 °C and 1 atm.
state = ChemicalState(system)
state.temperature(25.0, "celsius")
state.pressure(1.0, "atm")
state.set("H2O", 1.0, "kg")
state.set("Na+", 1.0, "mol")
state.set("Cl-", 1.0, "mol")
state.set("Calcite", 10.0, "mol") # plenty of Calcite to ensure saturation levels!
Finally, we set the conditions of this system at the chemical equilibrium state of interest (using EquilibriumConditions) and then use EquilibriumSolver to compute this state of equilibrium:
conditions = EquilibriumConditions(specs)
conditions.temperature(30.0, "celsius")
conditions.pressure(5.0, "atm")
conditions.set("pH", 7.0) # remember the symbol "pH" introduced before? you're using it here!
solver = EquilibriumSolver(specs)
solver.solve(state, conditions)
aprops = AqueousProps(state)
print(aprops)
+---------------------------------+-------------+-------+
| Property | Value | Unit |
+---------------------------------+-------------+-------+
| Temperature | 303.1500 | K |
| Pressure | 5.0663 | bar |
| Ionic Strength (Effective) | 1.0127 | molal |
| Ionic Strength (Stoichiometric) | 1.0144 | molal |
| pH | 7.0000 | |
| pE | -2.4134 | |
| Eh | -0.1452 | V |
| Charge Molality | -1.1091e-16 | molal |
| Element Molality: | | |
| :: C | 1.0525e-02 | molal |
| :: Na | 1.0001e+00 | molal |
| :: Cl | 1.0001e+00 | molal |
| :: Ca | 4.8575e-03 | molal |
| Species Molality: | | |
| :: CO3-2 | 1.4603e-05 | molal |
| :: H+ | 1.3489e-07 | molal |
| :: CO2 | 8.8708e-04 | molal |
| :: (CO2)2 | 2.1409e-08 | molal |
| :: HCO3- | 7.8249e-03 | molal |
| :: CH4 | 1.0001e-16 | molal |
| :: Ca+2 | 4.7329e-03 | molal |
| :: CaCO3 | 4.5248e-06 | molal |
| :: CaHCO3+ | 1.2005e-04 | molal |
| :: CaOH+ | 2.3414e-09 | molal |
| :: Cl- | 1.0001e+00 | molal |
| :: H2 | 1.0001e-16 | molal |
| :: Na+ | 9.9841e-01 | molal |
| :: NaCO3- | 5.8000e-05 | molal |
| :: NaHCO3 | 1.6156e-03 | molal |
| :: OH- | 2.4708e-07 | molal |
| :: NaOH | 1.0001e-16 | molal |
| :: O2 | 1.0001e-16 | molal |
| Saturation Indices: | | |
| :: CH4(g) | -13.7565 | - |
| :: CO2(g) | -1.8192 | - |
| :: H2(g) | -9.9372 | - |
| :: H2O(g) | -1.7645 | - |
| :: O2(g) | -63.9100 | - |
| :: Aragonite :: CaCO3 | -0.1495 | - |
| :: Calcite :: CaCO3 | 0.0000 | - |
| :: Halite :: NaCl | -1.9301 | - |
| Saturation Ratios: | | |
| :: CH4(g) | 1.7518e-14 | - |
| :: CO2(g) | 1.5162e-02 | - |
| :: H2(g) | 1.1556e-10 | - |
| :: H2O(g) | 1.7198e-02 | - |
| :: O2(g) | 1.2304e-64 | - |
| :: Aragonite :: CaCO3 | 7.0878e-01 | - |
| :: Calcite :: CaCO3 | 1.0000e+00 | - |
| :: Halite :: NaCl | 1.1745e-02 | - |
+---------------------------------+-------------+-------+
The aqueous properties displayed above should show that the prescribed temperature, pressure, and pH values are satisfied.
You may be wondering how much CO2 was transferred to the system for this pH value to be attained. We show this next:
tCO2 = state.equilibrium().titrantAmount("CO2") # the amount of CO2 titrated into the system
print(f"Amount of CO2 titrated into the system was {tCO2 * 1e+3} mmol")
Amount of CO2 titrated into the system was 5.666852117092365 mmol
That’s it. You should now be able to adapt these tutorials and define the constraints you want (if not already covered in the EquilibriumSpecs class).