solubilityccs.py

from acidwatch_api.models.base import (
    BaseAdapter,
    BaseParameters,
    Parameter,
    RunResult,
)
from acidwatch_api.models.datamodel import TextResult, Phase

from solubilityccs import Fluid, ModelResults  # type: ignore
from solubilityccs.neqsim_functions import get_co2_parameters  # type: ignore

DESCRIPTION: str = """Solubility model detects acid formation risks in CO2 streams. 
    
It uses the SRK-CPA (Soave-Redlich-Kwong Cubic Plus Association) equation of state to calculate fugacity coefficients and activity models to determine component activities in multiphase systems.

The model currently supports the following chemical systems:
CO₂-water (binary system)
CO₂-water-H₂SO₄ (ternary system with sulfuric acid)
CO₂-water-HNO₃ (ternary system with nitric acid)"""


class SolubilityCCSParameters(BaseParameters):
    flow_rate: float = Parameter(
        10,
        label="Flow rate",
        min=0.01,
        max=100,
        unit="Mt/year",
        description="Flow rate in Mt/year",
    )


class SolubilityCCSAdapter(BaseAdapter):
    model_id = "solubilityccs"
    display_name = "Solubility CCS"
    description = DESCRIPTION
    valid_substances = ["H2O", "H2SO4", "HNO3"]
    parameters: SolubilityCCSParameters
    category = "PhaseEquilibrium"

    async def run(self) -> RunResult:
        # Get concentrations (mole fractions)
        h2o = self.concentrations.get("H2O", 0.0)
        h2so4 = self.concentrations.get("H2SO4", 0.0)
        hno3 = self.concentrations.get("HNO3", 0.0)
        # Conditions.temperature is in Celsius; solubilityccs expects Kelvin.
        temp = self.conditions.temperature + 273
        pres = self.conditions.pressure
        flow_rate = self.parameters.flow_rate

        co2 = 1e6 - (h2o + h2so4 + hno3)
        fluid = Fluid()
        fluid.add_component("CO2", co2)

        if h2o > 0:
            fluid.add_component("H2O", h2o)
        if h2so4 > 0:
            fluid.add_component("H2SO4", h2so4)
        elif hno3 > 0:
            fluid.add_component("HNO3", hno3)
        fluid.set_temperature(temp)
        fluid.set_pressure(pres)
        fluid.set_flow_rate(flow_rate * 1e6 * 1000 / (365 * 24), "kg/hr")
        fluid.calc_vapour_pressure()
        fluid.flash_activity()

        co2_properties = get_co2_parameters(pres, temp)
        results_obj = ModelResults(fluid, co2_properties=co2_properties)
        table = results_obj.generate_table()

        phases = self._extract_phases(fluid)

        return phases, TextResult(data=table, label="Solubility Output")

    def _extract_phases(self, fluid: Fluid) -> list[Phase]:
        co2_rich_phase = fluid.phases[0]
        co2_rich_fraction = fluid.betta

        co2_rich_concs: dict[str, float | int] = {}
        for component, fraction in zip(
            co2_rich_phase.components, co2_rich_phase.fractions
        ):
            if component != "CO2":
                co2_rich_concs[component] = fraction * 1e6

        phases = [
            Phase(
                kind="co2-rich",
                fraction=co2_rich_fraction,
                concentrations=co2_rich_concs,
            )
        ]

        if co2_rich_fraction < 1.0 and len(fluid.phases) > 1:
            liquid_phase = fluid.phases[1]
            aqueous_concs: dict[str, float | int] = {}
            for component, fraction in zip(
                liquid_phase.components, liquid_phase.fractions
            ):
                if component != "CO2":
                    aqueous_concs[component] = fraction * 1e6

            phases.append(
                Phase(
                    kind="aqueous",
                    fraction=1.0 - co2_rich_fraction,
                    concentrations=aqueous_concs,
                )
            )

        return phases