Modeling of a VLL system. #377

legofernando1999 · 2024-04-05T21:13:22Z

legofernando1999
Apr 5, 2024

Description

Hello, I am doing flash calculations for the following oil,

For instance, using the PR EOS at a temperature and pressure of 105°F and 1295 psi, I should get three phases with this composition

but instead I'm getting a single liquid phase.
Ignoring for the moment the BIPs and focusing on qualitative agreement, why am I not getting three phases?
Please let me know if something needs clarification.

Example

Here is what I did:

import reaktoro as rkt

db = rkt.SupcrtDatabase("supcrtbl-organics")

# Define species

# ==== C2-C3 ====
db.addElement(rkt.Element()
              .withName("C2-C3")
              .withSymbol("C2-C3")
              .withMolarMass(0.0372))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C2-C3(l)")
              .withFormula("C2-C3")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C2-C3"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C2-C3(l)")
crprops.setTemperature(344.148, "K")
crprops.setPressure(44.992, "bar")
crprops.setAcentricFactor(0.131)
rkt.CriticalProps.append(crprops)

# ==== C4-C6 ====
db.addElement(rkt.Element()
              .withName("C4-C6")
              .withSymbol("C4-C6")
              .withMolarMass(0.0695))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C4-C6(l)")
              .withFormula("C4-C6")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C4-C6"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C4-C6(l)")
crprops.setTemperature(463.222, "K")
crprops.setPressure(33.996, "bar")
crprops.setAcentricFactor(0.24)
rkt.CriticalProps.append(crprops)

# ==== C7-C15 ====
db.addElement(rkt.Element()
              .withName("C7-C15")
              .withSymbol("C7-C15")
              .withMolarMass(0.14096))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C7-C15(l)")
              .withFormula("C7-C15")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C7-C15"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C7-C15(l)")
crprops.setTemperature(605.694, "K")
crprops.setPressure(21.749, "bar")
crprops.setAcentricFactor(0.618)
rkt.CriticalProps.append(crprops)

# ==== C16-C27 ====
db.addElement(rkt.Element()
              .withName("C16-C27")
              .withSymbol("C16-C27")
              .withMolarMass(0.28099))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C16-C27(l)")
              .withFormula("C16-C27")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C16-C27"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C16-C27(l)")
crprops.setTemperature(751.017, "K")
crprops.setPressure(16.541, "bar")
crprops.setAcentricFactor(0.957)
rkt.CriticalProps.append(crprops)

# ==== C28+ ====
db.addElement(rkt.Element()
              .withName("C28+")
              .withSymbol("C28+")
              .withMolarMass(0.51962))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C28+(l)")
              .withFormula("C28+")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C28+"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C28+(l)")
crprops.setTemperature(942.478, "K")
crprops.setPressure(16.418, "bar")
crprops.setAcentricFactor(1.268)
rkt.CriticalProps.append(crprops)

# Duplicate species to have gaseous and liquid phases

co2_l = db.species("CO2(g)").clone()
co2_l = co2_l.withName("CO2(l)")
co2_l = co2_l.withAggregateState(rkt.AggregateState.Liquid)
db.addSpecies(co2_l)

ch4_l = db.species("CH4(g)").clone()
ch4_l = ch4_l.withName("CH4(l)")
ch4_l = ch4_l.withAggregateState(rkt.AggregateState.Liquid)
db.addSpecies(ch4_l)

c2c3_g = db.species("C2-C3(l)").clone()
c2c3_g = c2c3_g.withName("C2-C3(g)")
c2c3_g = c2c3_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(c2c3_g)

c4c6_g = db.species("C4-C6(l)").clone()
c4c6_g = c4c6_g.withName("C4-C6(g)")
c4c6_g = c4c6_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(c4c6_g)

c7c15_g = db.species("C7-C15(l)").clone()
c7c15_g = c7c15_g.withName("C7-C15(g)")
c7c15_g = c7c15_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(c7c15_g)

c16c27_g = db.species("C16-C27(l)").clone()
c16c27_g = c16c27_g.withName("C16-C27(g)")
c16c27_g = c16c27_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(c16c27_g)

c28plus_g = db.species("C28+(l)").clone()
c28plus_g = c28plus_g.withName("C28+(g)")
c28plus_g = c28plus_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(c28plus_g)

# Set up phases in chemical system
gas = rkt.GaseousPhase("CO2(g) CH4(g) C2-C3(g) C4-C6(g) C7-C15(g) C16-C27(g) C28+(g)")
gas.set(rkt.ActivityModelPengRobinson())

liquid_1 = rkt.LiquidPhase("CO2(l) CH4(l) C2-C3(l) C4-C6(l) C7-C15(l) C16-C27(l) C28+(l)")
liquid_1.set(rkt.ActivityModelPengRobinson())

liquid_2 = rkt.LiquidPhase("CO2(l) CH4(l) C2-C3(l) C4-C6(l) C7-C15(l) C16-C27(l) C28+(l)")
liquid_2.set(rkt.ActivityModelPengRobinson())

# Define chemical system
phases = rkt.Phases(db)
phases.add(gas)
phases.add(liquid_1)
phases.add(liquid_2)
system = rkt.ChemicalSystem(phases)

oil_comps = ["CO2(g)", "CH4(g)", "C2-C3(g)", "C4-C6(g)", "C7-C15(g)", "C16-C27(g)", "C28+(g)",
             "CO2(l)", "CH4(l)", "C2-C3(l)", "C4-C6(l)", "C7-C15(l)", "C16-C27(l)", "C28+(l)",
             "CO2(l)!", "CH4(l)!", "C2-C3(l)!", "C4-C6(l)!", "C7-C15(l)!", "C16-C27(l)!", "C28+(l)!"]
oil_moles = [0.74337, 0.15894, 0.07527, 0.02086, 0.00155, 6.2515e-6, 9.0351e-10,
             0.53522, 0.05410, 0.10873, 0.10614, 0.12903, 0.04750, 0.01928,
             0.55099, 0.05592, 0.10957, 0.10546, 0.12236, 0.04173, 0.01399]
oil_init_amounts = dict(zip(oil_comps, oil_moles))

# Create a chemical state object
state = rkt.ChemicalState(system)
state.setTemperature(105., "F")
state.pressure(1295., "psi")
for comp in oil_comps:
    state.set(comp, oil_init_amounts[comp], "mol")

print(state)  # Non-equilibrated state
print()

rkt.equilibrate(state)

print(state)  # State after equilibration

This code produces the following output:

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   313.7056 |    K |
| Pressure        |    89.2871 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: H            | 1.0758e+00 |  mol |
| :: C            | 2.0985e+00 |  mol |
| :: O            | 3.6592e+00 |  mol |
| :: C2-C3        | 2.9357e-01 |  mol |
| :: C4-C6        | 2.3246e-01 |  mol |
| :: C7-C15       | 2.5294e-01 |  mol |
| :: C16-C27      | 8.9236e-02 |  mol |
| :: C28+         | 3.3270e-02 |  mol |
| Species Amount: |            |      |
| :: CO2(g)       | 5.3522e-01 |  mol |
| :: CH4(g)       | 5.4100e-02 |  mol |
| :: C2-C3(g)     | 1.0873e-01 |  mol |
| :: C4-C6(g)     | 1.0614e-01 |  mol |
| :: C7-C15(g)    | 1.2903e-01 |  mol |
| :: C16-C27(g)   | 4.7500e-02 |  mol |
| :: C28+(g)      | 1.9280e-02 |  mol |
| :: CO2(l)       | 7.4337e-01 |  mol |
| :: CH4(l)       | 1.5894e-01 |  mol |
| :: C2-C3(l)     | 7.5270e-02 |  mol |
| :: C4-C6(l)     | 2.0860e-02 |  mol |
| :: C7-C15(l)    | 1.5500e-03 |  mol |
| :: C16-C27(l)   | 6.2515e-06 |  mol |
| :: C28+(l)      | 9.0351e-10 |  mol |
| :: CO2(l)!      | 5.5099e-01 |  mol |
| :: CH4(l)!      | 5.5920e-02 |  mol |
| :: C2-C3(l)!    | 1.0957e-01 |  mol |
| :: C4-C6(l)!    | 1.0546e-01 |  mol |
| :: C7-C15(l)!   | 1.2236e-01 |  mol |
| :: C16-C27(l)!  | 4.1730e-02 |  mol |
| :: C28+(l)!     | 1.3990e-02 |  mol |
+-----------------+------------+------+

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   313.7056 |    K |
| Pressure        |    89.2871 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: H            | 1.0758e+00 |  mol |
| :: C            | 2.0985e+00 |  mol |
| :: O            | 3.6592e+00 |  mol |
| :: C2-C3        | 2.9357e-01 |  mol |
| :: C4-C6        | 2.3246e-01 |  mol |
| :: C7-C15       | 2.5294e-01 |  mol |
| :: C16-C27      | 8.9236e-02 |  mol |
| :: C28+         | 3.3270e-02 |  mol |
| Species Amount: |            |      |
| :: CO2(g)       | 5.4992e-15 |  mol |
| :: CH4(g)       | 8.0842e-16 |  mol |
| :: C2-C3(g)     | 8.8239e-16 |  mol |
| :: C4-C6(g)     | 6.9871e-16 |  mol |
| :: C7-C15(g)    | 7.6026e-16 |  mol |
| :: C16-C27(g)   | 2.6822e-16 |  mol |
| :: C28+(g)      | 1.0000e-16 |  mol |
| :: CO2(l)       | 5.4992e-15 |  mol |
| :: CH4(l)       | 8.0842e-16 |  mol |
| :: C2-C3(l)     | 8.8239e-16 |  mol |
| :: C4-C6(l)     | 6.9871e-16 |  mol |
| :: C7-C15(l)    | 7.6026e-16 |  mol |
| :: C16-C27(l)   | 2.6822e-16 |  mol |
| :: C28+(l)      | 1.0000e-16 |  mol |
| :: CO2(l)!      | 1.8296e+00 |  mol |
| :: CH4(l)!      | 2.6896e-01 |  mol |
| :: C2-C3(l)!    | 2.9357e-01 |  mol |
| :: C4-C6(l)!    | 2.3246e-01 |  mol |
| :: C7-C15(l)!   | 2.5294e-01 |  mol |
| :: C16-C27(l)!  | 8.9236e-02 |  mol |
| :: C28+(l)!     | 3.3270e-02 |  mol |
+-----------------+------------+------+

What is your operating system?

Windows

Additional Information

No response

allanleal · 2024-04-09T16:58:22Z

allanleal
Apr 9, 2024
Maintainer

I see above that you are zeroing out the standard state Gibbs energies of the new species (the ones you are defining) while using CO2(g) and CH4(g) from the SUPCRT database. This is not thermodynamically consistent. For your use case above, you want to define all species and set their default Gibbs energies to zero.

Your chemical species must also be defined in terms of components (e.g., the species CO2(g) and CO2(l) would be defined in terms of a fictitious “element” CO2 rather than the chemical elements C and O). If you don't do this, chemical reactions between the species will be possible and you won't have the same underlying mathematical problem as a flash calculation (in which the chemical potential of the species are equated).

Reaktoro needs an interface to enable users to set up VLLE and other phase equilibrium problems more easily.

5 replies

legofernando1999 Apr 10, 2024
Author

Hi @allanleal, I did as you told me. I defined CO2 and CH4 and set their default Gibbs energies to zero.

# ==== CO2 ====
db.addElement(rkt.Element()
              .withName("CO2")
              .withSymbol("CO2")
              .withMolarMass(0.04401))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("CO2(l)")
              .withFormula("CO2")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("CO2"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("CO2(l)")
crprops.setTemperature(304.21, "K")
crprops.setPressure(73.83046125, "bar")
crprops.setAcentricFactor(0.2236)
rkt.CriticalProps.overwrite(crprops)

# ==== CH4 ====
db.addElement(rkt.Element()
              .withName("CH4")
              .withSymbol("CH4")
              .withMolarMass(0.016043))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("CH4(l)")
              .withFormula("CH4")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("CH4"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("CH4(l)")
crprops.setTemperature(190.564, "K")
crprops.setPressure(45.99040425, "bar")
crprops.setAcentricFactor(0.0115)
rkt.CriticalProps.overwrite(crprops)

Then I added them to the database as gas

co2_g = db.species("CO2(l)").clone()
co2_g = co2_g.withName("CO2(gas)")
co2_g = co2_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(co2_g)

ch4_g = db.species("CH4(l)").clone()
ch4_g = ch4_g.withName("CH4(gas)")
ch4_g = ch4_g.withAggregateState(rkt.AggregateState.Gas)
db.addSpecies(ch4_g)

After that, I redefined my chemical system

# Set up phases in chemical system
gas = rkt.GaseousPhase("CO2(gas) CH4(gas) C2-C3(g) C4-C6(g) C7-C15(g) C16-C27(g) C28+(g)")
gas.set(rkt.ActivityModelPengRobinson())

liquid_1 = rkt.LiquidPhase("CO2(l) CH4(l) C2-C3(l) C4-C6(l) C7-C15(l) C16-C27(l) C28+(l)")
liquid_1.set(rkt.ActivityModelPengRobinson())

liquid_2 = rkt.LiquidPhase("CO2(l) CH4(l) C2-C3(l) C4-C6(l) C7-C15(l) C16-C27(l) C28+(l)")
liquid_2.set(rkt.ActivityModelPengRobinson())

# Define chemical system
phases = rkt.Phases(db)
phases.add(gas)
phases.add(liquid_1)
phases.add(liquid_2)
system = rkt.ChemicalSystem(phases)

oil_comps = ["CO2(gas)", "CH4(gas)", "C2-C3(g)", "C4-C6(g)", "C7-C15(g)", "C16-C27(g)", "C28+(g)",
             "CO2(l)", "CH4(l)", "C2-C3(l)", "C4-C6(l)", "C7-C15(l)", "C16-C27(l)", "C28+(l)",
             "CO2(l)!", "CH4(l)!", "C2-C3(l)!", "C4-C6(l)!", "C7-C15(l)!", "C16-C27(l)!", "C28+(l)!"]
oil_moles = [0.74337, 0.15894, 0.07527, 0.02086, 0.00155, 6.2515e-6, 9.0351e-10,
             0.53522, 0.05410, 0.10873, 0.10614, 0.12903, 0.04750, 0.01928,
             0.55099, 0.05592, 0.10957, 0.10546, 0.12236, 0.04173, 0.01399]

Finally, I ran the code and got the following output

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   313.7056 |    K |
| Pressure        |    89.2871 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: CH4          | 2.6896e-01 |  mol |
| :: C2-C3        | 2.9357e-01 |  mol |
| :: CO2          | 1.8296e+00 |  mol |
| :: C4-C6        | 2.3246e-01 |  mol |
| :: C7-C15       | 2.5294e-01 |  mol |
| :: C16-C27      | 8.9236e-02 |  mol |
| :: C28+         | 3.3270e-02 |  mol |
| Species Amount: |            |      |
| :: CO2(gas)     | 5.4992e-15 |  mol |
| :: CH4(gas)     | 8.0842e-16 |  mol |
| :: C2-C3(g)     | 8.8239e-16 |  mol |
| :: C4-C6(g)     | 6.9871e-16 |  mol |
| :: C7-C15(g)    | 7.6026e-16 |  mol |
| :: C16-C27(g)   | 2.6822e-16 |  mol |
| :: C28+(g)      | 1.0000e-16 |  mol |
| :: CO2(l)       | 5.4992e-15 |  mol |
| :: CH4(l)       | 8.0842e-16 |  mol |
| :: C2-C3(l)     | 8.8239e-16 |  mol |
| :: C4-C6(l)     | 6.9871e-16 |  mol |
| :: C7-C15(l)    | 7.6026e-16 |  mol |
| :: C16-C27(l)   | 2.6822e-16 |  mol |
| :: C28+(l)      | 1.0000e-16 |  mol |
| :: CO2(l)!      | 1.8296e+00 |  mol |
| :: CH4(l)!      | 2.6896e-01 |  mol |
| :: C2-C3(l)!    | 2.9357e-01 |  mol |
| :: C4-C6(l)!    | 2.3246e-01 |  mol |
| :: C7-C15(l)!   | 2.5294e-01 |  mol |
| :: C16-C27(l)!  | 8.9236e-02 |  mol |
| :: C28+(l)!     | 3.3270e-02 |  mol |
+-----------------+------------+------+

allanleal Apr 11, 2024
Maintainer

Thanks for the updates. At the moment, equilibrium calculations in which all phases (V, L1, L2) are set with a common thermodynamic model (e.g., Peng-Robinson set to the vapor, CO2-rich-liquid, and oil phases) may converge to a wrong solution because the same model/EOS brings no differentiation among the phases. Further work in Reaktoro is needed to ensure that Peng-Robinson assigned to the V-phase behaves differently than that for the L1 and L2 phases. The current implementation of these cubic equations of state in Reaktoro is such that: given (T, P, x), the model outputs V and other properties (e.g., fugacity coefficients) along with a flag indicating what is the stable physical state of the phase: Gas, Liquid, or Supercritical. So, the model is not currently coded to produce liquid-like or gas-like properties, despite the phases being stable/unstable. It always produces properties for the most stable state.

Currently, to get these calculations working, one needs to make some restrictions on the final equilibrium amounts of some species. For example, in the code below, I could get three stable phases (but this required an ad-hoc restriction on the amounts of CO2(l) and CO2(oil)). I also added a bip model. Note: The final result does not compare with the table above. Maybe my feed composition is not right. Would need to check this later.

Output for the Python code below

+-----------------+------------+------+
| Property        |      Value | Unit |
+-----------------+------------+------+
| Temperature     |   313.7056 |    K |
| Pressure        |    89.2871 |  bar |
| Charge:         | 0.0000e+00 |  mol |
| Element Amount: |            |      |
| :: CH4          | 7.3050e-02 |  mol |
| :: C2-C3        | 1.0297e-01 |  mol |
| :: CO2          | 5.7529e-01 |  mol |
| :: C4-C6        | 9.0880e-02 |  mol |
| :: C7-C15       | 1.0509e-01 |  mol |
| :: C16-C27      | 3.7920e-02 |  mol |
| :: C28+         | 1.4800e-02 |  mol |
| Species Amount: |            |      |
| :: CO2(gas)     | 6.5290e-02 |  mol |
| :: CH4(gas)     | 8.1885e-03 |  mol |
| :: C2-C3(gas)   | 1.3615e-02 |  mol |
| :: C4-C6(gas)   | 1.3161e-02 |  mol |
| :: C7-C15(gas)  | 2.0221e-02 |  mol |
| :: C16-C27(gas) | 1.1320e-02 |  mol |
| :: C28+(gas)    | 9.9121e-03 |  mol |
| :: CO2(l)       | 1.0000e-02 |  mol |
| :: CH4(l)       | 1.8735e-03 |  mol |
| :: C2-C3(l)     | 1.0581e-03 |  mol |
| :: C4-C6(l)     | 3.6883e-04 |  mol |
| :: C7-C15(l)    | 5.2113e-05 |  mol |
| :: C16-C27(l)   | 5.8565e-07 |  mol |
| :: C28+(l)      | 2.8946e-10 |  mol |
| :: CO2(oil)     | 5.0000e-01 |  mol |
| :: CH4(oil)     | 6.2988e-02 |  mol |
| :: C2-C3(oil)   | 8.8297e-02 |  mol |
| :: C4-C6(oil)   | 7.7350e-02 |  mol |
| :: C7-C15(oil)  | 8.4817e-02 |  mol |
| :: C16-C27(oil) | 2.6600e-02 |  mol |
| :: C28+(oil)    | 4.8879e-03 |  mol |
+-----------------+------------+------+

Mole Fractions
CO2(gas)            : 0.460739
CH4(gas)            : 0.057785
C2-C3(gas)          : 0.096079
C4-C6(gas)          : 0.092874
C7-C15(gas)         : 0.142695
C16-C27(gas)        : 0.079880
C28+(gas)           : 0.069948
CO2(l)              : 0.748884
CH4(l)              : 0.140307
C2-C3(l)            : 0.079242
C4-C6(l)            : 0.027621
C7-C15(l)           : 0.003903
C16-C27(l)          : 0.000044
C28+(l)             : 0.000000
CO2(oil)            : 0.591758
CH4(oil)            : 0.074547
C2-C3(oil)          : 0.104501
C4-C6(oil)          : 0.091545
C7-C15(oil)         : 0.100382
C16-C27(oil)        : 0.031481
C28+(oil)           : 0.005785

import reaktoro as rkt


# Define species
db = rkt.Database()

# ==== CO2 ====
db.addElement(rkt.Element()
              .withName("CO2")
              .withSymbol("CO2")
              .withMolarMass(0.04401))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("CO2(l)")
              .withFormula("CO2")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("CO2"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("CO2(l)")
crprops.setTemperature(304.21, "K")
crprops.setPressure(73.83046125, "bar")
crprops.setAcentricFactor(0.2236)
rkt.CriticalProps.overwrite(crprops)

# ==== CH4 ====
db.addElement(rkt.Element()
              .withName("CH4")
              .withSymbol("CH4")
              .withMolarMass(0.016043))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("CH4(l)")
              .withFormula("CH4")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("CH4"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("CH4(l)")
crprops.setTemperature(190.564, "K")
crprops.setPressure(45.99040425, "bar")
crprops.setAcentricFactor(0.0115)
rkt.CriticalProps.overwrite(crprops)

# ==== C2-C3 ====
db.addElement(rkt.Element()
              .withName("C2-C3")
              .withSymbol("C2-C3")
              .withMolarMass(0.0372))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C2-C3(l)")
              .withFormula("C2-C3")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C2-C3"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C2-C3(l)")
crprops.setTemperature(344.148, "K")
crprops.setPressure(44.992, "bar")
crprops.setAcentricFactor(0.131)
rkt.CriticalProps.append(crprops)

# ==== C4-C6 ====
db.addElement(rkt.Element()
              .withName("C4-C6")
              .withSymbol("C4-C6")
              .withMolarMass(0.0695))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C4-C6(l)")
              .withFormula("C4-C6")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C4-C6"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C4-C6(l)")
crprops.setTemperature(463.222, "K")
crprops.setPressure(33.996, "bar")
crprops.setAcentricFactor(0.24)
rkt.CriticalProps.append(crprops)

# ==== C7-C15 ====
db.addElement(rkt.Element()
              .withName("C7-C15")
              .withSymbol("C7-C15")
              .withMolarMass(0.14096))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C7-C15(l)")
              .withFormula("C7-C15")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C7-C15"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C7-C15(l)")
crprops.setTemperature(605.694, "K")
crprops.setPressure(21.749, "bar")
crprops.setAcentricFactor(0.618)
rkt.CriticalProps.append(crprops)

# ==== C16-C27 ====
db.addElement(rkt.Element()
              .withName("C16-C27")
              .withSymbol("C16-C27")
              .withMolarMass(0.28099))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C16-C27(l)")
              .withFormula("C16-C27")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C16-C27"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C16-C27(l)")
crprops.setTemperature(751.017, "K")
crprops.setPressure(16.541, "bar")
crprops.setAcentricFactor(0.957)
rkt.CriticalProps.append(crprops)

# ==== C28+ ====
db.addElement(rkt.Element()
              .withName("C28+")
              .withSymbol("C28+")
              .withMolarMass(0.51962))  # kg/mol

db.addSpecies(rkt.Species()
              .withName("C28+(l)")
              .withFormula("C28+")
              .withCharge(0.)
              .withElements(rkt.ElementalComposition([(db.elements().get("C28+"), 1.0)]))
              .withAggregateState(rkt.AggregateState.Liquid)
              .withStandardGibbsEnergy(0.))

# Specify critical properties
crprops = rkt.SubstanceCriticalProps("C28+(l)")
crprops.setTemperature(942.478, "K")
crprops.setPressure(16.418, "bar")
crprops.setAcentricFactor(1.268)
rkt.CriticalProps.append(crprops)

# Duplicate species to have gaseous and liquid phases
db.addSpecies(db.species("CO2(l)").clone().withName("CO2(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("CH4(l)").clone().withName("CH4(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("C2-C3(l)").clone().withName("C2-C3(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("C4-C6(l)").clone().withName("C4-C6(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("C7-C15(l)").clone().withName("C7-C15(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("C16-C27(l)").clone().withName("C16-C27(gas)").withAggregateState(rkt.AggregateState.Gas))
db.addSpecies(db.species("C28+(l)").clone().withName("C28+(gas)").withAggregateState(rkt.AggregateState.Gas))

db.addSpecies(db.species("CO2(l)").clone().withName("CO2(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("CH4(l)").clone().withName("CH4(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("C2-C3(l)").clone().withName("C2-C3(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("C4-C6(l)").clone().withName("C4-C6(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("C7-C15(l)").clone().withName("C7-C15(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("C16-C27(l)").clone().withName("C16-C27(oil)").withAggregateState(rkt.AggregateState.Liquid))
db.addSpecies(db.species("C28+(l)").clone().withName("C28+(oil)").withAggregateState(rkt.AggregateState.Liquid))

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

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

    iCO2    = index_substance("CO2")
    iCH4    = index_substance("CH4")
    iC2C3   = index_substance("C2-C3")
    iC4C6   = index_substance("C4-C6")
    iC7C15  = index_substance("C7-C15")
    iC16C27 = index_substance("C16-C27")
    iC28    = index_substance("C28+")

    size = len(substances)

    bip = rkt.CubicEOS.Bip(size)
    bip.k[iCO2, iCO2]    = bip.k[iCO2, iCO2]    = 0.000
    bip.k[iCO2, iCH4]    = bip.k[iCH4, iCO2]    = 0.055
    bip.k[iCO2, iC2C3]   = bip.k[iC2C3, iCO2]   = 0.055
    bip.k[iCO2, iC4C6]   = bip.k[iC4C6, iCO2]   = 0.055
    bip.k[iCO2, iC7C15]  = bip.k[iC7C15, iCO2]  = 0.105
    bip.k[iCO2, iC16C27] = bip.k[iC16C27, iCO2] = 0.105
    bip.k[iCO2, iC28]    = bip.k[iC28, iCO2]    = 0.105

    # Define the function that will be called whenever the BIP parameters need to be updated
    def evalfn(args):
        return bip

    return evalfn

# Set up phases in chemical system
gas = rkt.GaseousPhase("CO2(gas) CH4(gas) C2-C3(gas) C4-C6(gas) C7-C15(gas) C16-C27(gas) C28+(gas)")
gas.set(rkt.ActivityModelPengRobinson(bip_model))

liquid_1 = rkt.LiquidPhase("CO2(l) CH4(l) C2-C3(l) C4-C6(l) C7-C15(l) C16-C27(l) C28+(l)")
liquid_1.set(rkt.ActivityModelPengRobinson(bip_model))

liquid_2 = rkt.LiquidPhase("CO2(oil) CH4(oil) C2-C3(oil) C4-C6(oil) C7-C15(oil) C16-C27(oil) C28+(oil)")
liquid_2.set(rkt.ActivityModelPengRobinson(bip_model))

# Define chemical system
phases = rkt.Phases(db)
phases.add(gas)
phases.add(liquid_1)
phases.add(liquid_2)
system = rkt.ChemicalSystem(phases)


# Create a chemical state object
state = rkt.ChemicalState(system)
state.setTemperature(105., "F")
state.pressure(1295., "psi")

state.set("CO2(l)", 0.57529, "mol")
state.set("CH4(l)", 0.07305, "mol")
state.set("C2-C3(l)", 0.10297, "mol")
state.set("C4-C6(l)", 0.09088, "mol")
state.set("C7-C15(l)", 0.10509, "mol")
state.set("C16-C27(l)", 0.03792, "mol")
state.set("C28+(l)", 0.01480, "mol")

restrictions = rkt.EquilibriumRestrictions(system)
restrictions.cannotDecreaseBelow("CO2(l)", 0.01, "mol")
restrictions.cannotDecreaseBelow("CO2(oil)", 0.50, "mol")

res = rkt.equilibrate(state, restrictions)

assert res.succeeded()

print(state)  # State after equilibration

props = rkt.ChemicalProps(state)
for x, species in zip(props.speciesMoleFractions(), system.species()):
    print(f"{species.name():<20}: {x:.6f}")

legofernando1999 Apr 11, 2024
Author

Hi @allanleal, thank you for the explanation and the code example, it was very helpful.

Just one more question, how does Reaktoro handle the setting of species amounts? If I initialize the amounts for each species in each phase, does Reaktoro use this information in the equilibrate() method for initializing the independent variables (x)? Or do you take all of the species amounts set by the user and calculate the overall composition and then use some other procedure to initialize x?

Also, relevant to the initialization question, if I want to run a series of flash calculations where only temperature and/or pressure is changing, can I use a previously converged solution from equilibrate() to initialize my next call to equilibrate at a slightly different set of (T, P)?

allanleal Apr 12, 2024
Maintainer

how does Reaktoro handle the setting of species amounts?

You can set the amounts of each species individually using the method ChemicalState.set.

If I initialize the amounts for each species in each phase, does Reaktoro use this information in the equilibrate() method for initializing the independent variables (x)?

Yes, it will use this information as an initial guess for the calculation. The combination of set methods plus setting T and P produces an initial chemical state in disequilibrium, which you can equilibrate subsequently. The x notation I used before was for vector of mole fractions. In a chemical state, you set n, the vector of species amounts.

Or do you take all of the species amounts set by the user and calculate the overall composition and then use some other procedure to initialize x?

No, the initial amounts you set will be the initial guess.

Also, relevant to the initialization question, if I want to run a series of flash calculations where only temperature and/or pressure is changing, can I use a previously converged solution from equilibrate() to initialize my next call to equilibrate at a slightly different set of (T, P)?

Yes, you can use the previously computed chemical state (in equilibrium) as the initial guess for a subsequent calculation.

Let me know if you need further assistance.

legofernando1999 Apr 12, 2024
Author

I see, I don't think I need further assistance for now. Thank you very much for your time.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Modeling of a VLL system. #377

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 1 comment 5 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Modeling of a VLL system. #377

legofernando1999 Apr 5, 2024

Description

Example

What is your operating system?

Additional Information

Replies: 1 comment · 5 replies

allanleal Apr 9, 2024 Maintainer

legofernando1999 Apr 10, 2024 Author

allanleal Apr 11, 2024 Maintainer

legofernando1999 Apr 11, 2024 Author

allanleal Apr 12, 2024 Maintainer

legofernando1999 Apr 12, 2024 Author

legofernando1999
Apr 5, 2024

Replies: 1 comment 5 replies

allanleal
Apr 9, 2024
Maintainer

legofernando1999 Apr 10, 2024
Author

allanleal Apr 11, 2024
Maintainer

legofernando1999 Apr 11, 2024
Author

allanleal Apr 12, 2024
Maintainer

legofernando1999 Apr 12, 2024
Author