Defining materials
Contents
Defining materials#
Written by Allan Leal (ETH Zurich) on Jan 31st, 2022.
Last revised on Mar 26th, 2023 by Allan Leal.
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!
When performing chemical equilibrium or chemical kinetics calculations, an initial chemical state is required. From this state, often in disequilibrium, you are interested in:
immediately find a chemical equilibrium state that satisfies the desired constraints; or
compute a sequence of states describing the evolution of the system as it reacts over a period of time.
In most cases, you should be able to construct an initial state in terms of given amounts of chemical species in the system using the ChemicalState class. However, this may not always be convenient or possible for you. In this tutorial, we show an alternative way to define an initial chemical state using the Material class. First, however, we give you an example that shows some disadvantages of working only with the ChemicalState class as a motivation to use Material.
An initial chemical state not so conveniently defined#
Consider an aqueous solution with the following molal composition: 1 molal NaCl, 0.1 molal CaCl2, 0.05 molal MgCl2 and 0.1 molal CO2. Consider a rock with the following mass composition: 80%kg quartz (SiO2) and 20%kg calcite (CaCO3). Mix 1 kg of this aqueous solution with 10 g of rock and find the chemical equilibrium state of the system (i.e., solution and minerals) at 60°C and 15 bar.
As usual, we will start with the definition of a suitable chemical system to model this problem, which is done next.
from reaktoro import *
db = PhreeqcDatabase("phreeqc.dat")
solution = AqueousPhase()
solution.set(ActivityModelPhreeqc(db))
minerals = MineralPhases("Halite Calcite Dolomite Quartz")
system = ChemicalSystem(db, solution, minerals)
Note
In this system, the aqueous species were selected automatically based on the chemical elements composing the explicitly listed minerals (i.e., all aqueous species with elements H, O, Na, Cl, Ca, Mg, and Si). If you are curious, copy the code above and print the species in the chemical system!
for species in system.species():
print(species.name())
Let’s now create our initial chemical containing 1 kg of aqueous solution and 10 g of rock:
state = ChemicalState(system)
# Set initial brine composition
state.set("H2O" , 1.00, "kg")
state.set("Na+" , 1.00, "mol")
state.set("Ca+2", 0.10, "mol")
state.set("Mg+2", 0.05, "mol")
state.set("Cl-" , 1.30, "mol")
state.set("CO2" , 0.10, "mol")
# Set initial rock composition
state.set("Quartz", 80.0, "g")
state.set("Calcite", 20.0, "g")
# Scale fluid and solid masses to 1 kg and 10 g respectively
state.scaleFluidMass(1.0, "kg")
state.scaleSolidMass(10.0, "g")
While defining this initial chemical state, note that:
we needed to specify amounts for ions
Na+,Ca+2, andMg+2that reflect the desired molality of NaCl, CaCl2, and MgCl2.we needed to specify an amount for
Cl-that produces zero electric charge in the solution.we needed to provide species names exactly how they exist in the database instead of chemical formulas (e.g., names
QuartzandCalciteinstead of formulasSiO2andCaCO3).we needed to scale the mass of fluid and solids in the chemical state.
Let’s now equilibrate this state at 60 °C and 15 bar:
state.temperature(60.0, "celsius")
state.pressure(15.0, "bar")
equilibrate(state)
print(state)
+---------------------------+-------------+------+
| Property | Value | Unit |
+---------------------------+-------------+------+
| Temperature | 333.1500 | K |
| Pressure | 15.0000 | bar |
| Charge: | -3.0335e-16 | mol |
| Element Amount: | | |
| :: H | 1.0291e+02 | mol |
| :: C | 1.1269e-01 | mol |
| :: O | 5.1968e+01 | mol |
| :: Na | 9.2704e-01 | mol |
| :: Mg | 4.6352e-02 | mol |
| :: Si | 1.3315e-01 | mol |
| :: Cl | 1.2052e+00 | mol |
| :: Ca | 1.1269e-01 | mol |
| Species Amount: | | |
| :: CO3-2 | 3.4766e-07 | mol |
| :: H+ | 1.1599e-05 | mol |
| :: H2O | 5.1447e+01 | mol |
| :: CO2 | 8.2495e-02 | mol |
| :: (CO2)2 | 4.8663e-04 | mol |
| :: HCO3- | 1.0269e-02 | mol |
| :: CH4 | 1.0000e-16 | mol |
| :: Ca+2 | 9.7795e-02 | mol |
| :: CaCO3 | 4.0323e-06 | mol |
| :: CaHCO3+ | 4.1385e-03 | mol |
| :: CaOH+ | 4.6041e-10 | mol |
| :: Cl- | 1.2052e+00 | mol |
| :: H2 | 1.0000e-16 | mol |
| :: H4SiO4 | 1.8399e-04 | mol |
| :: H2SiO4-2 | 3.7469e-15 | mol |
| :: H3SiO4- | 2.0334e-08 | mol |
| :: Mg+2 | 4.3906e-02 | mol |
| :: MgCO3 | 8.4359e-07 | mol |
| :: MgHCO3+ | 2.4445e-03 | mol |
| :: MgOH+ | 1.2928e-07 | mol |
| :: Na+ | 9.2543e-01 | mol |
| :: NaCO3- | 3.9490e-06 | mol |
| :: NaHCO3 | 1.6087e-03 | mol |
| :: OH- | 1.7999e-08 | mol |
| :: NaOH | 1.0000e-16 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: Halite :: NaCl | 1.0000e-16 | mol |
| :: Calcite :: CaCO3 | 1.0748e-02 | mol |
| :: Dolomite :: CaMg(CO3)2 | 1.2258e-15 | mol |
| :: Quartz :: SiO2 | 1.3296e-01 | mol |
+---------------------------+-------------+------+
We managed to get the problem defined and solved, but the coding process could be simpler as shown in the next section.
Simplifying with class Material#
Let’s now use class Material to obtain the same result with more convenience:
brine = Material(system)
brine.add("H2O" , 1.00, "kg")
brine.add("NaCl" , 1.00, "mol")
brine.add("CaCl2", 0.10, "mol")
brine.add("MgCl2", 0.05, "mol")
brine.add("CO2" , 0.10, "mol")
rock = Material(system)
rock.add("SiO2", 80.0, "g")
rock.add("CaCO3", 20.0, "g")
mix = brine(1.0, "kg") + rock(10.0, "g")
state = mix.equilibrate(60, "celsius", 15, "bar")
print(state)
+---------------------------+-------------+------+
| Property | Value | Unit |
+---------------------------+-------------+------+
| Temperature | 333.1500 | K |
| Pressure | 15.0000 | bar |
| Charge: | -5.2987e-17 | mol |
| Element Amount: | | |
| :: H | 1.0291e+02 | mol |
| :: C | 1.1269e-01 | mol |
| :: O | 5.1968e+01 | mol |
| :: Na | 9.2704e-01 | mol |
| :: Mg | 4.6352e-02 | mol |
| :: Si | 1.3315e-01 | mol |
| :: Cl | 1.2052e+00 | mol |
| :: Ca | 1.1269e-01 | mol |
| Species Amount: | | |
| :: CO3-2 | 3.4766e-07 | mol |
| :: H+ | 1.1599e-05 | mol |
| :: H2O | 5.1447e+01 | mol |
| :: CO2 | 8.2495e-02 | mol |
| :: (CO2)2 | 4.8663e-04 | mol |
| :: HCO3- | 1.0269e-02 | mol |
| :: CH4 | 1.0000e-16 | mol |
| :: Ca+2 | 9.7795e-02 | mol |
| :: CaCO3 | 4.0323e-06 | mol |
| :: CaHCO3+ | 4.1385e-03 | mol |
| :: CaOH+ | 4.6041e-10 | mol |
| :: Cl- | 1.2052e+00 | mol |
| :: H2 | 1.0000e-16 | mol |
| :: H4SiO4 | 1.8399e-04 | mol |
| :: H2SiO4-2 | 3.7469e-15 | mol |
| :: H3SiO4- | 2.0334e-08 | mol |
| :: Mg+2 | 4.3907e-02 | mol |
| :: MgCO3 | 8.4359e-07 | mol |
| :: MgHCO3+ | 2.4445e-03 | mol |
| :: MgOH+ | 1.2928e-07 | mol |
| :: Na+ | 9.2543e-01 | mol |
| :: NaCO3- | 3.9490e-06 | mol |
| :: NaHCO3 | 1.6087e-03 | mol |
| :: OH- | 1.7999e-08 | mol |
| :: NaOH | 1.0000e-16 | mol |
| :: O2 | 1.0000e-16 | mol |
| :: Halite :: NaCl | 1.0000e-16 | mol |
| :: Calcite :: CaCO3 | 1.0748e-02 | mol |
| :: Dolomite :: CaMg(CO3)2 | 1.2258e-15 | mol |
| :: Quartz :: SiO2 | 1.3296e-01 | mol |
+---------------------------+-------------+------+
The code above should (hopefully!) be intuitive enough and self-explanatory. The table printed above should also be identical (or extremely close!) to the previous one.
That’s it; you now have an extra Reaktoro component in your arsenal when performing your chemical modeling. The Material class should provide that extra level of convenience if things get a little more confusing if only the ChemicalState class is used.
Attention
You may have to use the Material class more carefully if chemical kinetics is considered. In this case, it is critical that species undergoing chemical reactions controlled by kinetic rates must have their initial amounts explicitly specified in the ChemicalState object, using their names instead of their formulas to avoid ambiguity. For example, suppose the minerals calcite and quartz are to react kinetically with the aqueous solution (and the aqueous species react instantaneously at equilibrium). In that case, the initial amounts of these minerals must be given in the ChemicalState object. The code below should give you an idea of how to proceed in this scenario:
brine = Material(system)
brine.add("H2O" , 1.00, "kg")
brine.add("NaCl" , 1.00, "mol")
brine.add("CaCl2", 0.10, "mol")
brine.add("MgCl2", 0.05, "mol")
brine.add("CO2" , 0.10, "mol")
state = brine(1.0, "kg").equilibrate(60, "celsius", 15, "bar")
state.set("Quartz", 80.0, "g")
state.set("Calcite", 20.0, "g")
state.scaleSolidMass(10.0, "g")
# use the above chemical state as the initial condition
# for your chemical kinetics calculations below