search

equinor / neqsim-python

NeqSim is a library for calculation of fluid behavior, phase equilibrium and process simulation. This project is a Python interface to NeqSim.
https://equinor.github.io/neqsimhome/
Apache License 2.0
58 stars 15 forks source link

ERROR neqsim.thermo.phase.PhaseEos - Failed to solve for molarVolume within the iteration limit. #241

Closed sface-admin closed 3 months ago

sface-admin commented 9 months ago

I often get this error where composition has pseudo components. Is this particularly an error due to python envelope? The snip of example code below.

from neqsim.thermo.thermoTools import * import matplotlib.pyplot as plt

fluid1 = fluid("srk") # create a fluid using the SRK-EoS fluid1.setTemperature(80, "C") fluid1.setPressure(50, "bara") fluid1.addComponent("nitrogen", 0) fluid1.addComponent("CO2", 1.005) fluid1.addComponent("H2S", 0) fluid1.addComponent("methane", 24.62) fluid1.addComponent("ethane", 2.704) fluid1.addComponent("propane", 2.706) fluid1.addComponent("i-butane", 1.309) fluid1.addComponent("n-butane", 2.031) fluid1.addComponent("i-pentane", 2.047) fluid1.addComponent("n-pentane", 1.49) fluid1.addTBPfraction("C6", 3.176, 85.3/1000.0, 0.6649) fluid1.addTBPfraction("C7", 5.331, 93.2/1000.0, 0.7266) fluid1.addTBPfraction("C8", 7.796, 105.1/1000.0, 0.7607) fluid1.addTBPfraction("C9", 5.179, 119.4/1000.0, 0.7812) fluid1.addPlusFraction("C10+", 40.616, 268/1000.0, 0.865) fluid1.setMixingRule("classic")
fluid1.setMultiPhaseCheck(True) fluid1.useVolumeCorrection(True) TPflash(fluid1) fluid1.getInterphaseProperties().setInterfacialTensionModel("gas", "oil", "Parachor") fluid1.initProperties() fluid1.display() thermoOps = jNeqSim.thermodynamicOperations.ThermodynamicOperations(fluid1) thermoOps.calcPTphaseEnvelope() plt.plot(list(thermoOps.getOperation().get("dewT2")), list(thermoOps.getOperation().get("dewP2")), label="DewPoint") plt.plot(list(thermoOps.getOperation().get("bubT")), list(thermoOps.getOperation().get("bubP")), label="BubblePoint") plt.show()

EvenSol commented 9 months ago

Thanks for pointing at this.

It is a weakness in the phase envelope algorithm not accepting components with mole fraction 0. We will fix this in a later release, and mark this as a bug.

To run the script without any error, you can skip adding components where number of moles are 0. The script could be:

from neqsim.thermo.thermoTools import *
import matplotlib.pyplot as plt

fluid1 = fluid("srk") # create a fluid using the SRK-EoS
fluid1.setTemperature(80, "C")
fluid1.setPressure(50, "bara")
fluid1.addComponent("CO2", 1.005)
fluid1.addComponent("methane", 24.62)
fluid1.addComponent("ethane", 2.704)
fluid1.addComponent("propane", 2.706)
fluid1.addComponent("i-butane", 1.309)
fluid1.addComponent("n-butane", 2.031)
fluid1.addComponent("i-pentane", 2.047)
fluid1.addComponent("n-pentane", 1.49)
fluid1.addTBPfraction("C6", 3.176, 85.3/1000.0, 0.6649)
fluid1.addTBPfraction("C7", 5.331, 93.2/1000.0, 0.7266)
fluid1.addTBPfraction("C8", 7.796, 105.1/1000.0, 0.7607)
fluid1.addTBPfraction("C9", 5.179, 119.4/1000.0, 0.7812)
fluid1.addPlusFraction("C10+", 40.616, 268/1000.0, 0.865)
fluid1.setMixingRule("classic")
fluid1.setMultiPhaseCheck(True)
fluid1.useVolumeCorrection(True)
TPflash(fluid1)
fluid1.getInterphaseProperties().setInterfacialTensionModel("gas", "oil", "Parachor")
fluid1.initProperties()
fluid1.display()
thermoOps = jNeqSim.thermodynamicOperations.ThermodynamicOperations(fluid1)
thermoOps.calcPTphaseEnvelope()
plt.plot(list(thermoOps.getOperation().get("dewT")), list(thermoOps.getOperation().get("dewP")), label="DewPoint")
plt.plot(list(thermoOps.getOperation().get("bubT")), list(thermoOps.getOperation().get("bubP")), label="BubblePoint")
plt.show()
sface-admin commented 8 months ago

Thanks, that does solve for this composition. However I get similar error for another gas composition, where all components are non zero. Notice that error starts with addition of C7 pseudo onwards. See the script below,

fluid1 = fluid("srk")
fluid1.setTemperature(40, "C") fluid1.setPressure(50, "bara") fluid1.addComponent("nitrogen", 0.88) fluid1.addComponent("CO2", 5.7) fluid1.addComponent("methane", 86.89) fluid1.addComponent("ethane", 3.59) fluid1.addComponent("propane", 1.25) fluid1.addComponent("i-butane", 0.19) fluid1.addComponent("n-butane", 0.35) fluid1.addComponent("i-pentane", 0.12) fluid1.addComponent("n-pentane", 0.12) fluid1.addTBPfraction("C6", 0.15, 86/1000.0, 0.672) fluid1.addTBPfraction("C7", 0.2, 96/1000.0, 0.737) fluid1.addTBPfraction("C8", 0.22, 106/1000.0, 0.767) fluid1.addTBPfraction("C9", 0.13, 121/1000.0, 0.783) fluid1.addPlusFraction("C10+", 0.21, 172/1000.0, 0.818) fluid1.setMixingRule("classic") fluid1.setMultiPhaseCheck(True) fluid1.useVolumeCorrection(True) initPhysicalProperties() thermoOps = jNeqSim.thermodynamicOperations.ThermodynamicOperations(fluid1) thermoOps.calcPTphaseEnvelope() plt.plot(list(thermoOps.getOperation().get("dewT")), list(thermoOps.getOperation().get("dewP")), label="DewPoint") plt.plot(list(thermoOps.getOperation().get("bubT")), list(thermoOps.getOperation().get("bubP")), label="BubblePoint") plt.show()

EvenSol commented 3 months ago

All this seems to work in the new version and using the method:

thermoOps = jNeqSim.thermodynamicOperations.ThermodynamicOperations(fluid1) thermoOps.calcPTphaseEnvelope2()

See the tests in the notebook: https://colab.research.google.com/drive/1tuSXN5Myk97SSxxJJ21vo4G8lwILT9K2#scrollTo=iKZ_5EKuf09H