Chemical equilibrium with fixed pH and charge balance
Chemical equilibrium with fixed pH and charge balance#
Written by Allan Leal (ETH Zurich) on Jan 7th, 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!
As we saw in the previous tutorial, imposing the pH of a solution in a chemical equilibrium calculation leaves the system open to H+. As a result, the final solution may not have zero electrical charge.
So, in this tutorial we demonstrate how to perform chemical equilibrium calculations with fixed pH values while simultaneously satisfying charge balance.
Let’s create a chemical system (using ChemicalSystem) with just a single aqueous phase containing a few selected species. For this phase, we use a Debye–Hückel activity model. The required thermodynamic data for the species is extracted from the SupcrtDatabase object:
from reaktoro import *
db = SupcrtDatabase("supcrtbl")
solution = AqueousPhase("H2O(aq) H+ OH- Na+ Cl- HCO3- CO2(aq) CO3-2")
solution.set(ActivityModelDebyeHuckel())
system = ChemicalSystem(db, solution)
We want to perform a chemical equilibrium calculation in which the following properties are constrained:
temperature,
pressure,
pH, and
electric charge.
To attain a desired charge value in the solution, we will need to make the system open to a anion (negative ion) to counterbalance the amount of H+ entering/leaving the system. We will use Cl- for this.
Let’s now create a specialized and optimized EquilibriumSolver object to handle equilibrium problems with these specifications (note the use of EquilibriumSpecs below to achieve this):
specs = EquilibriumSpecs(system)
specs.temperature()
specs.pressure()
specs.pH()
specs.charge()
specs.openTo("Cl-")
solver = EquilibriumSolver(specs)
We should now create an initial chemical state for our chemical system:
state = ChemicalState(system)
state.temperature(30.0, "celsius")
state.pressure(1.0, "bar")
state.set("H2O(aq)", 1.00, "kg")
state.set("Na+" , 0.01, "mol")
state.set("Cl-" , 0.01, "mol")
state.set("CO2(aq)", 0.10, "mol")
print("INITIAL STATE")
print(state)
INITIAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 303.1500 | K |
| Pressure | 1.0000 | bar |
| Charge: | -3.0000e-16 | mol |
| Element Amount: | | |
| :: H | 1.1102e+02 | mol |
| :: C | 1.0000e-01 | mol |
| :: O | 5.5708e+01 | mol |
| :: Na | 1.0000e-02 | mol |
| :: Cl | 1.0000e-02 | mol |
| Species Amount: | | |
| :: H2O(aq) | 5.5508e+01 | mol |
| :: H+ | 1.0000e-16 | mol |
| :: OH- | 1.0000e-16 | mol |
| :: Na+ | 1.0000e-02 | mol |
| :: Cl- | 1.0000e-02 | mol |
| :: HCO3- | 1.0000e-16 | mol |
| :: CO2(aq) | 1.0000e-01 | mol |
| :: CO3-2 | 1.0000e-16 | mol |
+-----------------+-------------+------+
Let’s now define the temperature, pressure, pH, and charge conditions we want to impose at equilibrium for this aqueous solution currently in disequilibrium. We do this by using an object of class EquilibriumConditions:
conditions = EquilibriumConditions(specs)
conditions.temperature(50.0, "celsius")
conditions.pressure(10.0, "bar")
conditions.pH(2.0)
conditions.charge(0.0)
We have everything we need now to perform the equilibrium calculation:
an equilibrium solver (the
solverobject of classEquilibriumSolver),an initial chemical state to equilibrate (the
stateobject of classChemicalState), andthe conditions we are imposing for this equilibrium state (the
conditionsobject of classEquilibriumConditions).
The code below will perform the equilibrium calculation and modify the state object to reflect the computed chemical equilibrium state:
result = solver.solve(state, conditions)
print("FINAL STATE")
print(state)
FINAL STATE
+-----------------+-------------+------+
| Property | Value | Unit |
+-----------------+-------------+------+
| Temperature | 323.1500 | K |
| Pressure | 10.0000 | bar |
| Charge: | -2.4653e-18 | mol |
| Element Amount: | | |
| :: H | 1.1103e+02 | mol |
| :: C | 1.0000e-01 | mol |
| :: O | 5.5708e+01 | mol |
| :: Na | 1.0000e-02 | mol |
| :: Cl | 2.1127e-02 | mol |
| Species Amount: | | |
| :: H2O(aq) | 5.5508e+01 | mol |
| :: H+ | 1.1133e-02 | mol |
| :: OH- | 6.1358e-12 | mol |
| :: Na+ | 1.0000e-02 | mol |
| :: Cl- | 2.1127e-02 | mol |
| :: HCO3- | 6.1507e-06 | mol |
| :: CO2(aq) | 9.9994e-02 | mol |
| :: CO3-2 | 6.0702e-14 | mol |
+-----------------+-------------+------+
It’s always advisable to verify if the calculation succeeded:
print("Successful computation!" if result.succeeded() else "Computation has failed!")
Successful computation!
We can see from the table above that electric charge in the computed equilibrium state is numerically zero (i.e., very close to machine precision, 2.220446e-16).
Let’s check the aqueous properties of this equilibrium state to verify if we obtained the desired pH:
aprops = AqueousProps(state)
print("AQUEOUS PROPERTIES AT EQUILIBRIUM")
print(aprops)
AQUEOUS PROPERTIES AT EQUILIBRIUM
+---------------------------------+-------------+-------+
| Property | Value | Unit |
+---------------------------------+-------------+-------+
| Temperature | 323.1500 | K |
| Pressure | 10.0000 | bar |
| Ionic Strength (Effective) | 0.0211 | molal |
| Ionic Strength (Stoichiometric) | 0.0211 | molal |
| pH | 2.0000 | |
| pE | 1.3421 | |
| Eh | 0.0861 | V |
| Charge Molality | -2.4653e-18 | molal |
| Element Molality: | | |
| :: C | 1.0000e-01 | molal |
| :: Na | 1.0000e-02 | molal |
| :: Cl | 2.1127e-02 | molal |
| Species Molality: | | |
| :: H+ | 1.1133e-02 | molal |
| :: OH- | 6.1358e-12 | molal |
| :: Na+ | 1.0000e-02 | molal |
| :: Cl- | 2.1127e-02 | molal |
| :: HCO3- | 6.1507e-06 | molal |
| :: CO2(aq) | 9.9994e-02 | molal |
| :: CO3-2 | 6.0702e-14 | molal |
| Saturation Indices: | | |
| :: CH4(g) | -5.4132 | - |
| :: CO(g) | -9.9715 | - |
| :: CO2(g) | -0.2822 | - |
| :: Ethylene(g) :: C2H4 | -24.4885 | - |
| :: H2(g) | -7.1497 | - |
| :: H2O(g) | -1.9187 | - |
| :: O2(g) | -64.0586 | - |
| :: Phenol(g) :: C6H5OH | -46.3912 | - |
| :: m-Cresol(g) :: CH3C6H4(OH) | -51.6682 | - |
| :: o-Cresol(g) :: CH3C6H4(OH) | -52.2246 | - |
| :: p-Cresol(g) :: CH3C6H4(OH) | -53.2650 | - |
| :: Diamond :: C | -0.4748 | - |
| :: Graphite :: C | 0.0238 | - |
| :: Halite :: NaCl | -5.3824 | - |
| Saturation Ratios: | | |
| :: CH4(g) | 3.8617e-06 | - |
| :: CO(g) | 1.0677e-10 | - |
| :: CO2(g) | 5.2215e-01 | - |
| :: Ethylene(g) :: C2H4 | 3.2474e-25 | - |
| :: H2(g) | 7.0845e-08 | - |
| :: H2O(g) | 1.2058e-02 | - |
| :: O2(g) | 8.7381e-65 | - |
| :: Phenol(g) :: C6H5OH | 4.0622e-47 | - |
| :: m-Cresol(g) :: CH3C6H4(OH) | 2.1470e-52 | - |
| :: o-Cresol(g) :: CH3C6H4(OH) | 5.9621e-53 | - |
| :: p-Cresol(g) :: CH3C6H4(OH) | 5.4329e-54 | - |
| :: Diamond :: C | 3.3514e-01 | - |
| :: Graphite :: C | 1.0563e+00 | - |
| :: Halite :: NaCl | 4.1454e-06 | - |
+---------------------------------+-------------+-------+
Great! We can now see that pH is exactly 2, as we enforced, and temperature and pressure are 50 °C and 10 bar (but displayed in the table as 323.15 K and 106 Pa respectively).
Let’s now check how much H+ and Cl- entered/leaved the system (compared to the initial state) so that both pH and charge conditions could be reached:
print(f"Amount of H+ that entered the system: {state.equilibrium().implicitTitrantAmount('H+')} mol")
print(f"Amount of Cl- that entered the system: {state.equilibrium().explicitTitrantAmount('Cl-')} mol")
Amount of H+ that entered the system: 0.011127042666330142 mol
Amount of Cl- that entered the system: 0.011127042666329844 mol
We can also use the convenient titrantAmount method which will search for the titrant name among both implicit and explicit titrants:
print(f"Amount of H+ that entered the system: {state.equilibrium().titrantAmount('H+')} mol")
print(f"Amount of Cl- that entered the system: {state.equilibrium().titrantAmount('Cl-')} mol")
Amount of H+ that entered the system: 0.011127042666330142 mol
Amount of Cl- that entered the system: 0.011127042666329844 mol
That’s it! The results shown above demonstrate that our constraints were satisfied at equilibrium and that both H+ and Cl- entered the system with a numerically equal amount.