Specification of binary interaction parameters.

[legofernando1999]

Description

Hello, I am doing some equilibrium calculations using the PR EOS and I need to specify the BIPs. How can I do that?

Example

No response

What is your operating system?

Windows

Additional Information

No response

[allanleal]

Hi @legofernando1999 - below is an example that defines the Søreide & Whitson (1992) BIP model in Python and sets it to Reaktoro's Peng-Robinson model. It's based on the C++ implementation in Reaktoro:

from reaktoro import *
import numpy as np


def bip_model_soreide_whitson(specieslist: SpeciesList):
    substances = [x.formula().str() for x in specieslist]

    def index_substance(formula: str):
        try: return substances.index(formula) 
        except: return -1  

    iH2O    = index_substance("H2O")
    iCO2    = index_substance("CO2")
    iH2S    = index_substance("H2S")
    iCH4    = index_substance("CH4")
    iN2     = index_substance("N2")
    iC2H6   = index_substance("C2H6")
    iC3H8   = index_substance("C3H8")
    inC4H10 = index_substance("nC4H10")

    kH2O_CO2 = 0.1896  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/CO₂.
    kH2O_N2 = 0.4778  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/N₂.
    kH2O_CH4 = 0.4850  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/CH₄.
    kH2O_C2H6 = 0.4920  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/C₂H₆.
    kH2O_C3H8 = 0.5525  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/C₃H₈.
    kH2O_nC4H10 = 0.5525  # The binary interaction parameter kᵢⱼ for the substance pair H₂O/n-C₄H₁₀.
    kH2O_H2S_a1 = 0.19031  # The coefficient a₁ when computing the binary interaction parameter kᵢⱼ for the substance pair H₂O/H₂S.
    kH2O_H2S_a2 = -0.05965  # The coefficient a₂ when computing the binary interaction parameter kᵢⱼ for the substance pair H₂O/H₂S.

    size = len(substances)
    bip = CubicEOS.Bip(size)

    k, kT, kTT = np.zeros((size, size)), np.zeros((size, size)), np.zeros((size, size))

    # Define the function that will be called whenever the BIP parameters need to be updated
    def evalfn(args):
        T = args.T  # the temperature (in K)

        Tcr = args.Tcr  # the critical temperatures of the substances in the phase (in K)
        Pcr = args.Pcr  # the critical pressures of the substances in the phase (in Pa)
        omega = args.omega  # the acentric factors of the substances in the phase

        if iH2O < size:
            if  iCO2   != -1: k[iH2O,    iCO2] = k[   iCO2, iH2O] = kH2O_CO2
            if   iN2   != -1: k[iH2O,     iN2] = k[    iN2, iH2O] = kH2O_N2
            if  iH2S   != -1: k[iH2O,    iH2S] = k[   iH2S, iH2O] = kH2O_H2S_a1 + kH2O_H2S_a2 * float(T/Tcr[iH2S])  # For now, convert autodiff.real to float to avoid error with numpy
            if  iCH4   != -1: k[iH2O,    iCH4] = k[   iCH4, iH2O] = kH2O_CH4
            if iC2H6   != -1: k[iH2O,   iC2H6] = k[  iC2H6, iH2O] = kH2O_C2H6
            if iC3H8   != -1: k[iH2O,   iC3H8] = k[  iC3H8, iH2O] = kH2O_C3H8
            if inC4H10 != -1: k[iH2O, inC4H10] = k[inC4H10, iH2O] = kH2O_nC4H10

        bip.set(k, kT, kTT)  # This is not available in current version of Reaktoro 2.11.1

        return bip

    return evalfn


db = PhreeqcDatabase("phreeqc.dat")

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

gases = GaseousPhase("CO2(g) H2O(g)")
gases.set(ActivityModelPengRobinson78(bip_model_soreide_whitson))

system = ChemicalSystem(db, solution, gases)

state = ChemicalState(system)
state.temperature(25.0, "celsius")
state.pressure(1.0, "bar")
state.set("H2O"   , 1.0 , "kg")
state.set("CO2(g)", 10.0, "mol")
state.set("Na+"   , 4.00, "mol")
state.set("Cl-"   , 4.00, "mol")

solver = EquilibriumSolver(system)
solver.solve(state)

print(state)

This code produces the following output:

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   298.1500 |    K |
| Pressure        |     1.0000 |  bar |
| Charge:         | 1.2369e-16 |  mol |
| Element Amount: |            |      |
| :: H            | 1.1101e+02 |  mol |
| :: C            | 1.0000e+01 |  mol |
| :: O            | 7.5506e+01 |  mol |
| :: Na           | 4.0000e+00 |  mol |
| :: Cl           | 4.0000e+00 |  mol |
| Species Amount: |            |      |
| :: CO3-2        | 1.9273e-10 |  mol |
| :: H+           | 2.1727e-04 |  mol |
| :: H2O          | 5.5224e+01 |  mol |
| :: CO2          | 1.2975e-02 |  mol |
| :: (CO2)2       | 7.8382e-06 |  mol |
| :: HCO3-        | 1.3556e-04 |  mol |
| :: CH4          | 1.0000e-16 |  mol |
| :: Cl-          | 4.0000e+00 |  mol |
| :: H2           | 1.0000e-16 |  mol |
| :: Na+          | 3.9999e+00 |  mol |
| :: NaCO3-       | 1.0939e-09 |  mol |
| :: NaHCO3       | 8.1708e-05 |  mol |
| :: OH-          | 1.1399e-10 |  mol |
| :: NaOH         | 1.0000e-16 |  mol |
| :: O2           | 1.0000e-16 |  mol |
| :: CO2(g)       | 9.9868e+00 |  mol |
| :: H2O(g)       | 2.8219e-01 |  mol |
+-----------------+------------+------+

and it is equivalent if we had defined the gas phase with:

gases.set(ActivityModelPengRobinsonSoreideWhitson())

Note the comment above:

bip.set(k, kT, kTT)  # This is not available in current version of Reaktoro 2.11.1

I'm preparing a release for v2.12.0 which should contain this set method. Let me know if you need clarifications.

• Søreide, I., & Whitson, H. C. (1992). Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2S with pure water and NaCI brine. Fluid Phase Equilibria, 77(2), 217–240. https://doi.org/10.1016/0378-3812(92)85105-H

[legofernando1999]

Hi @allanleal, thanks for the example, it was very helpful.
One more question, since bip.set(k, kT, kTT) is not available in the current version, is there another way for setting the BIPs?

[allanleal]

Right now, only using C++. Without any changes in the current codebase, pybind11 tries to convert k, a C++ Eigen::Matrix with autodiff::real numbers, into a numpy array. Numpy complains that it does not know autodiff::real.

[legofernando1999]

I see. Thank you for your time.