co2ch4_solubility/solCode_mixture.py at main · sinmec/co2ch4_solubility · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
from reaktoro import *
import CoolProp.CoolProp as CP
import pandas as pd

# Database
db = PhreeqcDatabase("phreeqc.dat")
aqueousphase = AqueousPhase(speciate("H O C Na Cl"))
aqueousphase.set(ActivityModelPitzer())
gaseousphase = GaseousPhase("CO2(g) CH4(g)")
gaseousphase.set(ActivityModelPengRobinson())
phases = Phases(db)
phases.add(aqueousphase)
phases.add(gaseousphase)
system = ChemicalSystem(phases)
solver = EquilibriumSolver(system)

def compute_mole_fractions(T, P, n_NaCl, y_co2g, y_ch4g):
    # Molar mass [kg/mol]
    mm_CO2 = CP.PropsSI('MOLARMASS', 'T', T+273.15, 'P', P*1e5, 'CO2')
    mm_CH4 = CP.PropsSI('MOLARMASS', 'T', T+273.15, 'P', P*1e5, 'CH4')
    mm_H2O = CP.PropsSI('MOLARMASS', 'T', T+273.15, 'P', P*1e5, 'Water')

    # Initial quantities (mol)
    mH2O = 0.25/100
    mNa = n_NaCl * (mH2O / 55.5)
    mCl = n_NaCl * (mH2O / 55.5)
    mCO2 = 0 #2.2595 / (100 * 0.044) * (mH2O / 55.5)
    n0CO2g = (1 - mH2O - mNa - mCl - mCO2) * y_co2g
    n0CH4g = (1 - mH2O - mNa - mCl - mCO2) * y_ch4g

    # Chemical state
    state = ChemicalState(system)
    state.temperature(T, "celsius")
    state.pressure(P, "bar")
    state.set("H2O", mH2O, "mol")
    state.set("CO2(g)", n0CO2g, "mol")
    state.set("CH4(g)", n0CH4g, "mol")
    state.set("Na+", mNa, "mol")
    state.set("Cl-", mCl, "mol")
    #state.set("CO2", mCO2, "mol")

    res = solver.solve(state)
    assert res.succeeded(), "Equilibrium failed!"

    # Properties of aqueous phase
    aqprops = AqueousProps(state)
    molality_CO2 = aqprops.speciesMolality("CO2") + aqprops.speciesMolality("(CO2)2") + aqprops.speciesMolality("CO3-2") + aqprops.speciesMolality("HCO3-")
    molality_CH4 = aqprops.speciesMolality("CH4")
    molality_NaCl = aqprops.speciesMolality("Na+") + aqprops.speciesMolality("Cl-")
    molality_H2O = aqprops.speciesMolality("H2O")
    molality_total = molality_CO2 + molality_CH4 + molality_NaCl + molality_H2O

    # Molar fractions
    x_CO2 = float(molality_CO2 / molality_total)
    x_CH4 = float(molality_CH4 / molality_total)
    x_NaCl = float(molality_NaCl / molality_total)
    x_H2O = 1.0 - x_CO2 - x_CH4 - x_NaCl

    # Solubulity: mass and molal
    mass_CO2 = mm_CO2 * x_CO2
    mass_H2O = mm_H2O * x_H2O
    sol_CO2 = (mass_CO2 / mass_H2O) * 100
    molal_CO2 = x_CO2 / mass_H2O

    return {
        "Temperature [°C]": T,
        "Pressure [bar]": P,
        "x_CO2 [%]": x_CO2 * 100,
        "x_CH4 [%]": x_CH4 * 100,
        "molal [mol CO2 / kg H2O]": molal_CO2,
        "Solubility CO2 [kg/100kg H2O]": sol_CO2
    }


def run_from_table(tabela: pd.DataFrame):
    """Executa a rotina para todos os casos da tabela."""
    results = []
    for _, row in tabela.iterrows():
        res = compute_mole_fractions(row["T"], row["P"], row["n_NaCl"], row["y_co2g"], row["y_ch4g"])
        results.append(res)
    return results