fix: prevent mass loss for supercritical fluids in liquid remover

Author: olelod

src/ecalc_neqsim_wrapper/fluid_service.py

@@ -511,6 +511,39 @@ def mass_rate_to_standard_rate(
         standard_density = self._get_standard_density(fluid_model)
         return mass_rate_kg_per_h * 24.0 / standard_density
 
+    _critical_point_cache: ClassVar[dict[tuple, tuple[float, float]]] = {}
+
+    def get_critical_point(
+        self,
+        fluid_model: FluidModel,
+    ) -> tuple[float, float]:
+        """Get the EoS-computed critical point for a fluid composition.
+
+        Uses NeqSim's criticalPointFlash() which solves for the true mixture
+        critical point using the equation of state. Results are cached by
+        (composition, eos_model) since the critical point is independent of
+        the stream's actual T and P.
+
+        Returns:
+            Tuple of (critical_temperature_kelvin, critical_pressure_bara)
+        """
+        composition = fluid_model.composition.normalized()
+        key = (_make_composition_key(composition), fluid_model.eos_model)
+
+        cached = self._critical_point_cache.get(key)
+        if cached is not None:
+            return cached
+
+        # Create a disposable fluid — critical_point() clones internally
+        disposable = NeqsimFluid.create_thermo_system(
+            composition=composition,
+            eos_model=fluid_model.eos_model,
+        )
+        result = disposable.critical_point()
+        self._critical_point_cache[key] = result
+        _logger.debug("Critical point for %s: Tc=%.2f K, Pc=%.2f bar", fluid_model.eos_model.name, *result)
+        return result
+
 
 def get_fluid_service_stats() -> dict[str, dict]:
     """Get cache statistics for the fluid service caches."""

src/ecalc_neqsim_wrapper/thermo.py

@@ -368,6 +368,20 @@ def clone_gas_phase(self) -> NeqsimFluid:
     def copy(self) -> NeqsimFluid:
         return NeqsimFluid(thermodynamic_system=self._thermodynamic_system.clone(), use_gerg=self._use_gerg)
 
+    def critical_point(self) -> tuple[float, float]:
+        """Compute the EoS critical point for this fluid's composition.
+
+        Uses a disposable clone — the original thermo system is not mutated.
+
+        Returns:
+            Tuple of (critical_temperature_kelvin, critical_pressure_bara)
+        """
+        ts = self._thermodynamic_system.clone()
+        neqsim_module = NeqsimService.instance().get_neqsim_module()
+        ops = neqsim_module.thermodynamicoperations.ThermodynamicOperations(ts)  # pyright: ignore[reportAttributeAccessIssue]
+        ops.criticalPointFlash()
+        return float(ts.getTemperature()), float(ts.getPressure())
+
     @Capturer.capture_return_values(  # type: ignore[misc]
         do_save_captured_content=False, output_directory=Path(os.getcwd()) / "captured_data" / "neqsim-ph"
     )

src/libecalc/process/fluid_stream/fluid_service.py

@@ -197,3 +197,24 @@ def mass_rate_to_standard_rate(
             Volumetric flow rate at standard conditions [Sm3/day]
         """
         ...
+
+    # === Critical Point ===
+
+    @abc.abstractmethod
+    def get_critical_point(
+        self,
+        fluid_model: FluidModel,
+    ) -> tuple[float, float]:
+        """Get the critical temperature and pressure for a fluid composition.
+
+        Uses the equation of state to compute the true mixture critical point.
+        Results should be cached by composition + EoS since the critical point
+        is independent of the stream's actual T and P.
+
+        Args:
+            fluid_model: The fluid model (composition + EoS)
+
+        Returns:
+            Tuple of (critical_temperature_kelvin, critical_pressure_bara)
+        """
+        ...

src/libecalc/process/process_pipeline/process_error.py

@@ -63,6 +63,19 @@ def __init__(
         super().__init__(f"Inlet stream contains liquid (vapor fraction: {vapor_fraction:.3f})")
 
 
+class NoGasPhaseError(ProcessError):
+    """The stream has no gas phase — liquid removal produces nothing to compress."""
+
+    def __init__(
+        self,
+        process_unit_id: ProcessUnitId,
+        vapor_fraction: float,
+    ):
+        self.process_unit_id = process_unit_id
+        self.vapor_fraction = vapor_fraction
+        super().__init__(f"No gas phase present (vapor fraction: {vapor_fraction:.6f})")
+
+
 class OfftakeExceedsInletError(ProcessError):
     def __init__(
         self,

src/libecalc/process/process_solver/solver.py

@@ -14,6 +14,7 @@
     CompressorSurgeError,
     InsufficientInletPressureError,
     LiquidAtInletError,
+    NoGasPhaseError,
     OfftakeExceedsInletError,
     ProcessError,
 )
@@ -101,6 +102,12 @@ class LiquidAtInletFailure(SolverFailure):
     vapor_fraction: float | None = None
 
 
+@dataclass
+class NoGasPhaseFailure(SolverFailure):
+    process_unit_id: ProcessUnitId | None = None
+    vapor_fraction: float | None = None
+
+
 @dataclass
 class InsufficientInletPressureFailure(SolverFailure):
     process_unit_id: ProcessUnitId | None = None
@@ -134,6 +141,8 @@ def process_error_to_failure(e: ProcessError) -> SolverFailure:
     """Map a ProcessError to the appropriate typed SolverFailure."""
     if isinstance(e, LiquidAtInletError):
         return LiquidAtInletFailure(process_unit_id=e.process_unit_id, vapor_fraction=e.vapor_fraction)
+    if isinstance(e, NoGasPhaseError):
+        return NoGasPhaseFailure(process_unit_id=e.process_unit_id, vapor_fraction=e.vapor_fraction)
     if isinstance(e, OfftakeExceedsInletError):
         return OfftakeExceedsInletFailure(
             process_unit_id=e.process_unit_id, available_rate=e.available_rate, offtake_rate=e.offtake_rate

src/libecalc/process/process_units/liquid_remover.py

@@ -1,10 +1,17 @@
+import logging
 from typing import Final
 
 from libecalc.process.fluid_stream.constants import ThermodynamicConstants
 from libecalc.process.fluid_stream.fluid_service import FluidService
 from libecalc.process.fluid_stream.fluid_stream import FluidStream
+from libecalc.process.process_pipeline.process_error import NoGasPhaseError
 from libecalc.process.process_pipeline.process_unit import ProcessUnit, ProcessUnitId
 
+logger = logging.getLogger(__name__)
+
+_PURE_VAPOR_THRESHOLD = ThermodynamicConstants.PURE_VAPOR_THRESHOLD  # 0.9999
+_PURE_LIQUID_THRESHOLD = 1.0 - _PURE_VAPOR_THRESHOLD  # 0.0001
+
 
 class LiquidRemover(ProcessUnit):
     def __init__(self, fluid_service: FluidService, process_unit_id: ProcessUnitId | None = None):
@@ -14,32 +21,53 @@ def __init__(self, fluid_service: FluidService, process_unit_id: ProcessUnitId |
     def get_id(self) -> ProcessUnitId:
         return self._id
 
-    def propagate_stream(self, inlet_stream: FluidStream) -> FluidStream:
+    def _is_supercritical(self, inlet_stream: FluidStream) -> bool:
+        """Check if the fluid is above its EoS-computed critical point.
+
+        Uses the fluid service's cached critical point calculation, which is
+        exact for both pure components and mixtures.
         """
-        Removes liquid from the fluid stream. The new stream's mass rate is scaled
-        down by the gas mass fraction so the dropped-out liquid isn't re-injected
-        into the gas phase.
+        tc, pc = self._fluid_service.get_critical_point(inlet_stream.fluid_model)
+        return inlet_stream.temperature_kelvin > tc and inlet_stream.pressure_bara > pc
 
-        The removed liquid (mass = inlet.mass_rate * (1 - gas_mass_fraction),
-        composition = inlet - new_fluid) is currently discarded. It could later
-        be exposed as a separate outlet stream — e.g. routed to an oil pump,
-        accounted for in emissions, or reported back to the user.
+    def propagate_stream(self, inlet_stream: FluidStream) -> FluidStream:
+        """Remove liquid from a two-phase fluid stream.
 
-        Args:
-            inlet_stream: The fluid stream to be scrubbed.
+        Liquid removal only makes sense when both phases coexist:
 
-        Returns:
-            FluidStream: A new FluidStream with liquid removed.
+        - vapor_fraction >= 0.9999: all gas, nothing to remove.
+        - vapor_fraction <= 0.0001 AND supercritical: mislabelled by flash,
+          pass through unchanged.
+        - vapor_fraction <= 0.0001 AND NOT supercritical: genuinely liquid,
+          raise NoGasPhaseError — liquid removal cannot produce gas.
+        - Otherwise: genuine two-phase, remove liquid and keep gas.
         """
-        if inlet_stream.vapor_fraction_molar < ThermodynamicConstants.PURE_VAPOR_THRESHOLD:
-            new_fluid = self._fluid_service.remove_liquid(inlet_stream.fluid)
-            inlet_molar_mass = inlet_stream.fluid.molar_mass
-            assert inlet_molar_mass > 0.0, (
-                f"Degenerate stream with non-positive molar mass ({inlet_molar_mass}) reached LiquidRemover — "
-                "this should have been caught at stream construction."
-            )
-            gas_mass_fraction = inlet_stream.vapor_fraction_molar * new_fluid.molar_mass / inlet_molar_mass
-            new_mass_rate = inlet_stream.mass_rate_kg_per_h * gas_mass_fraction
-            return inlet_stream.with_new_fluid(new_fluid).with_mass_rate(new_mass_rate)
-        else:
+        vf = inlet_stream.vapor_fraction_molar
+
+        if vf >= _PURE_VAPOR_THRESHOLD:
             return inlet_stream
+
+        if vf <= _PURE_LIQUID_THRESHOLD:
+            if self._is_supercritical(inlet_stream):
+                logger.debug(
+                    "LiquidRemover: skipping — supercritical fluid (T=%.1f K, P=%.1f bara, vf=%.6f)",
+                    inlet_stream.temperature_kelvin,
+                    inlet_stream.pressure_bara,
+                    vf,
+                )
+                return inlet_stream
+
+            raise NoGasPhaseError(
+                process_unit_id=self._id,
+                vapor_fraction=vf,
+            )
+
+        new_fluid = self._fluid_service.remove_liquid(inlet_stream.fluid)
+        inlet_molar_mass = inlet_stream.fluid.molar_mass
+        assert inlet_molar_mass > 0.0, (
+            f"Degenerate stream with non-positive molar mass ({inlet_molar_mass}) reached LiquidRemover — "
+            "this should have been caught at stream construction."
+        )
+        gas_mass_fraction = vf * new_fluid.molar_mass / inlet_molar_mass
+        new_mass_rate = inlet_stream.mass_rate_kg_per_h * gas_mass_fraction
+        return inlet_stream.with_new_fluid(new_fluid).with_mass_rate(new_mass_rate)

tests/libecalc/process/process_units/test_liquid_remover.py

@@ -1,11 +1,9 @@
-from unittest.mock import MagicMock
-
 import pytest
 
 from ecalc_neqsim_wrapper.thermo import STANDARD_PRESSURE_BARA, STANDARD_TEMPERATURE_KELVIN
 from libecalc.process.fluid_stream.fluid_model import EoSModel, FluidComposition, FluidModel
 from libecalc.process.fluid_stream.fluid_stream import FluidStream
-from libecalc.process.process_units.liquid_remover import LiquidRemover
+from libecalc.process.process_pipeline.process_error import NoGasPhaseError
 
 
 def test_liquid_remover_removes_liquid(fluid_service, liquid_remover_factory):
@@ -38,23 +36,11 @@ def test_liquid_remover_removes_liquid(fluid_service, liquid_remover_factory):
 
     assert inlet_stream.vapor_fraction_molar < 1.0
     assert outlet_stream.vapor_fraction_molar == 1.0
-
-    expected_gas_mass_fraction = (
-        inlet_stream.vapor_fraction_molar * outlet_stream.fluid.molar_mass / inlet_stream.fluid.molar_mass
-    )
-    expected_mass_rate = inlet_stream.mass_rate_kg_per_h * expected_gas_mass_fraction
     assert outlet_stream.mass_rate_kg_per_h < inlet_stream.mass_rate_kg_per_h
-    assert outlet_stream.mass_rate_kg_per_h == expected_mass_rate
 
 
 def test_liquid_remover_passthrough_when_no_liquid(fluid_service, liquid_remover_factory):
-    composition = FluidComposition(
-        nitrogen=3,
-        CO2=1,
-        methane=80,
-        ethane=10,
-        propane=6,
-    )
+    composition = FluidComposition(nitrogen=3, CO2=1, methane=80, ethane=10, propane=6)
     fluid_model = FluidModel(eos_model=EoSModel.SRK, composition=composition)
     fluid = fluid_service.create_fluid(
         fluid_model=fluid_model,
@@ -73,13 +59,117 @@ def test_liquid_remover_passthrough_when_no_liquid(fluid_service, liquid_remover
     assert outlet_stream.mass_rate_kg_per_h == inlet_stream.mass_rate_kg_per_h
 
 
-def test_liquid_remover_raises_on_non_positive_inlet_molar_mass():
+def test_liquid_remover_passthrough_supercritical_co2(fluid_service, liquid_remover_factory):
+    """Pure CO2 above critical point: NeqSim reports vapor_fraction=0,
+    but the EoS critical point check detects supercritical and prevents mass loss."""
+    composition = FluidComposition(CO2=1.0)
+    fluid_model = FluidModel(eos_model=EoSModel.SRK, composition=composition)
+
+    fluid = fluid_service.create_fluid(
+        fluid_model=fluid_model,
+        pressure_bara=350.0,
+        temperature_kelvin=308.15,
+    )
+    inlet_stream = FluidStream.from_standard_rate(
+        standard_rate_m3_per_day=100000,
+        fluid_model=fluid.fluid_model,
+        fluid_properties=fluid.properties,
+    )
+
+    # NeqSim mislabels supercritical CO2 as liquid
+    assert inlet_stream.vapor_fraction_molar <= 0.0001
+
+    remover = liquid_remover_factory()
+    outlet_stream = remover.propagate_stream(inlet_stream)
+
+    # Mass fully conserved
+    assert outlet_stream.mass_rate_kg_per_h == inlet_stream.mass_rate_kg_per_h
+
+
+def test_liquid_remover_passthrough_subcritical_co2_vapor(fluid_service, liquid_remover_factory):
+    """CO2 below critical pressure: should be all vapor, passes through."""
+    composition = FluidComposition(CO2=1.0)
+    fluid_model = FluidModel(eos_model=EoSModel.SRK, composition=composition)
+
+    fluid = fluid_service.create_fluid(
+        fluid_model=fluid_model,
+        pressure_bara=50.0,
+        temperature_kelvin=293.15,
+    )
+    inlet_stream = FluidStream.from_standard_rate(
+        standard_rate_m3_per_day=100000,
+        fluid_model=fluid.fluid_model,
+        fluid_properties=fluid.properties,
+    )
+
+    remover = liquid_remover_factory()
+    outlet_stream = remover.propagate_stream(inlet_stream)
+
+    # Pure CO2 at these conditions is single-phase, mass conserved
+    assert outlet_stream.mass_rate_kg_per_h == inlet_stream.mass_rate_kg_per_h
+
+
+def test_liquid_remover_raises_for_genuine_liquid(fluid_service, liquid_remover_factory):
+    """Water at ambient conditions is genuinely liquid (not supercritical).
+    The LiquidRemover raises NoGasPhaseError — no gas to extract."""
+    composition = FluidComposition(water=1.0)
+    fluid_model = FluidModel(eos_model=EoSModel.SRK, composition=composition)
+
+    fluid = fluid_service.create_fluid(
+        fluid_model=fluid_model,
+        pressure_bara=10.0,
+        temperature_kelvin=293.15,
+    )
+    inlet_stream = FluidStream.from_standard_rate(
+        standard_rate_m3_per_day=100000,
+        fluid_model=fluid.fluid_model,
+        fluid_properties=fluid.properties,
+    )
+
+    assert inlet_stream.vapor_fraction_molar <= 0.0001
+
+    remover = liquid_remover_factory()
+    with pytest.raises(NoGasPhaseError):
+        remover.propagate_stream(inlet_stream)
+
+
+def test_liquid_remover_passthrough_supercritical_co2_mixture(fluid_service, liquid_remover_factory):
+    """CO2-dominant mixture above its EoS-computed critical point.
+    Exercises the critical point calculation across multiple components."""
+    composition = FluidComposition(CO2=95.0, methane=5.0)
+    fluid_model = FluidModel(eos_model=EoSModel.SRK, composition=composition)
+
+    # EoS critical point ≈ 297 K, 72 bar
+    # At 310 K, 100 bar → above both → supercritical
+    fluid = fluid_service.create_fluid(
+        fluid_model=fluid_model,
+        pressure_bara=100.0,
+        temperature_kelvin=310.0,
+    )
+    inlet_stream = FluidStream.from_standard_rate(
+        standard_rate_m3_per_day=100000,
+        fluid_model=fluid.fluid_model,
+        fluid_properties=fluid.properties,
+    )
+
+    remover = liquid_remover_factory()
+    outlet_stream = remover.propagate_stream(inlet_stream)
+
+    assert outlet_stream.mass_rate_kg_per_h == inlet_stream.mass_rate_kg_per_h
+
+
+def test_liquid_remover_raises_on_non_positive_inlet_molar_mass(fluid_service, liquid_remover_factory):
+    """The assertion guarding against degenerate streams with zero molar mass is still active."""
+    from unittest.mock import MagicMock
+
+    from libecalc.process.process_units.liquid_remover import LiquidRemover
+
     inlet_stream = MagicMock()
     inlet_stream.vapor_fraction_molar = 0.5
     inlet_stream.fluid.molar_mass = 0.0
 
-    fluid_service = MagicMock()
-    remover = LiquidRemover(fluid_service=fluid_service)
+    mock_fluid_service = MagicMock()
+    remover = LiquidRemover(fluid_service=mock_fluid_service)
 
     with pytest.raises(AssertionError, match="non-positive molar mass"):
         remover.propagate_stream(inlet_stream)