Chemical equilibrium with fixed phase amount#

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!

Consider an aqueous solution at 60 °C obtained by mixing 1 kg of H2O, 1 mol of NaCl and 1 mol of CO2. What must the pressure of this aqueous solution be so that it is saturated with CO2? In other words, what is the pressure that causes a bubble of CO2 to form (the bubble pressure)?

In this tutorial, we formulate a chemical equilibrium problem to calculate the CO2 bubble/saturation pressure of this aqueous saline solution. In this problem, we consider a chemical system with aqueous and gaseous phases. Our equilibrium problem will consider a constraint that imposes the amount of the gas phase to be a minuscule value (to model the limit condition in which a bubble is formed).

Let’s create our ChemicalSystem object with these phases:

from reaktoro import *

db = PhreeqcDatabase("pitzer.dat")

solution = AqueousPhase(speciate("H O Na Cl C"))
solution.set(ActivityModelPitzer())

gases = GaseousPhase("CO2(g) H2O(g)")
gases.set(ActivityModelPengRobinsonPhreeqcOriginal())

system = ChemicalSystem(db, solution, gases)

Next, we create an initial chemical state for this system representing our aqueous solution at 60 °C (with composition 1 molal NaCl and 1 molal CO2):

state = ChemicalState(system)
state.temperature(60.0, "celsius")
state.set("H2O", 1.0, "kg")
state.set("Na+", 1.0, "mol")
state.set("Cl-", 1.0, "mol")
state.set("CO2", 1.0, "mol")

print("INITIAL STATE")
print(state)

INITIAL STATE
+-----------------+-------------+------+
| Property        |       Value | Unit |
+-----------------+-------------+------+
| Temperature     |    333.1500 |    K |
| Pressure        |      1.0000 |  bar |
| Charge:         | -1.0000e-16 |  mol |
| Element Amount: |             |      |
| :: H            |  1.1101e+02 |  mol |
| :: C            |  1.0000e+00 |  mol |
| :: O            |  5.7506e+01 |  mol |
| :: Na           |  1.0000e+00 |  mol |
| :: Cl           |  1.0000e+00 |  mol |
| Species Amount: |             |      |
| :: H+           |  1.0000e-16 |  mol |
| :: H2O          |  5.5506e+01 |  mol |
| :: CO3-2        |  1.0000e-16 |  mol |
| :: CO2          |  1.0000e+00 |  mol |
| :: Cl-          |  1.0000e+00 |  mol |
| :: HCO3-        |  1.0000e-16 |  mol |
| :: Na+          |  1.0000e+00 |  mol |
| :: OH-          |  1.0000e-16 |  mol |
| :: CO2(g)       |  1.0000e-16 |  mol |
| :: H2O(g)       |  1.0000e-16 |  mol |
+-----------------+-------------+------+

We need now to create an EquilibriumSolver that can solve equilibrium problems with the following constrained properties:

  • temperature and

  • amount of gaseous phase.

These constraint specifications are provided to the EquilibriumSpecs object below, which is then used to construct our EquilibriumSolver object:

specs = EquilibriumSpecs(system)
specs.temperature()
specs.phaseAmount("GaseousPhase")

solver = EquilibriumSolver(specs)

The next step is to create an object EquilibriumConditions in which we specify the values of temperature and amount of gas phase (which should be a small amount). We do this in the next code block, which also defines the lower and upper limits for pressure (to avoid negative and unrealistic large pressure values during the equilibrium calculation!):

conditions = EquilibriumConditions(specs)
conditions.temperature(60.0, "celsius")
conditions.phaseAmount("GaseousPhase", 1.0, "umol")  # umol = 1e-6 moles
conditions.setLowerBoundPressure(1.0, "bar")
conditions.setUpperBoundPressure(1000.0, "bar")

We now have everything needed to perform the calculation:

solver.solve(state, conditions)

print("FINAL STATE")
print(state)

FINAL STATE
+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   333.1500 |    K |
| Pressure        |   149.2481 |  bar |
| Charge:         | 2.3903e-17 |  mol |
| Element Amount: |            |      |
| :: H            | 1.1101e+02 |  mol |
| :: C            | 1.0000e+00 |  mol |
| :: O            | 5.7506e+01 |  mol |
| :: Na           | 1.0000e+00 |  mol |
| :: Cl           | 1.0000e+00 |  mol |
| Species Amount: |            |      |
| :: H+           | 1.2655e-03 |  mol |
| :: H2O          | 5.5505e+01 |  mol |
| :: CO3-2        | 5.2102e-10 |  mol |
| :: CO2          | 9.9873e-01 |  mol |
| :: Cl-          | 1.0000e+00 |  mol |
| :: HCO3-        | 1.2655e-03 |  mol |
| :: Na+          | 1.0000e+00 |  mol |
| :: OH-          | 1.6776e-10 |  mol |
| :: CO2(g)       | 9.9293e-07 |  mol |
| :: H2O(g)       | 7.0668e-09 |  mol |
+-----------------+------------+------+

The CO2 saturation pressure for our aqueous solution is printed next:

P = state.pressure() * 1e-5  # convert pressure from Pa to bar
print("Computed saturation pressure is", P, "bar")

Computed saturation pressure is 149.248 bar

And with that, we conclude this tutorial, which demonstrates how Reaktoro can perform chemical equilibrium calculations with given temperature and phase amounts to calculate the CO2 saturation/bubble pressure of an aqueous electrolyte solution.