Activity of Solids

[ghost]

Description

Hello everyone,

I am looking at a simple carbonation reaction at equilibrium (see attached file), where CO2(g) reacts with CaO(s) to form CaCO3(s). If my understanding is correct, the solids have an activity of unity. Hence, the resulting CO2(g) amount should be independent of added CaO(s), as long as enough is present to satisfy the stoichiometry. However, when I run this example in Reaktoro, I find a dependency of the amount of added CaO(s) and the activity of CaCO3(s) is not unity.
Is that a fault in Reaktoro or is there a fault in my understanding?

Example

from reaktoro import *

db = SupcrtDatabase("supcrtbl")
gases = GaseousPhase("CO2(g)")
solid_phase = SolidPhase("Lime Calcite")
system = ChemicalSystem(db, gases, solid_phase)

state = ChemicalState(system)
comp_vec = ["CO2(g)", "Lime"]
n0_vec = [100, 100] # change to n0_vec = [100, 100.5] and the result changes
n0 = dict(zip(comp_vec, n0_vec))
state.temperature(650, "celsius")
state.pressure(0.1013, "MPa")

for comp in comp_vec:
    state.set(comp , n0[comp], "mol")
    
solver = EquilibriumSolver(system)
result = solver.solve(state)  # equilibrate the `state` object!

print(state)
print("\nResult: N_CO2 = {}".format(round(float(state.speciesAmount("CO2(g)")), 2)))

What is your operating system?

Windows

Additional Information

No response

[allanleal]

Hi @Git-User-JW , you are creating a solid solution phase (with ideal thermodynamic model; activity of end member equal its mole fraction).

Instead of

solidphase = SolidPhase("Lime Calcite")

use

solidphases = MineralPhases("Lime Calcite")  # Need to implement SolidPhases! For now, use MineralPhases. No difference. 

The following complete script for convenience:

from reaktoro import *

db = SupcrtDatabase("supcrtbl")
gases = GaseousPhase("CO2(g)")
# solid_phase = SolidPhase("Lime Calcite")
solid_phases = MineralPhases("Lime Calcite")
system = ChemicalSystem(db, gases, solid_phases)

state = ChemicalState(system)
comp_vec = ["CO2(g)", "Lime"]
n0_vec = [100, 100] # change to n0_vec = [100, 100.5] and the result changes
n0 = dict(zip(comp_vec, n0_vec))
state.temperature(650, "celsius")
state.pressure(0.1013, "MPa")

for comp in comp_vec:
    state.set(comp , n0[comp], "mol")

solver = EquilibriumSolver(system)
result = solver.solve(state) # equilibrate the state object!

assert result.succeeded()

print(state)
print("\nResult: N_CO2 = {}".format(round(float(state.speciesAmount("CO2(g)")), 2)))

props = ChemicalProps(state)

print(f'a["Lime"]    = {props.speciesActivity("Lime")}')
print(f'a["Calcite"] = {props.speciesActivity("Calcite")}')
print(f'a["CO2(g)"]  = {props.speciesActivity("CO2(g)")}')

This prints:

+---------------------+------------+------+
| Property            |      Value | Unit |
+---------------------+------------+------+
| Temperature         |   923.1500 |    K |
| Pressure            |     1.0130 |  bar |
| Charge:             | 0.0000e+00 |  mol |
| Element Amount:     |            |      |
| :: C                | 1.0000e+02 |  mol |
| :: O                | 3.0000e+02 |  mol |
| :: Ca               | 1.0000e+02 |  mol |
| Species Amount:     |            |      |
| :: CO2(g)           | 1.0000e-16 |  mol |
| :: Lime :: CaO      | 1.0000e-16 |  mol |
| :: Calcite :: CaCO3 | 1.0000e+02 |  mol |
+---------------------+------------+------+

Result: N_CO2 = 0.0
a["Lime"]    = 1.0
a["Calcite"] = 1.0
a["CO2(g)"]  = 1.0130000000000001

[ghost]

Thank you, I really appreciate the swift reply!
That resolved my confusion with the activity. However, for some reason the result doesn't agree with my hand calculations. On that note, the value for the equilibrium constant was calculated with:

db = SupcrtDatabase("supcrtbl")

T = 650+273.15
p = 1.013
rxn = db.reaction("CO2(g) + Lime = Calcite")
rprops = rxn.props(T, "K", p, "bar")
K = np.exp(-rprops.dG0 / (8.314 * T))
print(K)

giving:

101.09

After some testing, I noticed that the result for N_CO2 now equals zero for t < 890°C, and 100 for t > 890°C. Instead, I would expect a gradual change. Is there something else I am missing?

[allanleal]

Could you please provide your expected K value at that T and P? The value you got is based on supcrtbl database. Different databases can produce different values. For your problem, maybe the NASA database is more appropriate (need to double check; on the phone now). As to the temperature "threshold" above, I'll check later. Do you also have expected behavior for CO2(g)?

[allanleal]

See plot below of the log(K) of the reaction. As you can see, around 890 °C, the sign of log(K) changes. That should explain your observation above. Remember that you are calculating chemical equilibrium and not kinetics. If you want to see gradual changes in the amounts of the species, you'll need to solve a different problem (e.g., a chemical kinetics description or chemical equilibrium in which enthalpy of the system is gradually increased).

This is the code to generate this plot (in a Jupyter notebook):

from reaktoro import *
from reaktplot import *
import numpy as np

db = SupcrtDatabase("supcrtbl")

rnx = db.reaction("CO2(g) + Lime = Calcite")

Ts = np.linspace(100.0, 1000.0, 100)
lgKs = [rnx.props(T, "°C", 1.0, "bar").lgK for T in Ts]

fig = Figure()
fig.drawLineWithMarkers(Ts, lgKs, name="lnK")
fig.title("logK of CO2(g) + CaO(s) = CaCO3(s)")
fig.xaxisTitle("Temperature [°C]")
fig.yaxisTitle("log(K)")
fig.show()

[ghost]

Ah that helped a lot and resolved my problem. I noticed that I had to add an inert component in my example to get the gradual change of nCO2 over the temperature that I expected. Above 890°C, the reverse reaction would occur, so nCO2 of 100mol is not surprising. For future reference, I will attach the corrected code and hand calculation. Thank you @allanleal.

Updated hand calculation:

nCO2 at equilibrium for different temperatures:

generated with:

from reaktoro import *
from reaktplot import *
import numpy as np

db = SupcrtDatabase("supcrtbl")
gases = GaseousPhase("CO2(g) N2(g)")
solid_phases = MineralPhases("Lime Calcite")
system = ChemicalSystem(db, gases, solid_phases)

state = ChemicalState(system)
comp_vec = ["CO2(g)", "N2(g)", "Lime"]
n0_vec = [100, 100, 100] # Same result for n0_Lime > 100 --> check
n0 = dict(zip(comp_vec, n0_vec))

T_vec = np.linspace(100.0, 1000.0, 100)
nCO2_vec = []

for T in T_vec:
    
    state = ChemicalState(system)
    state.temperature(T, "celsius")
    state.pressure(0.1013, "MPa")
    for comp in comp_vec:
        state.set(comp , n0[comp], "mol")
        
    result = solver.solve(state)  # equilibrate the `state` object!
    assert result.succeeded()
    
    nCO2_vec.append(round(float(state.speciesAmount("CO2(g)")), 2))
    
fig = Figure()
fig.drawLineWithMarkers(T_vec, nCO2_vec, name="CO2")
fig.title("nCO2_eq of CO2(g) + CaO(s) = CaCO3(s)")
fig.xaxisTitle("Temperature [°C]")
fig.yaxisTitle("$n_{CO_2, eq}$")
fig.show()
```

[allanleal]

Super interesting. Thanks for posting this!

[ghost]

Just one follow-up question to that, since you suggested to try the NASA database (which I would prefer). While the database file has the information for CaO and CaCO3 in their solid states, I cant define them as mineral phases. Is there a way to do that?

The following code gives the error that no such species could be found:

db = NasaDatabase("nasa-cea")

gases = GaseousPhase(speciate("C H O"))
# condensed_phase = CondensedPhase("CaO(cd) CaCO3(cd)")
solid_phases = MineralPhases("CaO(cd) CaCO3(cd)")

system = ChemicalSystem(db, gases, solid_phases)

[allanleal]

Try CondensedPhases. I need to organize this better!