📚 Documentation

VLE 🌡️

VLE

Bases: Equilibria

Model for vapor-liquid equilibrium (VLE) calculations for multi-component systems using the Raoult's law, modified Raoult's law and more.

Source code in pyThermoFlash/docs/vle.py

class VLE(Equilibria):
    '''
    Model for vapor-liquid equilibrium (VLE) calculations for multi-component systems using the Raoult's law, modified Raoult's law and more.
    '''

    def __init__(self,
                 components: List[str],
                 model_source: Optional[Dict] = None,
                 **kwargs):
        '''
        Initialize the Partition class.
        '''
        self.__model_source = model_source
        self.__datasource = {
        } if model_source is None else model_source['datasource']
        self.__equationsource = {
        } if model_source is None else model_source['equationsource']

        # NOTE: init class

        Equilibria.__init__(self, components, model_source)

    @property
    def model_source(self) -> Dict:
        '''
        Get the model source property.

        Returns
        -------
        dict
            The model source dictionary.
        '''
        # NOTE: check if model source is valid
        if self.__model_source is None:
            return {}
        return self.__model_source

    @property
    def datasource(self) -> Dict:
        '''
        Get the datasource property.

        Returns
        -------
        dict
            The datasource dictionary.
        '''
        # NOTE: check if model source is valid
        if self.__datasource is None:
            return {}
        return self.__datasource

    @property
    def equationsource(self) -> Dict:
        '''
        Get the equationsource property.

        Returns
        -------
        dict
            The equationsource dictionary.
        '''
        # NOTE: check if model source is valid
        if self.__equationsource is None:
            return {}
        return self.__equationsource

    def __set_models(self,
                     equilibrium_model: str,
                     fugacity_model: Optional[str],
                     activity_model: Optional[str]):
        '''
        Set the models for the calculation.

        Parameters
        ----------
        equilibrium_model : str
            The equilibrium model to use for the calculation.
        fugacity_model : str, optional
            The fugacity model to use for the calculation.
        activity_model : str, optional
            The activity coefficient model to use for the calculation.

        Returns
        -------
        dict
            A dictionary containing the models used for the calculation.
            - equilibrium_model : str
                The equilibrium model used for the calculation.
            - fugacity_model : str, optional
                The fugacity model used for the calculation.
            - activity_model : str, optional
                The activity coefficient model used for the calculation.
        '''
        try:
            # NOTE: set based on equilibrium model
            if equilibrium_model == 'raoult':
                fugacity_model_ = None
                activity_model_ = None
            elif equilibrium_model == 'modified-raoult':
                fugacity_model_ = None
                activity_model_ = activity_model
            elif equilibrium_model == 'fugacity-ratio':
                fugacity_model_ = fugacity_model
                activity_model_ = None
            else:
                raise ValueError(
                    f"Invalid equilibrium model: {equilibrium_model}.")

            # res
            return {
                "equilibrium_model": equilibrium_model,
                "fugacity_model": fugacity_model_,
                "activity_model": activity_model_
            }

        except Exception as e:
            raise Exception(
                f"Error in setting models: {e}")

    def bubble_pressure(self,
                        inputs: Dict[str, Any],
                        equilibrium_model: Literal[
                            'raoult', 'modified-raoult'
                        ] = 'raoult',
                        fugacity_model: Optional[Literal[
                            'vdW', 'PR', 'RK', 'SRK'
                        ]] = None,
                        activity_model: Optional[Literal[
                            'NRTL', 'UNIQUAC'
                        ]] = None,
                        message: Optional[str] = None,
                        **kwargs):
        '''
        The Bubble-Pressure (BP) calculation determines the pressure at which the first bubble of vapor forms when a liquid mixture is heated at a constant temperature. It is used to find the pressure for a given temperature at which the liquid will begin to vaporize.

        Parameters
        ----------
        inputs : dict
            Dictionary containing the input parameters for the calculation.
            - mole_fraction : dict
                Dictionary of component names and their respective mole fractions.
            - temperature : float
                Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
        equilibrium_model : str, optional
            The equilibrium model to use for the calculation. Default is 'raoult'.
        fugacity_model : str, optional
            The fugacity model to use for the calculation. Default is 'SRK'.
        activity_model : str, optional
            The activity coefficient model to use for the calculation. Default is 'NRTL'.
        message : str, optional
            Message to display during the calculation. Default is None.
        **kwargs : dict, optional
            Additional parameters for the model.
            - activity_inputs : dict, model parameters.
            - activity_coefficients : dict, activity coefficients for each component.

        Returns
        -------
        res : dict
            Dictionary containing the results of the calculation.
            - bubble_pressure: bubble pressure [Pa]
            - temperature: temperature [K]
            - feed_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]
            - equilibrium_model: equilibrium model used for the calculation
            - fugacity_model: fugacity model used for the calculation
            - activity_model: activity coefficient model used for the calculation
            - components: list of components used in the calculation
            - message: message displayed during the calculation
            - computation_time: computation time [s]

        Notes
        -----
        - Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
        - Mole fractions must be between 0 and 1.
        - The function will raise a ValueError if the inputs are not valid.
        - The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
        - The function will return the bubble pressure in Pascals.
        - The activity coefficient for each component can be directly defined in the activity_inputs dictionary.


        For each component, the model parameters can be provided for activity models like NRTL or UNIQUAC:

        ```python
        # activity inputs
        activity_inputs = {
            'alpha': alpha,
            'a_ij': a_ij,
            'b_ij': b_ij,
            'c_ij': c_ij,
            'd_ij': d_ij
        }
        ```

        Or the activity coefficients can be provided directly:

        # calculated activity coefficients (for bubble-pressure calculation)
        ```python
        activity_coefficients = {
            'benzene': 1.0,
            'toluene': 1.1
        }
        ```
        '''
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check inputs
            # check inputs
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")

            # check keys
            if 'mole_fraction' not in inputs or 'temperature' not in inputs:
                raise ValueError(
                    "Inputs must contain 'mole_fraction' and 'temperature' keys.")

            # NOTE: set inputs - dict
            mole_fractions = inputs['mole_fraction']
            # NOTE: temperature [K]
            temperature = pycuc.convert_from_to(
                inputs['temperature'][0], inputs['temperature'][1], 'K')

            # check mole fractions
            if not isinstance(mole_fractions, dict):
                raise ValueError("Mole fractions must be a dictionary.")

            if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be numeric.")

            if not all(0 <= v <= 1 for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be between 0 and 1.")

            # check temperature
            if not isinstance(temperature, (int, float)):
                raise ValueError("Temperature must be numeric.")

            # NOTE: components
            components = self.components

            # mole fractions based on components id
            mole_fractions = [mole_fractions[component]
                              for component in components]

            # SECTION: extract source
            Source_ = Source(self.model_source)

            # NOTE: vapor pressure source
            VaPr_comp = {}

            # looping through components
            for component in components:
                # NOTE: equation source
                # antoine equations [Pa]
                VaPr_eq = Source_.eq_extractor(component, 'VaPr')

                # NOTE: args
                VaPr_args = VaPr_eq.args
                # check args (SI)
                VaPr_args_required = Source_.check_args(component, VaPr_args)

                # build args
                _VaPr_args = Source_.build_args(component, VaPr_args_required)

                # NOTE: update P and T
                _VaPr_args['T'] = temperature

                # NOTE: execute
                _VaPr_res = VaPr_eq.cal(**_VaPr_args)
                # extract
                _VaPr_value = _VaPr_res['value']
                _VaPr_unit = _VaPr_res['unit']
                # unit conversion
                # NOTE: unit conversion
                _unit_block = f"{_VaPr_unit} => Pa"
                _VaPr = pycuc.to(_VaPr_value, _unit_block)
                # set
                VaPr_comp[component] = {
                    "value": _VaPr,
                    "unit": "Pa"}

            # SECTION: set models
            # check equilibrium model
            res_ = self.__set_models(
                equilibrium_model=equilibrium_model,
                fugacity_model=fugacity_model,
                activity_model=activity_model
            )

            # NOTE: parameters
            params = {
                "components": components,
                "mole_fraction": np.array(mole_fractions),
                "temperature": {
                    "value": temperature,
                    "unit": "K"
                },
                "vapor_pressure": VaPr_comp,
                "equilibrium_model": equilibrium_model,
                "fugacity_model": res_['fugacity_model'],
                "activity_model": res_['activity_model'],
            }

            # res
            res = self._BP(params, **kwargs)

            # NOTE: set message
            message = message if message is not None else "Bubble Pressure Calculation"
            # add
            res['message'] = message

            # NOTE: components
            res['components'] = components

            # NOTE: models
            res["equilibrium_model"] = equilibrium_model
            res["fugacity_model"] = res_['fugacity_model']
            res["activity_model"] = res_['activity_model']

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # returns
            return res
        except Exception as e:
            raise ValueError(
                f"Error in bubble_pressure calculation: {e}")

    def dew_pressure(self,
                     inputs: Dict[str, Any],
                     equilibrium_model: Literal[
                         'raoult', 'modified-raoult'
                     ] = 'raoult',
                     fugacity_model: Optional[Literal[
                         'vdW', 'PR', 'RK', 'SRK'
                     ]] = None,
                     activity_model: Optional[Literal[
                         'NRTL', 'UNIQUAC'
                     ]] = None,
                     message: Optional[str] = None,
                     **kwargs):
        '''
        The dew-point pressure (DP) calculation determines the pressure at which the first drop of liquid condenses when a vapor mixture is cooled at a constant temperature. It is used to find the pressure at which vapor will begin to condense.

        Parameters
        ----------
        inputs : dict
            Dictionary containing the input parameters for the calculation.
            - mole_fraction : dict
                Dictionary of component names and their respective mole fractions.
            - temperature : float
                Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
        equilibrium_model : str, optional
            The equilibrium model to use for the calculation. Default is 'raoult'.
        fugacity_model : str, optional
            The fugacity model to use for the calculation. Default is 'SRK'.
        activity_model : str, optional
            The activity coefficient model to use for the calculation. Default is 'NRTL'.
        message : str, optional
            Message to display during the calculation. Default is None.
        **kwargs : dict, optional
            Additional parameters for the model.
            - activity_inputs : dict, model parameters.
            - activity_coefficients : dict, activity coefficients for each component.

        Returns
        -------
        res : dict
            Dictionary containing the results of the calculation.
            - dew_pressure: dew pressure [Pa]
            - temperature: temperature [K]
            - feed_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]
            - equilibrium_model: equilibrium model used for the calculation
            - fugacity_model: fugacity model used for the calculation
            - activity_model: activity coefficient model used for the calculation
            - components: list of components used in the calculation
            - message: message displayed during the calculation
            - computation_time: computation time [s]

        Notes
        -----
        - Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
        - Mole fractions must be between 0 and 1.
        - The function will raise a ValueError if the inputs are not valid.
        - The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
        - The function will return the dew pressure in Pascals.
        '''
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check inputs
            # check inputs
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")

            # check keys
            if 'mole_fraction' not in inputs or 'temperature' not in inputs:
                raise ValueError(
                    "Inputs must contain 'mole_fraction' and 'temperature' keys.")

            # NOTE: set inputs - dict
            mole_fractions = inputs['mole_fraction']
            # NOTE: temperature [K]
            temperature = pycuc.convert_from_to(
                inputs['temperature'][0], inputs['temperature'][1], 'K')

            # check mole fractions
            if not isinstance(mole_fractions, dict):
                raise ValueError("Mole fractions must be a dictionary.")

            if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be numeric.")

            if not all(0 <= v <= 1 for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be between 0 and 1.")

            # check temperature
            if not isinstance(temperature, (int, float)):
                raise ValueError("Temperature must be numeric.")

            # NOTE: components
            components = self.components

            # mole fractions based on components id
            mole_fractions = [mole_fractions[component]
                              for component in components]

            # SECTION: extract source
            Source_ = Source(self.model_source)

            # NOTE: vapor pressure source
            VaPr_comp = {}

            # looping through components
            for component in components:
                # NOTE: equation source
                # antoine equations [Pa]
                VaPr_eq = Source_.eq_extractor(component, 'VaPr')

                # NOTE: args
                VaPr_args = VaPr_eq.args
                # check args (SI)
                VaPr_args_required = Source_.check_args(component, VaPr_args)

                # build args
                _VaPr_args = Source_.build_args(component, VaPr_args_required)

                # NOTE: update P and T
                _VaPr_args['T'] = temperature

                # NOTE: execute
                _VaPr_res = VaPr_eq.cal(**_VaPr_args)
                # extract
                _VaPr_value = _VaPr_res['value']
                _VaPr_unit = _VaPr_res['unit']
                # unit conversion
                # NOTE: unit conversion
                _unit_block = f"{_VaPr_unit} => Pa"
                _VaPr = pycuc.to(_VaPr_value, _unit_block)
                # set
                VaPr_comp[component] = {
                    "value": _VaPr,
                    "unit": "Pa",
                    "equation": VaPr_eq,
                    "args": _VaPr_args,
                    "return": VaPr_eq.returns
                }

            # SECTION: set models
            # check equilibrium model
            res_ = self.__set_models(
                equilibrium_model=equilibrium_model,
                fugacity_model=fugacity_model,
                activity_model=activity_model
            )

            # NOTE: parameters
            params = {
                "components": components,
                "mole_fraction": np.array(mole_fractions),
                "temperature": {
                    "value": temperature,
                    "unit": "K"
                },
                "vapor_pressure": VaPr_comp,
                "equilibrium_model": equilibrium_model,
                "fugacity_model": res_['fugacity_model'],
                "activity_model": res_['activity_model'],
            }

            # res
            res = self._DP(params, **kwargs)

            # NOTE: set message
            message = message if message is not None else "Dew Pressure Calculation"
            # add
            res['message'] = message

            # NOTE: components
            res['components'] = components

            # NOTE: models
            res["equilibrium_model"] = equilibrium_model
            res["fugacity_model"] = res_['fugacity_model']
            res["activity_model"] = res_['activity_model']

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # returns
            return res
        except Exception as e:
            raise ValueError(
                f"Error in bubble_pressure calculation: {e}")

    def bubble_temperature(self,
                           inputs: Dict[str, Any],
                           equilibrium_model: Literal[
                               'raoult', 'modified-raoult'
                           ] = 'raoult',
                           fugacity_model: Literal[
                               'vdW', 'PR', 'RK', 'SRK'
                           ] = 'SRK',
                           activity_model: Literal[
                               'NRTL', 'UNIQUAC']
                           = 'NRTL',
                           solver_method: Literal[
                               'root', 'least-squares', 'fsolve'
                           ] = 'root',
                           message: Optional[str] = None,
                           **kwargs):
        '''
        The `bubble-point temperature` (BT) calculation determines the temperature at which the first bubble of vapor forms when a liquid mixture is heated at a constant pressure. It helps identify the temperature at which the liquid will start vaporizing.

        Parameters
        ----------
        inputs : dict
            Dictionary containing the input parameters for the calculation.
            - mole_fraction : dict
                Dictionary of component names and their respective mole fractions.
            - temperature : float
                Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
        equilibrium_model : str, optional
            The equilibrium model to use for the calculation. Default is 'raoult'.
        fugacity_model : str, optional
            The fugacity model to use for the calculation. Default is 'SRK'.
        activity_model : str, optional
            The activity coefficient model to use for the calculation. Default is 'NRTL'.
        solver_method : str, optional
            The solver method to use for the calculation. Default is 'root'.
            - 'root' : Root-finding algorithm.
            - 'least-squares' : Least-squares optimization algorithm.
            - 'fsolve' : SciPy's fsolve function.
        message : str, optional
            Message to display during the calculation. Default is None.
        **kwargs : dict, optional
            Additional parameters for the model.
            - `guess_temperature` : float
                Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

        Returns
        -------
        res : dict
            Dictionary containing the results of the calculation.
            - bubble_temperature: bubble temperature [K]
            - pressure: pressure [Pa]
            - feed_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]
            - equilibrium_model: equilibrium model used for the calculation
            - fugacity_model: fugacity model used for the calculation
            - activity_model: activity coefficient model used for the calculation
            - components: list of components used in the calculation
            - message: message displayed during the calculation
            - solver_method: solver method used for the calculation
            - computation_time: computation time [s]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Pressure (P) and mole fraction of the components in the liquid phase (xi).
        then zi = xi.
        - Computed Information: Temperature (T) and mole fraction in the vapor phase (yi).
        '''
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check inputs
            # check inputs
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")

            # check keys
            if 'mole_fraction' not in inputs or 'pressure' not in inputs:
                raise ValueError(
                    "Inputs must contain 'mole_fraction' and 'temperature' keys.")

            # NOTE: set inputs - dict
            mole_fractions = inputs['mole_fraction']
            # NOTE: pressure [Pa]
            pressure = pycuc.convert_from_to(
                inputs['pressure'][0], inputs['pressure'][1], 'Pa')

            # check mole fractions
            if not isinstance(mole_fractions, dict):
                raise ValueError("Mole fractions must be a dictionary.")

            if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be numeric.")

            if not all(0 <= v <= 1 for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be between 0 and 1.")

            # check temperature
            if not isinstance(pressure, (int, float)):
                raise ValueError("pressure must be numeric.")

            # NOTE: components
            components = self.components

            # mole fractions based on components id
            mole_fractions = [mole_fractions[component]
                              for component in components]

            # SECTION: check for activity model

            # SECTION: extract source
            Source_ = Source(self.model_source)

            # NOTE: vapor pressure source
            VaPr_comp = {}

            # looping through components
            for component in components:
                # NOTE: equation source
                VaPr_eq = None
                # antoine equations [Pa]
                VaPr_eq = Source_.eq_extractor(component, 'VaPr')

                # NOTE: args
                VaPr_args = VaPr_eq.args
                # check args (SI)
                VaPr_args_required = Source_.check_args(component, VaPr_args)

                # build args
                _VaPr_args = Source_.build_args(component, VaPr_args_required)

                # NOTE: update P and T
                _VaPr_args['T'] = None

                # set
                VaPr_comp[component] = {
                    "value": VaPr_eq,
                    "args": _VaPr_args,
                    "return": VaPr_eq.returns
                }

            # SECTION: set models
            # check equilibrium model
            res_ = self.__set_models(
                equilibrium_model=equilibrium_model,
                fugacity_model=fugacity_model,
                activity_model=activity_model
            )

            # NOTE: parameters
            params = {
                "components": components,
                "mole_fraction": np.array(mole_fractions),
                "pressure": {
                    "value": pressure,
                    "unit": "Pa"
                },
                "vapor_pressure": VaPr_comp,
                "equilibrium_model": equilibrium_model,
                "fugacity_model": res_['fugacity_model'],
                "activity_model": res_['activity_model'],
                "solver_method": solver_method
            }

            # res
            res = self._BT(params, **kwargs)

            # NOTE: set message
            message = message if message is not None else "Bubble Temperature Calculation"
            # add
            res['message'] = message

            # NOTE: components
            res['components'] = components

            # NOTE: models
            res["equilibrium_model"] = equilibrium_model
            res["fugacity_model"] = res_['fugacity_model']
            res["activity_model"] = res_['activity_model']
            res["solver_method"] = solver_method

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # returns
            return res
        except Exception as e:
            raise Exception(
                f"Error in bubble_temperature calculation: {e}")

    def dew_temperature(self,
                        inputs: Dict[str, Any],
                        equilibrium_model: Literal[
                            'raoult', 'modified-raoult'
                        ] = 'raoult',
                        fugacity_model: Literal[
                            'vdW', 'PR', 'RK', 'SRK'
                        ] = 'SRK',
                        activity_model: Literal[
                            'NRTL', 'UNIQUAC']
                        = 'NRTL',
                        solver_method: Literal[
                            'root', 'least-squares', 'fsolve'
                        ] = 'least-squares',
                        message: Optional[str] = None,
                        **kwargs):
        '''
        The `dew-point temperature` (DT) calculation determines the temperature at which the first drop of liquid condenses when a vapor mixture is cooled at a constant pressure. It identifies the temperature at which vapor will start to condense.

        Parameters
        ----------
        inputs : dict
            Dictionary containing the input parameters for the calculation.
            - mole_fraction : dict
                Dictionary of component names and their respective mole fractions.
            - temperature : float
                Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
        equilibrium_model : str, optional
            The equilibrium model to use for the calculation. Default is 'raoult'.
        fugacity_model : str, optional
            The fugacity model to use for the calculation. Default is 'SRK'.
        activity_model : str, optional
            The activity coefficient model to use for the calculation. Default is 'NRTL'.
        solver_method : str, optional
            The solver method to use for the calculation. Default is 'root'.
            - 'root' : Root-finding algorithm.
            - 'least-squares' : Least-squares optimization algorithm.
            - 'fsolve' : SciPy's fsolve function.
        message : str, optional
            Message to display during the calculation. Default is None.
        **kwargs : dict, optional
            Additional parameters for the model.
            - `guess_temperature` : float
                Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

        Returns
        -------
        res : dict
            Dictionary containing the results of the calculation.
            - dew_temperature: dew temperature [K]
            - pressure: pressure [Pa]
            - feed_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]
            - equilibrium_model: equilibrium model used for the calculation
            - fugacity_model: fugacity model used for the calculation
            - activity_model: activity coefficient model used for the calculation
            - components: list of components used in the calculation
            - message: message displayed during the calculation
            - solver_method: solver method used for the calculation
            - computation_time: computation time [s]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi).
        then zi = yi.
        - Computed Information: Temperature (T) and mole fraction in the liquid phase (xi).
        '''
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check inputs
            # check inputs
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")

            # check keys
            if 'mole_fraction' not in inputs or 'pressure' not in inputs:
                raise ValueError(
                    "Inputs must contain 'mole_fraction' and 'temperature' keys.")

            # NOTE: set inputs - dict
            mole_fractions = inputs['mole_fraction']
            # NOTE: pressure [Pa]
            # set
            pressure_value = float(inputs['pressure'][0])
            pressure_unit = str(inputs['pressure'][1])

            pressure = pycuc.convert_from_to(
                pressure_value, pressure_unit, 'Pa')

            # check mole fractions
            if not isinstance(mole_fractions, dict):
                raise ValueError("Mole fractions must be a dictionary.")

            if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be numeric.")

            if not all(0 <= v <= 1 for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be between 0 and 1.")

            # check temperature
            if not isinstance(pressure, (int, float)):
                raise ValueError("pressure must be numeric.")

            # NOTE: components
            components = self.components

            # mole fractions based on components id
            mole_fractions = [mole_fractions[component]
                              for component in components]

            # SECTION: extract source
            Source_ = Source(self.model_source)

            # NOTE: vapor pressure source
            VaPr_comp = {}

            # looping through components
            for component in components:
                # NOTE: equation source
                VaPr_eq = None
                # antoine equations [Pa]
                VaPr_eq = Source_.eq_extractor(component, 'VaPr')

                # NOTE: args
                VaPr_args = VaPr_eq.args
                # check args (SI)
                VaPr_args_required = Source_.check_args(component, VaPr_args)

                # build args
                _VaPr_args = Source_.build_args(component, VaPr_args_required)

                # NOTE: update P and T
                _VaPr_args['T'] = None

                # set
                VaPr_comp[component] = {
                    "value": VaPr_eq,
                    "args": _VaPr_args,
                    "return": VaPr_eq.returns
                }

            # SECTION: set models
            # check equilibrium model
            res_ = self.__set_models(
                equilibrium_model=equilibrium_model,
                fugacity_model=fugacity_model,
                activity_model=activity_model
            )

            # NOTE: parameters
            params = {
                "components": components,
                "mole_fraction": np.array(mole_fractions),
                "pressure": {
                    "value": pressure,
                    "unit": "Pa"
                },
                "vapor_pressure": VaPr_comp,
                "equilibrium_model": equilibrium_model,
                "fugacity_model": res_['fugacity_model'],
                "activity_model": res_['activity_model'],
                "solver_method": solver_method
            }

            # res
            res = self._DT(params, **kwargs)

            # NOTE: set message
            message = message if message is not None else "Dew Temperature Calculation"
            # add
            res['message'] = message

            # NOTE: components
            res['components'] = components

            # NOTE: models
            res["equilibrium_model"] = equilibrium_model
            res["fugacity_model"] = res_['fugacity_model']
            res["activity_model"] = res_['activity_model']
            res["solver_method"] = solver_method

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # returns
            return res
        except Exception as e:
            raise Exception(
                f"Error in dew_temperature calculation: {e}")

    def flash_isothermal(self,
                         inputs: Dict[str, Any],
                         equilibrium_model: Literal[
                             'raoult', 'modified-raoult'
                         ] = 'raoult',
                         fugacity_model: Literal[
                             'vdW', 'PR', 'RK', 'SRK'
                         ] = 'SRK',
                         activity_model: Literal[
                             'NRTL', 'UNIQUAC']
                         = 'NRTL',
                         solver_method: Literal[
                             'least_squares', 'minimize'
                         ] = 'least_squares',
                         flash_checker: bool = False,
                         message: Optional[str] = None,
                         **kwargs):
        '''
        The `isothermal-flash` (IFL) calculation This calculation determines the vapor and liquid phase compositions and amounts at a specified temperature and pressure. The system is "flashed" isothermally, meaning the temperature is kept constant while the phase behavior is calculated for a mixture.

        Parameters
        ----------
        inputs : dict
            Dictionary containing the input parameters for the calculation.
            - mole_fraction : dict
                Dictionary of component names and their respective mole fractions.
            - temperature : float
                Flash temperature (in Kelvin, Celsius, Fahrenheit).
            - pressure : float
                Flash pressure (in Pascal, bar, atm).
        equilibrium_model : str, optional
            the equilibrium model to use for the calculation. Default is 'raoult'.
        fugacity_model : str, optional
            The fugacity model to use for the calculation. Default is 'SRK'.
        activity_model : str, optional
            The activity coefficient model to use for the calculation. Default is 'NRTL'.
        solver_method : str, optional
            The solver method to use for the calculation. Default is 'least_squares'.
            - 'minimize' : Minimize the objective function.
            - 'least-squares' : Least-squares optimization algorithm.
        flash_checker : bool, optional
            If True, check if the system is in a two-phase region. Default is False.
        message : str, optional
            Message to display during the calculation. Default is None.
        **kwargs : dict, optional
            Additional parameters for the model.
            - `guess_V_F_ratio`: initial guess for the vapor-to-liquid ratio (V/F), default is 0.5

        Returns
        -------
        res : dict
            Dictionary containing the results of the calculation.
            - vapor_to_liquid_ratio: V/F ratio [dimensionless]
            - liquid_to_vapor_ratio: L/F ratio [dimensionless]
            - feed_mole_fraction: liquid mole fraction (zi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - vapor_pressure: vapor pressure [Pa]
            - K_ratio: K-ratio [dimensionless]
            - temperature: temperature [K]
            - pressure: pressure [Pa]
            - equilibrium_model: equilibrium model used for the calculation
            - fugacity_model: fugacity model used for the calculation
            - activity_model: activity coefficient model used for the calculation
            - components: list of components used in the calculation
            - message: message displayed during the calculation
            - solver_method: solver method used for the calculation
            - flash_checker: flash checker used for the calculation
            - flash_checker_res: flash checker result
            - computation_time: computation time [s]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Temperature (T), Pressure (P), and feed mole fraction of the components (zi).
        - Computed Information: Vapor and liquid phase compositions (yi and xi), and the vapor-to-liquid ratio (V/F).

        Flash occurs when the bubble pressure is greater than the dew pressure, and the system is in a two-phase region. The calculation will return the phase compositions and flow rates based on the specified temperature and pressure.

        flash case:
        - P[bubble] > P[flash] > P[dew] results in the two phase feed
        - P[bubble] < P[flash] results in the liquid phase feed
        - P[dew] > P[flash] results in the vapor phase feed
        '''
        try:
            # ! Start timing
            start_time = time.time()

            # SECTION: check inputs
            # check inputs
            if not isinstance(inputs, dict):
                raise ValueError("Inputs must be a dictionary.")

            # check keys
            if 'mole_fraction' not in inputs or 'pressure' not in inputs or 'temperature' not in inputs:
                raise ValueError(
                    "Inputs must contain 'mole_fraction', 'pressure', and 'temperature' keys.")

            # NOTE: set inputs - dict
            mole_fractions = inputs['mole_fraction']
            # NOTE: pressure [Pa]
            pressure = pycuc.convert_from_to(
                inputs['pressure'][0], inputs['pressure'][1], 'Pa')
            # NOTE: temperature [K]
            temperature = pycuc.convert_from_to(
                inputs['temperature'][0], inputs['temperature'][1], 'K')

            # check mole fractions
            if not isinstance(mole_fractions, dict):
                raise ValueError("Mole fractions must be a dictionary.")

            if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be numeric.")

            if not all(0 <= v <= 1 for v in mole_fractions.values()):
                raise ValueError("Mole fractions must be between 0 and 1.")

            # check pressure
            if not isinstance(pressure, (int, float)):
                raise ValueError("pressure must be numeric.")
            # check temperature
            if not isinstance(temperature, (int, float)):
                raise ValueError("temperature must be numeric.")

            # NOTE: components
            components = self.components

            # mole fractions based on components id
            mole_fractions = [mole_fractions[component]
                              for component in components]

            # SECTION: extract source
            Source_ = Source(self.model_source)

            # NOTE: vapor pressure source
            VaPr_comp = {}

            # looping through components
            for component in components:
                # NOTE: equation source
                # antoine equations [Pa]
                VaPr_eq = Source_.eq_extractor(component, 'VaPr')

                # NOTE: args
                VaPr_args = VaPr_eq.args
                # check args (SI)
                VaPr_args_required = Source_.check_args(component, VaPr_args)

                # build args
                _VaPr_args = Source_.build_args(component, VaPr_args_required)

                # NOTE: update P and T
                _VaPr_args['T'] = temperature

                # set
                VaPr_comp[component] = {
                    "value": VaPr_eq,
                    "args": _VaPr_args,
                    "return": VaPr_eq.returns
                }

            # SECTION: set models
            # check equilibrium model
            res_ = self.__set_models(
                equilibrium_model=equilibrium_model,
                fugacity_model=fugacity_model,
                activity_model=activity_model
            )

            # NOTE: parameters
            params = {
                "components": components,
                "mole_fraction": np.array(mole_fractions),
                "pressure": {
                    "value": pressure,
                    "unit": "Pa"
                },
                "temperature": {
                    "value": temperature,
                    "unit": "K"
                },
                "vapor_pressure": VaPr_comp,
                "equilibrium_model": equilibrium_model,
                "fugacity_model": res_['fugacity_model'],
                "activity_model": res_['activity_model'],
                "solver_method": solver_method
            }

            # SECTION: flash checker
            flash_checker_res_ = None
            if flash_checker:
                # check
                flash_checker_res_ = self._flash_checker(
                    z_i=mole_fractions,
                    Pf=pressure,
                    Tf=temperature,
                    VaPr_comp=VaPr_comp,
                )

                # check
                if flash_checker_res_ is False:
                    raise ValueError(
                        "Flash calculation failed! The system is not in a two-phase region.")

            # SECTION: flash calculation
            res = self._IFL(params, **kwargs)

            # NOTE: set message
            message = message if message is not None else "Flash Isothermal Calculation"
            # add
            res['message'] = message

            # NOTE: components
            res['components'] = components

            # NOTE: models
            res["equilibrium_model"] = equilibrium_model
            res["flash_checker"] = flash_checker
            res["flash_checker_res"] = flash_checker_res_
            res["fugacity_model"] = res_['fugacity_model']
            res["activity_model"] = res_['activity_model']
            res["solver_method"] = solver_method

            # NOTE: set time
            # ! Stop timing
            end_time = time.time()
            computation_time = end_time - start_time
            # add to res
            res['computation_time'] = {
                "value": computation_time,
                "unit": "s"
            }

            # res
            return res
        except Exception as e:
            raise Exception(
                f"Error in flash_isothermal calculation: {e}")

datasource: Dict property

Get the datasource property.

Returns

dict The datasource dictionary.

equationsource: Dict property

Get the equationsource property.

Returns

dict The equationsource dictionary.

model_source: Dict property

Get the model source property.

Returns

dict The model source dictionary.

__init__(components, model_source=None, **kwargs)

Initialize the Partition class.

Source code in pyThermoFlash/docs/vle.py

def __init__(self,
             components: List[str],
             model_source: Optional[Dict] = None,
             **kwargs):
    '''
    Initialize the Partition class.
    '''
    self.__model_source = model_source
    self.__datasource = {
    } if model_source is None else model_source['datasource']
    self.__equationsource = {
    } if model_source is None else model_source['equationsource']

    # NOTE: init class

    Equilibria.__init__(self, components, model_source)

__set_models(equilibrium_model, fugacity_model, activity_model)

Set the models for the calculation.

Parameters

equilibrium_model : str The equilibrium model to use for the calculation. fugacity_model : str, optional The fugacity model to use for the calculation. activity_model : str, optional The activity coefficient model to use for the calculation.

Returns

dict A dictionary containing the models used for the calculation. - equilibrium_model : str The equilibrium model used for the calculation. - fugacity_model : str, optional The fugacity model used for the calculation. - activity_model : str, optional The activity coefficient model used for the calculation.

Source code in pyThermoFlash/docs/vle.py

def __set_models(self,
                 equilibrium_model: str,
                 fugacity_model: Optional[str],
                 activity_model: Optional[str]):
    '''
    Set the models for the calculation.

    Parameters
    ----------
    equilibrium_model : str
        The equilibrium model to use for the calculation.
    fugacity_model : str, optional
        The fugacity model to use for the calculation.
    activity_model : str, optional
        The activity coefficient model to use for the calculation.

    Returns
    -------
    dict
        A dictionary containing the models used for the calculation.
        - equilibrium_model : str
            The equilibrium model used for the calculation.
        - fugacity_model : str, optional
            The fugacity model used for the calculation.
        - activity_model : str, optional
            The activity coefficient model used for the calculation.
    '''
    try:
        # NOTE: set based on equilibrium model
        if equilibrium_model == 'raoult':
            fugacity_model_ = None
            activity_model_ = None
        elif equilibrium_model == 'modified-raoult':
            fugacity_model_ = None
            activity_model_ = activity_model
        elif equilibrium_model == 'fugacity-ratio':
            fugacity_model_ = fugacity_model
            activity_model_ = None
        else:
            raise ValueError(
                f"Invalid equilibrium model: {equilibrium_model}.")

        # res
        return {
            "equilibrium_model": equilibrium_model,
            "fugacity_model": fugacity_model_,
            "activity_model": activity_model_
        }

    except Exception as e:
        raise Exception(
            f"Error in setting models: {e}")

bubble_pressure(inputs, equilibrium_model='raoult', fugacity_model=None, activity_model=None, message=None, **kwargs)

The Bubble-Pressure (BP) calculation determines the pressure at which the first bubble of vapor forms when a liquid mixture is heated at a constant temperature. It is used to find the pressure for a given temperature at which the liquid will begin to vaporize.

Parameters

inputs : dict Dictionary containing the input parameters for the calculation. - mole_fraction : dict Dictionary of component names and their respective mole fractions. - temperature : float Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit). equilibrium_model : str, optional The equilibrium model to use for the calculation. Default is 'raoult'. fugacity_model : str, optional The fugacity model to use for the calculation. Default is 'SRK'. activity_model : str, optional The activity coefficient model to use for the calculation. Default is 'NRTL'. message : str, optional Message to display during the calculation. Default is None. **kwargs : dict, optional Additional parameters for the model. - activity_inputs : dict, model parameters. - activity_coefficients : dict, activity coefficients for each component.

Returns

res : dict Dictionary containing the results of the calculation. - bubble_pressure: bubble pressure [Pa] - temperature: temperature [K] - feed_mole_fraction: liquid mole fraction (xi) - vapor_mole_fraction: vapor mole fraction (yi) - liquid_mole_fraction: liquid mole fraction (xi) - vapor_pressure: vapor pressure [Pa] - activity_coefficient: activity coefficient [dimensionless] - K_ratio: K-ratio [dimensionless] - equilibrium_model: equilibrium model used for the calculation - fugacity_model: fugacity model used for the calculation - activity_model: activity coefficient model used for the calculation - components: list of components used in the calculation - message: message displayed during the calculation - computation_time: computation time [s]

Notes

• Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
• Mole fractions must be between 0 and 1.
• The function will raise a ValueError if the inputs are not valid.
• The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
• The function will return the bubble pressure in Pascals.
• The activity coefficient for each component can be directly defined in the activity_inputs dictionary.

For each component, the model parameters can be provided for activity models like NRTL or UNIQUAC:

# activity inputs
activity_inputs = {
    'alpha': alpha,
    'a_ij': a_ij,
    'b_ij': b_ij,
    'c_ij': c_ij,
    'd_ij': d_ij
}

Or the activity coefficients can be provided directly:

calculated activity coefficients (for bubble-pressure calculation)

activity_coefficients = {
    'benzene': 1.0,
    'toluene': 1.1
}

Source code in pyThermoFlash/docs/vle.py

def bubble_pressure(self,
                    inputs: Dict[str, Any],
                    equilibrium_model: Literal[
                        'raoult', 'modified-raoult'
                    ] = 'raoult',
                    fugacity_model: Optional[Literal[
                        'vdW', 'PR', 'RK', 'SRK'
                    ]] = None,
                    activity_model: Optional[Literal[
                        'NRTL', 'UNIQUAC'
                    ]] = None,
                    message: Optional[str] = None,
                    **kwargs):
    '''
    The Bubble-Pressure (BP) calculation determines the pressure at which the first bubble of vapor forms when a liquid mixture is heated at a constant temperature. It is used to find the pressure for a given temperature at which the liquid will begin to vaporize.

    Parameters
    ----------
    inputs : dict
        Dictionary containing the input parameters for the calculation.
        - mole_fraction : dict
            Dictionary of component names and their respective mole fractions.
        - temperature : float
            Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
    equilibrium_model : str, optional
        The equilibrium model to use for the calculation. Default is 'raoult'.
    fugacity_model : str, optional
        The fugacity model to use for the calculation. Default is 'SRK'.
    activity_model : str, optional
        The activity coefficient model to use for the calculation. Default is 'NRTL'.
    message : str, optional
        Message to display during the calculation. Default is None.
    **kwargs : dict, optional
        Additional parameters for the model.
        - activity_inputs : dict, model parameters.
        - activity_coefficients : dict, activity coefficients for each component.

    Returns
    -------
    res : dict
        Dictionary containing the results of the calculation.
        - bubble_pressure: bubble pressure [Pa]
        - temperature: temperature [K]
        - feed_mole_fraction: liquid mole fraction (xi)
        - vapor_mole_fraction: vapor mole fraction (yi)
        - liquid_mole_fraction: liquid mole fraction (xi)
        - vapor_pressure: vapor pressure [Pa]
        - activity_coefficient: activity coefficient [dimensionless]
        - K_ratio: K-ratio [dimensionless]
        - equilibrium_model: equilibrium model used for the calculation
        - fugacity_model: fugacity model used for the calculation
        - activity_model: activity coefficient model used for the calculation
        - components: list of components used in the calculation
        - message: message displayed during the calculation
        - computation_time: computation time [s]

    Notes
    -----
    - Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
    - Mole fractions must be between 0 and 1.
    - The function will raise a ValueError if the inputs are not valid.
    - The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
    - The function will return the bubble pressure in Pascals.
    - The activity coefficient for each component can be directly defined in the activity_inputs dictionary.


    For each component, the model parameters can be provided for activity models like NRTL or UNIQUAC:

    ```python
    # activity inputs
    activity_inputs = {
        'alpha': alpha,
        'a_ij': a_ij,
        'b_ij': b_ij,
        'c_ij': c_ij,
        'd_ij': d_ij
    }
    ```

    Or the activity coefficients can be provided directly:

    # calculated activity coefficients (for bubble-pressure calculation)
    ```python
    activity_coefficients = {
        'benzene': 1.0,
        'toluene': 1.1
    }
    ```
    '''
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check inputs
        # check inputs
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")

        # check keys
        if 'mole_fraction' not in inputs or 'temperature' not in inputs:
            raise ValueError(
                "Inputs must contain 'mole_fraction' and 'temperature' keys.")

        # NOTE: set inputs - dict
        mole_fractions = inputs['mole_fraction']
        # NOTE: temperature [K]
        temperature = pycuc.convert_from_to(
            inputs['temperature'][0], inputs['temperature'][1], 'K')

        # check mole fractions
        if not isinstance(mole_fractions, dict):
            raise ValueError("Mole fractions must be a dictionary.")

        if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be numeric.")

        if not all(0 <= v <= 1 for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be between 0 and 1.")

        # check temperature
        if not isinstance(temperature, (int, float)):
            raise ValueError("Temperature must be numeric.")

        # NOTE: components
        components = self.components

        # mole fractions based on components id
        mole_fractions = [mole_fractions[component]
                          for component in components]

        # SECTION: extract source
        Source_ = Source(self.model_source)

        # NOTE: vapor pressure source
        VaPr_comp = {}

        # looping through components
        for component in components:
            # NOTE: equation source
            # antoine equations [Pa]
            VaPr_eq = Source_.eq_extractor(component, 'VaPr')

            # NOTE: args
            VaPr_args = VaPr_eq.args
            # check args (SI)
            VaPr_args_required = Source_.check_args(component, VaPr_args)

            # build args
            _VaPr_args = Source_.build_args(component, VaPr_args_required)

            # NOTE: update P and T
            _VaPr_args['T'] = temperature

            # NOTE: execute
            _VaPr_res = VaPr_eq.cal(**_VaPr_args)
            # extract
            _VaPr_value = _VaPr_res['value']
            _VaPr_unit = _VaPr_res['unit']
            # unit conversion
            # NOTE: unit conversion
            _unit_block = f"{_VaPr_unit} => Pa"
            _VaPr = pycuc.to(_VaPr_value, _unit_block)
            # set
            VaPr_comp[component] = {
                "value": _VaPr,
                "unit": "Pa"}

        # SECTION: set models
        # check equilibrium model
        res_ = self.__set_models(
            equilibrium_model=equilibrium_model,
            fugacity_model=fugacity_model,
            activity_model=activity_model
        )

        # NOTE: parameters
        params = {
            "components": components,
            "mole_fraction": np.array(mole_fractions),
            "temperature": {
                "value": temperature,
                "unit": "K"
            },
            "vapor_pressure": VaPr_comp,
            "equilibrium_model": equilibrium_model,
            "fugacity_model": res_['fugacity_model'],
            "activity_model": res_['activity_model'],
        }

        # res
        res = self._BP(params, **kwargs)

        # NOTE: set message
        message = message if message is not None else "Bubble Pressure Calculation"
        # add
        res['message'] = message

        # NOTE: components
        res['components'] = components

        # NOTE: models
        res["equilibrium_model"] = equilibrium_model
        res["fugacity_model"] = res_['fugacity_model']
        res["activity_model"] = res_['activity_model']

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # returns
        return res
    except Exception as e:
        raise ValueError(
            f"Error in bubble_pressure calculation: {e}")

bubble_temperature(inputs, equilibrium_model='raoult', fugacity_model='SRK', activity_model='NRTL', solver_method='root', message=None, **kwargs)

The bubble-point temperature (BT) calculation determines the temperature at which the first bubble of vapor forms when a liquid mixture is heated at a constant pressure. It helps identify the temperature at which the liquid will start vaporizing.

Parameters

inputs : dict Dictionary containing the input parameters for the calculation. - mole_fraction : dict Dictionary of component names and their respective mole fractions. - temperature : float Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit). equilibrium_model : str, optional The equilibrium model to use for the calculation. Default is 'raoult'. fugacity_model : str, optional The fugacity model to use for the calculation. Default is 'SRK'. activity_model : str, optional The activity coefficient model to use for the calculation. Default is 'NRTL'. solver_method : str, optional The solver method to use for the calculation. Default is 'root'. - 'root' : Root-finding algorithm. - 'least-squares' : Least-squares optimization algorithm. - 'fsolve' : SciPy's fsolve function. message : str, optional Message to display during the calculation. Default is None. **kwargs : dict, optional Additional parameters for the model. - guess_temperature : float Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

Returns

res : dict Dictionary containing the results of the calculation. - bubble_temperature: bubble temperature [K] - pressure: pressure [Pa] - feed_mole_fraction: liquid mole fraction (xi) - vapor_mole_fraction: vapor mole fraction (yi) - liquid_mole_fraction: liquid mole fraction (xi) - vapor_pressure: vapor pressure [Pa] - activity_coefficient: activity coefficient [dimensionless] - K_ratio: K-ratio [dimensionless] - equilibrium_model: equilibrium model used for the calculation - fugacity_model: fugacity model used for the calculation - activity_model: activity coefficient model used for the calculation - components: list of components used in the calculation - message: message displayed during the calculation - solver_method: solver method used for the calculation - computation_time: computation time [s]

Notes

The summary of the calculation is as follows:

• Known Information: Pressure (P) and mole fraction of the components in the liquid phase (xi). then zi = xi.
• Computed Information: Temperature (T) and mole fraction in the vapor phase (yi).

Source code in pyThermoFlash/docs/vle.py

def bubble_temperature(self,
                       inputs: Dict[str, Any],
                       equilibrium_model: Literal[
                           'raoult', 'modified-raoult'
                       ] = 'raoult',
                       fugacity_model: Literal[
                           'vdW', 'PR', 'RK', 'SRK'
                       ] = 'SRK',
                       activity_model: Literal[
                           'NRTL', 'UNIQUAC']
                       = 'NRTL',
                       solver_method: Literal[
                           'root', 'least-squares', 'fsolve'
                       ] = 'root',
                       message: Optional[str] = None,
                       **kwargs):
    '''
    The `bubble-point temperature` (BT) calculation determines the temperature at which the first bubble of vapor forms when a liquid mixture is heated at a constant pressure. It helps identify the temperature at which the liquid will start vaporizing.

    Parameters
    ----------
    inputs : dict
        Dictionary containing the input parameters for the calculation.
        - mole_fraction : dict
            Dictionary of component names and their respective mole fractions.
        - temperature : float
            Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
    equilibrium_model : str, optional
        The equilibrium model to use for the calculation. Default is 'raoult'.
    fugacity_model : str, optional
        The fugacity model to use for the calculation. Default is 'SRK'.
    activity_model : str, optional
        The activity coefficient model to use for the calculation. Default is 'NRTL'.
    solver_method : str, optional
        The solver method to use for the calculation. Default is 'root'.
        - 'root' : Root-finding algorithm.
        - 'least-squares' : Least-squares optimization algorithm.
        - 'fsolve' : SciPy's fsolve function.
    message : str, optional
        Message to display during the calculation. Default is None.
    **kwargs : dict, optional
        Additional parameters for the model.
        - `guess_temperature` : float
            Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

    Returns
    -------
    res : dict
        Dictionary containing the results of the calculation.
        - bubble_temperature: bubble temperature [K]
        - pressure: pressure [Pa]
        - feed_mole_fraction: liquid mole fraction (xi)
        - vapor_mole_fraction: vapor mole fraction (yi)
        - liquid_mole_fraction: liquid mole fraction (xi)
        - vapor_pressure: vapor pressure [Pa]
        - activity_coefficient: activity coefficient [dimensionless]
        - K_ratio: K-ratio [dimensionless]
        - equilibrium_model: equilibrium model used for the calculation
        - fugacity_model: fugacity model used for the calculation
        - activity_model: activity coefficient model used for the calculation
        - components: list of components used in the calculation
        - message: message displayed during the calculation
        - solver_method: solver method used for the calculation
        - computation_time: computation time [s]

    Notes
    -----
    The summary of the calculation is as follows:

    - Known Information: Pressure (P) and mole fraction of the components in the liquid phase (xi).
    then zi = xi.
    - Computed Information: Temperature (T) and mole fraction in the vapor phase (yi).
    '''
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check inputs
        # check inputs
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")

        # check keys
        if 'mole_fraction' not in inputs or 'pressure' not in inputs:
            raise ValueError(
                "Inputs must contain 'mole_fraction' and 'temperature' keys.")

        # NOTE: set inputs - dict
        mole_fractions = inputs['mole_fraction']
        # NOTE: pressure [Pa]
        pressure = pycuc.convert_from_to(
            inputs['pressure'][0], inputs['pressure'][1], 'Pa')

        # check mole fractions
        if not isinstance(mole_fractions, dict):
            raise ValueError("Mole fractions must be a dictionary.")

        if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be numeric.")

        if not all(0 <= v <= 1 for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be between 0 and 1.")

        # check temperature
        if not isinstance(pressure, (int, float)):
            raise ValueError("pressure must be numeric.")

        # NOTE: components
        components = self.components

        # mole fractions based on components id
        mole_fractions = [mole_fractions[component]
                          for component in components]

        # SECTION: check for activity model

        # SECTION: extract source
        Source_ = Source(self.model_source)

        # NOTE: vapor pressure source
        VaPr_comp = {}

        # looping through components
        for component in components:
            # NOTE: equation source
            VaPr_eq = None
            # antoine equations [Pa]
            VaPr_eq = Source_.eq_extractor(component, 'VaPr')

            # NOTE: args
            VaPr_args = VaPr_eq.args
            # check args (SI)
            VaPr_args_required = Source_.check_args(component, VaPr_args)

            # build args
            _VaPr_args = Source_.build_args(component, VaPr_args_required)

            # NOTE: update P and T
            _VaPr_args['T'] = None

            # set
            VaPr_comp[component] = {
                "value": VaPr_eq,
                "args": _VaPr_args,
                "return": VaPr_eq.returns
            }

        # SECTION: set models
        # check equilibrium model
        res_ = self.__set_models(
            equilibrium_model=equilibrium_model,
            fugacity_model=fugacity_model,
            activity_model=activity_model
        )

        # NOTE: parameters
        params = {
            "components": components,
            "mole_fraction": np.array(mole_fractions),
            "pressure": {
                "value": pressure,
                "unit": "Pa"
            },
            "vapor_pressure": VaPr_comp,
            "equilibrium_model": equilibrium_model,
            "fugacity_model": res_['fugacity_model'],
            "activity_model": res_['activity_model'],
            "solver_method": solver_method
        }

        # res
        res = self._BT(params, **kwargs)

        # NOTE: set message
        message = message if message is not None else "Bubble Temperature Calculation"
        # add
        res['message'] = message

        # NOTE: components
        res['components'] = components

        # NOTE: models
        res["equilibrium_model"] = equilibrium_model
        res["fugacity_model"] = res_['fugacity_model']
        res["activity_model"] = res_['activity_model']
        res["solver_method"] = solver_method

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # returns
        return res
    except Exception as e:
        raise Exception(
            f"Error in bubble_temperature calculation: {e}")

dew_pressure(inputs, equilibrium_model='raoult', fugacity_model=None, activity_model=None, message=None, **kwargs)

The dew-point pressure (DP) calculation determines the pressure at which the first drop of liquid condenses when a vapor mixture is cooled at a constant temperature. It is used to find the pressure at which vapor will begin to condense.

Parameters

inputs : dict Dictionary containing the input parameters for the calculation. - mole_fraction : dict Dictionary of component names and their respective mole fractions. - temperature : float Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit). equilibrium_model : str, optional The equilibrium model to use for the calculation. Default is 'raoult'. fugacity_model : str, optional The fugacity model to use for the calculation. Default is 'SRK'. activity_model : str, optional The activity coefficient model to use for the calculation. Default is 'NRTL'. message : str, optional Message to display during the calculation. Default is None. **kwargs : dict, optional Additional parameters for the model. - activity_inputs : dict, model parameters. - activity_coefficients : dict, activity coefficients for each component.

Returns

res : dict Dictionary containing the results of the calculation. - dew_pressure: dew pressure [Pa] - temperature: temperature [K] - feed_mole_fraction: vapor mole fraction (yi) - liquid_mole_fraction: liquid mole fraction (xi) - vapor_pressure: vapor pressure [Pa] - activity_coefficient: activity coefficient [dimensionless] - K_ratio: K-ratio [dimensionless] - equilibrium_model: equilibrium model used for the calculation - fugacity_model: fugacity model used for the calculation - activity_model: activity coefficient model used for the calculation - components: list of components used in the calculation - message: message displayed during the calculation - computation_time: computation time [s]

Notes

• Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
• Mole fractions must be between 0 and 1.
• The function will raise a ValueError if the inputs are not valid.
• The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
• The function will return the dew pressure in Pascals.

Source code in pyThermoFlash/docs/vle.py

def dew_pressure(self,
                 inputs: Dict[str, Any],
                 equilibrium_model: Literal[
                     'raoult', 'modified-raoult'
                 ] = 'raoult',
                 fugacity_model: Optional[Literal[
                     'vdW', 'PR', 'RK', 'SRK'
                 ]] = None,
                 activity_model: Optional[Literal[
                     'NRTL', 'UNIQUAC'
                 ]] = None,
                 message: Optional[str] = None,
                 **kwargs):
    '''
    The dew-point pressure (DP) calculation determines the pressure at which the first drop of liquid condenses when a vapor mixture is cooled at a constant temperature. It is used to find the pressure at which vapor will begin to condense.

    Parameters
    ----------
    inputs : dict
        Dictionary containing the input parameters for the calculation.
        - mole_fraction : dict
            Dictionary of component names and their respective mole fractions.
        - temperature : float
            Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
    equilibrium_model : str, optional
        The equilibrium model to use for the calculation. Default is 'raoult'.
    fugacity_model : str, optional
        The fugacity model to use for the calculation. Default is 'SRK'.
    activity_model : str, optional
        The activity coefficient model to use for the calculation. Default is 'NRTL'.
    message : str, optional
        Message to display during the calculation. Default is None.
    **kwargs : dict, optional
        Additional parameters for the model.
        - activity_inputs : dict, model parameters.
        - activity_coefficients : dict, activity coefficients for each component.

    Returns
    -------
    res : dict
        Dictionary containing the results of the calculation.
        - dew_pressure: dew pressure [Pa]
        - temperature: temperature [K]
        - feed_mole_fraction: vapor mole fraction (yi)
        - liquid_mole_fraction: liquid mole fraction (xi)
        - vapor_pressure: vapor pressure [Pa]
        - activity_coefficient: activity coefficient [dimensionless]
        - K_ratio: K-ratio [dimensionless]
        - equilibrium_model: equilibrium model used for the calculation
        - fugacity_model: fugacity model used for the calculation
        - activity_model: activity coefficient model used for the calculation
        - components: list of components used in the calculation
        - message: message displayed during the calculation
        - computation_time: computation time [s]

    Notes
    -----
    - Temperature must be in Kelvin, Celsius, or Fahrenheit, and will be converted to Kelvin.
    - Mole fractions must be between 0 and 1.
    - The function will raise a ValueError if the inputs are not valid.
    - The default equilibrium model is 'raoult', and the default activity model is 'NRTL'.
    - The function will return the dew pressure in Pascals.
    '''
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check inputs
        # check inputs
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")

        # check keys
        if 'mole_fraction' not in inputs or 'temperature' not in inputs:
            raise ValueError(
                "Inputs must contain 'mole_fraction' and 'temperature' keys.")

        # NOTE: set inputs - dict
        mole_fractions = inputs['mole_fraction']
        # NOTE: temperature [K]
        temperature = pycuc.convert_from_to(
            inputs['temperature'][0], inputs['temperature'][1], 'K')

        # check mole fractions
        if not isinstance(mole_fractions, dict):
            raise ValueError("Mole fractions must be a dictionary.")

        if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be numeric.")

        if not all(0 <= v <= 1 for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be between 0 and 1.")

        # check temperature
        if not isinstance(temperature, (int, float)):
            raise ValueError("Temperature must be numeric.")

        # NOTE: components
        components = self.components

        # mole fractions based on components id
        mole_fractions = [mole_fractions[component]
                          for component in components]

        # SECTION: extract source
        Source_ = Source(self.model_source)

        # NOTE: vapor pressure source
        VaPr_comp = {}

        # looping through components
        for component in components:
            # NOTE: equation source
            # antoine equations [Pa]
            VaPr_eq = Source_.eq_extractor(component, 'VaPr')

            # NOTE: args
            VaPr_args = VaPr_eq.args
            # check args (SI)
            VaPr_args_required = Source_.check_args(component, VaPr_args)

            # build args
            _VaPr_args = Source_.build_args(component, VaPr_args_required)

            # NOTE: update P and T
            _VaPr_args['T'] = temperature

            # NOTE: execute
            _VaPr_res = VaPr_eq.cal(**_VaPr_args)
            # extract
            _VaPr_value = _VaPr_res['value']
            _VaPr_unit = _VaPr_res['unit']
            # unit conversion
            # NOTE: unit conversion
            _unit_block = f"{_VaPr_unit} => Pa"
            _VaPr = pycuc.to(_VaPr_value, _unit_block)
            # set
            VaPr_comp[component] = {
                "value": _VaPr,
                "unit": "Pa",
                "equation": VaPr_eq,
                "args": _VaPr_args,
                "return": VaPr_eq.returns
            }

        # SECTION: set models
        # check equilibrium model
        res_ = self.__set_models(
            equilibrium_model=equilibrium_model,
            fugacity_model=fugacity_model,
            activity_model=activity_model
        )

        # NOTE: parameters
        params = {
            "components": components,
            "mole_fraction": np.array(mole_fractions),
            "temperature": {
                "value": temperature,
                "unit": "K"
            },
            "vapor_pressure": VaPr_comp,
            "equilibrium_model": equilibrium_model,
            "fugacity_model": res_['fugacity_model'],
            "activity_model": res_['activity_model'],
        }

        # res
        res = self._DP(params, **kwargs)

        # NOTE: set message
        message = message if message is not None else "Dew Pressure Calculation"
        # add
        res['message'] = message

        # NOTE: components
        res['components'] = components

        # NOTE: models
        res["equilibrium_model"] = equilibrium_model
        res["fugacity_model"] = res_['fugacity_model']
        res["activity_model"] = res_['activity_model']

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # returns
        return res
    except Exception as e:
        raise ValueError(
            f"Error in bubble_pressure calculation: {e}")

dew_temperature(inputs, equilibrium_model='raoult', fugacity_model='SRK', activity_model='NRTL', solver_method='least-squares', message=None, **kwargs)

The dew-point temperature (DT) calculation determines the temperature at which the first drop of liquid condenses when a vapor mixture is cooled at a constant pressure. It identifies the temperature at which vapor will start to condense.

Parameters

inputs : dict Dictionary containing the input parameters for the calculation. - mole_fraction : dict Dictionary of component names and their respective mole fractions. - temperature : float Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit). equilibrium_model : str, optional The equilibrium model to use for the calculation. Default is 'raoult'. fugacity_model : str, optional The fugacity model to use for the calculation. Default is 'SRK'. activity_model : str, optional The activity coefficient model to use for the calculation. Default is 'NRTL'. solver_method : str, optional The solver method to use for the calculation. Default is 'root'. - 'root' : Root-finding algorithm. - 'least-squares' : Least-squares optimization algorithm. - 'fsolve' : SciPy's fsolve function. message : str, optional Message to display during the calculation. Default is None. **kwargs : dict, optional Additional parameters for the model. - guess_temperature : float Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

Returns

res : dict Dictionary containing the results of the calculation. - dew_temperature: dew temperature [K] - pressure: pressure [Pa] - feed_mole_fraction: vapor mole fraction (yi) - liquid_mole_fraction: liquid mole fraction (xi) - vapor_pressure: vapor pressure [Pa] - activity_coefficient: activity coefficient [dimensionless] - K_ratio: K-ratio [dimensionless] - equilibrium_model: equilibrium model used for the calculation - fugacity_model: fugacity model used for the calculation - activity_model: activity coefficient model used for the calculation - components: list of components used in the calculation - message: message displayed during the calculation - solver_method: solver method used for the calculation - computation_time: computation time [s]

Notes

The summary of the calculation is as follows:

• Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi). then zi = yi.
• Computed Information: Temperature (T) and mole fraction in the liquid phase (xi).

Source code in pyThermoFlash/docs/vle.py

def dew_temperature(self,
                    inputs: Dict[str, Any],
                    equilibrium_model: Literal[
                        'raoult', 'modified-raoult'
                    ] = 'raoult',
                    fugacity_model: Literal[
                        'vdW', 'PR', 'RK', 'SRK'
                    ] = 'SRK',
                    activity_model: Literal[
                        'NRTL', 'UNIQUAC']
                    = 'NRTL',
                    solver_method: Literal[
                        'root', 'least-squares', 'fsolve'
                    ] = 'least-squares',
                    message: Optional[str] = None,
                    **kwargs):
    '''
    The `dew-point temperature` (DT) calculation determines the temperature at which the first drop of liquid condenses when a vapor mixture is cooled at a constant pressure. It identifies the temperature at which vapor will start to condense.

    Parameters
    ----------
    inputs : dict
        Dictionary containing the input parameters for the calculation.
        - mole_fraction : dict
            Dictionary of component names and their respective mole fractions.
        - temperature : float
            Temperature at which to calculate the bubble pressure (in Kelvin, Celsius, Fahrenheit).
    equilibrium_model : str, optional
        The equilibrium model to use for the calculation. Default is 'raoult'.
    fugacity_model : str, optional
        The fugacity model to use for the calculation. Default is 'SRK'.
    activity_model : str, optional
        The activity coefficient model to use for the calculation. Default is 'NRTL'.
    solver_method : str, optional
        The solver method to use for the calculation. Default is 'root'.
        - 'root' : Root-finding algorithm.
        - 'least-squares' : Least-squares optimization algorithm.
        - 'fsolve' : SciPy's fsolve function.
    message : str, optional
        Message to display during the calculation. Default is None.
    **kwargs : dict, optional
        Additional parameters for the model.
        - `guess_temperature` : float
            Initial guess for the bubble-point temperature (in Kelvin, Celsius, Fahrenheit), default is 295 K.

    Returns
    -------
    res : dict
        Dictionary containing the results of the calculation.
        - dew_temperature: dew temperature [K]
        - pressure: pressure [Pa]
        - feed_mole_fraction: vapor mole fraction (yi)
        - liquid_mole_fraction: liquid mole fraction (xi)
        - vapor_pressure: vapor pressure [Pa]
        - activity_coefficient: activity coefficient [dimensionless]
        - K_ratio: K-ratio [dimensionless]
        - equilibrium_model: equilibrium model used for the calculation
        - fugacity_model: fugacity model used for the calculation
        - activity_model: activity coefficient model used for the calculation
        - components: list of components used in the calculation
        - message: message displayed during the calculation
        - solver_method: solver method used for the calculation
        - computation_time: computation time [s]

    Notes
    -----
    The summary of the calculation is as follows:

    - Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi).
    then zi = yi.
    - Computed Information: Temperature (T) and mole fraction in the liquid phase (xi).
    '''
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check inputs
        # check inputs
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")

        # check keys
        if 'mole_fraction' not in inputs or 'pressure' not in inputs:
            raise ValueError(
                "Inputs must contain 'mole_fraction' and 'temperature' keys.")

        # NOTE: set inputs - dict
        mole_fractions = inputs['mole_fraction']
        # NOTE: pressure [Pa]
        # set
        pressure_value = float(inputs['pressure'][0])
        pressure_unit = str(inputs['pressure'][1])

        pressure = pycuc.convert_from_to(
            pressure_value, pressure_unit, 'Pa')

        # check mole fractions
        if not isinstance(mole_fractions, dict):
            raise ValueError("Mole fractions must be a dictionary.")

        if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be numeric.")

        if not all(0 <= v <= 1 for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be between 0 and 1.")

        # check temperature
        if not isinstance(pressure, (int, float)):
            raise ValueError("pressure must be numeric.")

        # NOTE: components
        components = self.components

        # mole fractions based on components id
        mole_fractions = [mole_fractions[component]
                          for component in components]

        # SECTION: extract source
        Source_ = Source(self.model_source)

        # NOTE: vapor pressure source
        VaPr_comp = {}

        # looping through components
        for component in components:
            # NOTE: equation source
            VaPr_eq = None
            # antoine equations [Pa]
            VaPr_eq = Source_.eq_extractor(component, 'VaPr')

            # NOTE: args
            VaPr_args = VaPr_eq.args
            # check args (SI)
            VaPr_args_required = Source_.check_args(component, VaPr_args)

            # build args
            _VaPr_args = Source_.build_args(component, VaPr_args_required)

            # NOTE: update P and T
            _VaPr_args['T'] = None

            # set
            VaPr_comp[component] = {
                "value": VaPr_eq,
                "args": _VaPr_args,
                "return": VaPr_eq.returns
            }

        # SECTION: set models
        # check equilibrium model
        res_ = self.__set_models(
            equilibrium_model=equilibrium_model,
            fugacity_model=fugacity_model,
            activity_model=activity_model
        )

        # NOTE: parameters
        params = {
            "components": components,
            "mole_fraction": np.array(mole_fractions),
            "pressure": {
                "value": pressure,
                "unit": "Pa"
            },
            "vapor_pressure": VaPr_comp,
            "equilibrium_model": equilibrium_model,
            "fugacity_model": res_['fugacity_model'],
            "activity_model": res_['activity_model'],
            "solver_method": solver_method
        }

        # res
        res = self._DT(params, **kwargs)

        # NOTE: set message
        message = message if message is not None else "Dew Temperature Calculation"
        # add
        res['message'] = message

        # NOTE: components
        res['components'] = components

        # NOTE: models
        res["equilibrium_model"] = equilibrium_model
        res["fugacity_model"] = res_['fugacity_model']
        res["activity_model"] = res_['activity_model']
        res["solver_method"] = solver_method

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # returns
        return res
    except Exception as e:
        raise Exception(
            f"Error in dew_temperature calculation: {e}")

flash_isothermal(inputs, equilibrium_model='raoult', fugacity_model='SRK', activity_model='NRTL', solver_method='least_squares', flash_checker=False, message=None, **kwargs)

The isothermal-flash (IFL) calculation This calculation determines the vapor and liquid phase compositions and amounts at a specified temperature and pressure. The system is "flashed" isothermally, meaning the temperature is kept constant while the phase behavior is calculated for a mixture.

Parameters

inputs : dict Dictionary containing the input parameters for the calculation. - mole_fraction : dict Dictionary of component names and their respective mole fractions. - temperature : float Flash temperature (in Kelvin, Celsius, Fahrenheit). - pressure : float Flash pressure (in Pascal, bar, atm). equilibrium_model : str, optional the equilibrium model to use for the calculation. Default is 'raoult'. fugacity_model : str, optional The fugacity model to use for the calculation. Default is 'SRK'. activity_model : str, optional The activity coefficient model to use for the calculation. Default is 'NRTL'. solver_method : str, optional The solver method to use for the calculation. Default is 'least_squares'. - 'minimize' : Minimize the objective function. - 'least-squares' : Least-squares optimization algorithm. flash_checker : bool, optional If True, check if the system is in a two-phase region. Default is False. message : str, optional Message to display during the calculation. Default is None. **kwargs : dict, optional Additional parameters for the model. - guess_V_F_ratio: initial guess for the vapor-to-liquid ratio (V/F), default is 0.5

Returns

res : dict Dictionary containing the results of the calculation. - vapor_to_liquid_ratio: V/F ratio [dimensionless] - liquid_to_vapor_ratio: L/F ratio [dimensionless] - feed_mole_fraction: liquid mole fraction (zi) - liquid_mole_fraction: liquid mole fraction (xi) - vapor_mole_fraction: vapor mole fraction (yi) - vapor_pressure: vapor pressure [Pa] - K_ratio: K-ratio [dimensionless] - temperature: temperature [K] - pressure: pressure [Pa] - equilibrium_model: equilibrium model used for the calculation - fugacity_model: fugacity model used for the calculation - activity_model: activity coefficient model used for the calculation - components: list of components used in the calculation - message: message displayed during the calculation - solver_method: solver method used for the calculation - flash_checker: flash checker used for the calculation - flash_checker_res: flash checker result - computation_time: computation time [s]

Notes

The summary of the calculation is as follows:

• Known Information: Temperature (T), Pressure (P), and feed mole fraction of the components (zi).
• Computed Information: Vapor and liquid phase compositions (yi and xi), and the vapor-to-liquid ratio (V/F).

Flash occurs when the bubble pressure is greater than the dew pressure, and the system is in a two-phase region. The calculation will return the phase compositions and flow rates based on the specified temperature and pressure.

flash case: - P[bubble] > P[flash] > P[dew] results in the two phase feed - P[bubble] < P[flash] results in the liquid phase feed - P[dew] > P[flash] results in the vapor phase feed

Source code in pyThermoFlash/docs/vle.py

def flash_isothermal(self,
                     inputs: Dict[str, Any],
                     equilibrium_model: Literal[
                         'raoult', 'modified-raoult'
                     ] = 'raoult',
                     fugacity_model: Literal[
                         'vdW', 'PR', 'RK', 'SRK'
                     ] = 'SRK',
                     activity_model: Literal[
                         'NRTL', 'UNIQUAC']
                     = 'NRTL',
                     solver_method: Literal[
                         'least_squares', 'minimize'
                     ] = 'least_squares',
                     flash_checker: bool = False,
                     message: Optional[str] = None,
                     **kwargs):
    '''
    The `isothermal-flash` (IFL) calculation This calculation determines the vapor and liquid phase compositions and amounts at a specified temperature and pressure. The system is "flashed" isothermally, meaning the temperature is kept constant while the phase behavior is calculated for a mixture.

    Parameters
    ----------
    inputs : dict
        Dictionary containing the input parameters for the calculation.
        - mole_fraction : dict
            Dictionary of component names and their respective mole fractions.
        - temperature : float
            Flash temperature (in Kelvin, Celsius, Fahrenheit).
        - pressure : float
            Flash pressure (in Pascal, bar, atm).
    equilibrium_model : str, optional
        the equilibrium model to use for the calculation. Default is 'raoult'.
    fugacity_model : str, optional
        The fugacity model to use for the calculation. Default is 'SRK'.
    activity_model : str, optional
        The activity coefficient model to use for the calculation. Default is 'NRTL'.
    solver_method : str, optional
        The solver method to use for the calculation. Default is 'least_squares'.
        - 'minimize' : Minimize the objective function.
        - 'least-squares' : Least-squares optimization algorithm.
    flash_checker : bool, optional
        If True, check if the system is in a two-phase region. Default is False.
    message : str, optional
        Message to display during the calculation. Default is None.
    **kwargs : dict, optional
        Additional parameters for the model.
        - `guess_V_F_ratio`: initial guess for the vapor-to-liquid ratio (V/F), default is 0.5

    Returns
    -------
    res : dict
        Dictionary containing the results of the calculation.
        - vapor_to_liquid_ratio: V/F ratio [dimensionless]
        - liquid_to_vapor_ratio: L/F ratio [dimensionless]
        - feed_mole_fraction: liquid mole fraction (zi)
        - liquid_mole_fraction: liquid mole fraction (xi)
        - vapor_mole_fraction: vapor mole fraction (yi)
        - vapor_pressure: vapor pressure [Pa]
        - K_ratio: K-ratio [dimensionless]
        - temperature: temperature [K]
        - pressure: pressure [Pa]
        - equilibrium_model: equilibrium model used for the calculation
        - fugacity_model: fugacity model used for the calculation
        - activity_model: activity coefficient model used for the calculation
        - components: list of components used in the calculation
        - message: message displayed during the calculation
        - solver_method: solver method used for the calculation
        - flash_checker: flash checker used for the calculation
        - flash_checker_res: flash checker result
        - computation_time: computation time [s]

    Notes
    -----
    The summary of the calculation is as follows:

    - Known Information: Temperature (T), Pressure (P), and feed mole fraction of the components (zi).
    - Computed Information: Vapor and liquid phase compositions (yi and xi), and the vapor-to-liquid ratio (V/F).

    Flash occurs when the bubble pressure is greater than the dew pressure, and the system is in a two-phase region. The calculation will return the phase compositions and flow rates based on the specified temperature and pressure.

    flash case:
    - P[bubble] > P[flash] > P[dew] results in the two phase feed
    - P[bubble] < P[flash] results in the liquid phase feed
    - P[dew] > P[flash] results in the vapor phase feed
    '''
    try:
        # ! Start timing
        start_time = time.time()

        # SECTION: check inputs
        # check inputs
        if not isinstance(inputs, dict):
            raise ValueError("Inputs must be a dictionary.")

        # check keys
        if 'mole_fraction' not in inputs or 'pressure' not in inputs or 'temperature' not in inputs:
            raise ValueError(
                "Inputs must contain 'mole_fraction', 'pressure', and 'temperature' keys.")

        # NOTE: set inputs - dict
        mole_fractions = inputs['mole_fraction']
        # NOTE: pressure [Pa]
        pressure = pycuc.convert_from_to(
            inputs['pressure'][0], inputs['pressure'][1], 'Pa')
        # NOTE: temperature [K]
        temperature = pycuc.convert_from_to(
            inputs['temperature'][0], inputs['temperature'][1], 'K')

        # check mole fractions
        if not isinstance(mole_fractions, dict):
            raise ValueError("Mole fractions must be a dictionary.")

        if not all(isinstance(v, (int, float)) for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be numeric.")

        if not all(0 <= v <= 1 for v in mole_fractions.values()):
            raise ValueError("Mole fractions must be between 0 and 1.")

        # check pressure
        if not isinstance(pressure, (int, float)):
            raise ValueError("pressure must be numeric.")
        # check temperature
        if not isinstance(temperature, (int, float)):
            raise ValueError("temperature must be numeric.")

        # NOTE: components
        components = self.components

        # mole fractions based on components id
        mole_fractions = [mole_fractions[component]
                          for component in components]

        # SECTION: extract source
        Source_ = Source(self.model_source)

        # NOTE: vapor pressure source
        VaPr_comp = {}

        # looping through components
        for component in components:
            # NOTE: equation source
            # antoine equations [Pa]
            VaPr_eq = Source_.eq_extractor(component, 'VaPr')

            # NOTE: args
            VaPr_args = VaPr_eq.args
            # check args (SI)
            VaPr_args_required = Source_.check_args(component, VaPr_args)

            # build args
            _VaPr_args = Source_.build_args(component, VaPr_args_required)

            # NOTE: update P and T
            _VaPr_args['T'] = temperature

            # set
            VaPr_comp[component] = {
                "value": VaPr_eq,
                "args": _VaPr_args,
                "return": VaPr_eq.returns
            }

        # SECTION: set models
        # check equilibrium model
        res_ = self.__set_models(
            equilibrium_model=equilibrium_model,
            fugacity_model=fugacity_model,
            activity_model=activity_model
        )

        # NOTE: parameters
        params = {
            "components": components,
            "mole_fraction": np.array(mole_fractions),
            "pressure": {
                "value": pressure,
                "unit": "Pa"
            },
            "temperature": {
                "value": temperature,
                "unit": "K"
            },
            "vapor_pressure": VaPr_comp,
            "equilibrium_model": equilibrium_model,
            "fugacity_model": res_['fugacity_model'],
            "activity_model": res_['activity_model'],
            "solver_method": solver_method
        }

        # SECTION: flash checker
        flash_checker_res_ = None
        if flash_checker:
            # check
            flash_checker_res_ = self._flash_checker(
                z_i=mole_fractions,
                Pf=pressure,
                Tf=temperature,
                VaPr_comp=VaPr_comp,
            )

            # check
            if flash_checker_res_ is False:
                raise ValueError(
                    "Flash calculation failed! The system is not in a two-phase region.")

        # SECTION: flash calculation
        res = self._IFL(params, **kwargs)

        # NOTE: set message
        message = message if message is not None else "Flash Isothermal Calculation"
        # add
        res['message'] = message

        # NOTE: components
        res['components'] = components

        # NOTE: models
        res["equilibrium_model"] = equilibrium_model
        res["flash_checker"] = flash_checker
        res["flash_checker_res"] = flash_checker_res_
        res["fugacity_model"] = res_['fugacity_model']
        res["activity_model"] = res_['activity_model']
        res["solver_method"] = solver_method

        # NOTE: set time
        # ! Stop timing
        end_time = time.time()
        computation_time = end_time - start_time
        # add to res
        res['computation_time'] = {
            "value": computation_time,
            "unit": "s"
        }

        # res
        return res
    except Exception as e:
        raise Exception(
            f"Error in flash_isothermal calculation: {e}")

Source 🛠️

Source

Source class for the pyThermoFlash package.

Source code in pyThermoFlash/docs/source.py

class Source:
    '''
    Source class for the pyThermoFlash package.
    '''
    # NOTE: variables

    def __init__(self, model_source: Optional[Dict] = None, **kwargs):
        '''Initialize the Source class.'''
        # set
        self.model_source = model_source

        # NOTE: source
        # set
        self._datasource, self._equationsource = self.set_source(
            model_source=model_source)

    def __repr__(self):
        des = "Source class for the pyThermoFlash package."
        return des

    @property
    def datasource(self) -> Dict:
        '''
        Get the datasource property.

        Returns
        -------
        dict
            The datasource dictionary.
        '''
        # NOTE: check if model source is valid
        if self._datasource is None:
            return {}
        return self._datasource

    @property
    def equationsource(self) -> Dict:
        '''
        Get the equationsource property.

        Returns
        -------
        dict
            The equationsource dictionary.
        '''
        # NOTE: check if model source is valid
        if self._equationsource is None:
            return {}
        return self._equationsource

    def set_source(self, model_source: Dict):
        '''
        Set the model source.

        Parameters
        ----------
        model_source : dict
            The model source dictionary.
        '''
        # NOTE: source
        # datasource
        _datasource = {
        } if model_source is None else model_source[DATASOURCE]

        # equationsource
        _equationsource = {
        } if model_source is None else model_source[EQUATIONSOURCE]

        # res
        return _datasource, _equationsource

    def eq_extractor(self, component_name: str, prop_name: str):
        '''
        Extracts the equilibrium property from the model source.

        Parameters
        ----------
        component_name : str
            The name of the component.
        prop_name : str
            The name of the property to extract.

        Returns
        -------
        dict
            The extracted property.
        '''
        if self.equationsource is None:
            raise ValueError("Equation source is not defined.")

        # NOTE: check component
        if component_name not in self.equationsource.keys():
            raise ValueError(
                f"Component '{component_name}' not found in model source.")

        # NOTE: check property
        if prop_name not in self.equationsource[component_name].keys():
            raise ValueError(
                f"Property '{prop_name}' not found in model source registered for {component_name}.")

        return self.equationsource[component_name][prop_name]

    def data_extractor(self, component_name: str, prop_name: str):
        '''
        Extracts the data property from the model source.

        Parameters
        ----------
        component_name : str
            The name of the component.
        prop_name : str
            The name of the property to extract.

        Returns
        -------
        dict
            The extracted property.
        '''
        if self.datasource is None:
            return None

        # NOTE: check component
        if component_name not in self.datasource.keys():
            raise ValueError(
                f"Component '{component_name}' not found in model datasource.")

        # NOTE: check property
        if prop_name not in self.datasource[component_name].keys():
            raise ValueError(
                f"Property '{prop_name}' not found in model datasource registered for {component_name}.")

        return self.datasource[component_name][prop_name]

    def check_args(self, component_name: str, args):
        '''
        Checks equation args

        Parameters
        ----------
        component_name : str
            The name of the component.
        args : tuple
            equation args
        '''
        try:
            # required args
            required_args = []

            # datasource list
            datasource_component_list = list(
                self.datasource[component_name].keys())

            # NOTE: default args
            datasource_component_list.append("P")
            datasource_component_list.append("T")

            # check args within datasource
            for arg_key, arg_value in args.items():
                # symbol
                if arg_value['symbol'] in datasource_component_list:
                    # update
                    required_args.append(arg_value)
                else:
                    raise Exception('Args not in datasource!')

            # res
            return required_args

        except Exception as e:
            raise Exception('Finding args failed!, ', e)

    def build_args(self, component_name: str,
                   args,
                   ignore_symbols: List[str] = ["T", "P"]):
        '''
        Builds args

        Parameters
        ----------
        component_name : str
            The name of the component.
        args : tuple
            equation args
        ignore_symbols : list
            list of symbols to ignore, default is ["T", "P"]
        '''
        try:
            # res
            res = {}
            for arg in args:
                # symbol
                symbol = arg['symbol']

                # NOTE: check if symbol is in ignore symbols
                if ignore_symbols is not None:
                    # check in ignore symbols
                    if symbol not in ignore_symbols:
                        # check in component database
                        for key, value in self.datasource.items():
                            if symbol == key:
                                res[symbol] = value
                else:
                    # check in component database
                    for key, value in self.datasource.items():
                        if symbol == key:
                            res[symbol] = value
            return res
        except Exception as e:
            raise Exception('Building args failed!, ', e)

datasource: Dict property

Get the datasource property.

Returns

dict The datasource dictionary.

equationsource: Dict property

Get the equationsource property.

Returns

dict The equationsource dictionary.

__init__(model_source=None, **kwargs)

Initialize the Source class.

Source code in pyThermoFlash/docs/source.py

def __init__(self, model_source: Optional[Dict] = None, **kwargs):
    '''Initialize the Source class.'''
    # set
    self.model_source = model_source

    # NOTE: source
    # set
    self._datasource, self._equationsource = self.set_source(
        model_source=model_source)

build_args(component_name, args, ignore_symbols=['T', 'P'])

Builds args

Parameters

component_name : str The name of the component. args : tuple equation args ignore_symbols : list list of symbols to ignore, default is ["T", "P"]

Source code in pyThermoFlash/docs/source.py

def build_args(self, component_name: str,
               args,
               ignore_symbols: List[str] = ["T", "P"]):
    '''
    Builds args

    Parameters
    ----------
    component_name : str
        The name of the component.
    args : tuple
        equation args
    ignore_symbols : list
        list of symbols to ignore, default is ["T", "P"]
    '''
    try:
        # res
        res = {}
        for arg in args:
            # symbol
            symbol = arg['symbol']

            # NOTE: check if symbol is in ignore symbols
            if ignore_symbols is not None:
                # check in ignore symbols
                if symbol not in ignore_symbols:
                    # check in component database
                    for key, value in self.datasource.items():
                        if symbol == key:
                            res[symbol] = value
            else:
                # check in component database
                for key, value in self.datasource.items():
                    if symbol == key:
                        res[symbol] = value
        return res
    except Exception as e:
        raise Exception('Building args failed!, ', e)

check_args(component_name, args)

Checks equation args

Parameters

component_name : str The name of the component. args : tuple equation args

Source code in pyThermoFlash/docs/source.py

def check_args(self, component_name: str, args):
    '''
    Checks equation args

    Parameters
    ----------
    component_name : str
        The name of the component.
    args : tuple
        equation args
    '''
    try:
        # required args
        required_args = []

        # datasource list
        datasource_component_list = list(
            self.datasource[component_name].keys())

        # NOTE: default args
        datasource_component_list.append("P")
        datasource_component_list.append("T")

        # check args within datasource
        for arg_key, arg_value in args.items():
            # symbol
            if arg_value['symbol'] in datasource_component_list:
                # update
                required_args.append(arg_value)
            else:
                raise Exception('Args not in datasource!')

        # res
        return required_args

    except Exception as e:
        raise Exception('Finding args failed!, ', e)

data_extractor(component_name, prop_name)

Extracts the data property from the model source.

Parameters

component_name : str The name of the component. prop_name : str The name of the property to extract.

Returns

dict The extracted property.

Source code in pyThermoFlash/docs/source.py

def data_extractor(self, component_name: str, prop_name: str):
    '''
    Extracts the data property from the model source.

    Parameters
    ----------
    component_name : str
        The name of the component.
    prop_name : str
        The name of the property to extract.

    Returns
    -------
    dict
        The extracted property.
    '''
    if self.datasource is None:
        return None

    # NOTE: check component
    if component_name not in self.datasource.keys():
        raise ValueError(
            f"Component '{component_name}' not found in model datasource.")

    # NOTE: check property
    if prop_name not in self.datasource[component_name].keys():
        raise ValueError(
            f"Property '{prop_name}' not found in model datasource registered for {component_name}.")

    return self.datasource[component_name][prop_name]

eq_extractor(component_name, prop_name)

Extracts the equilibrium property from the model source.

Parameters

component_name : str The name of the component. prop_name : str The name of the property to extract.

Returns

dict The extracted property.

Source code in pyThermoFlash/docs/source.py

def eq_extractor(self, component_name: str, prop_name: str):
    '''
    Extracts the equilibrium property from the model source.

    Parameters
    ----------
    component_name : str
        The name of the component.
    prop_name : str
        The name of the property to extract.

    Returns
    -------
    dict
        The extracted property.
    '''
    if self.equationsource is None:
        raise ValueError("Equation source is not defined.")

    # NOTE: check component
    if component_name not in self.equationsource.keys():
        raise ValueError(
            f"Component '{component_name}' not found in model source.")

    # NOTE: check property
    if prop_name not in self.equationsource[component_name].keys():
        raise ValueError(
            f"Property '{prop_name}' not found in model source registered for {component_name}.")

    return self.equationsource[component_name][prop_name]

set_source(model_source)

Set the model source.

Parameters

model_source : dict The model source dictionary.

Source code in pyThermoFlash/docs/source.py

def set_source(self, model_source: Dict):
    '''
    Set the model source.

    Parameters
    ----------
    model_source : dict
        The model source dictionary.
    '''
    # NOTE: source
    # datasource
    _datasource = {
    } if model_source is None else model_source[DATASOURCE]

    # equationsource
    _equationsource = {
    } if model_source is None else model_source[EQUATIONSOURCE]

    # res
    return _datasource, _equationsource

Equilibria ⚖️

Equilibria

Phase Equilibria Calculations

Source code in pyThermoFlash/docs/equilibria.py

class Equilibria:
    '''
    Phase Equilibria Calculations
    '''

    def __init__(self,
                 components: List[str],
                 datasource: Optional[Dict] = None,
                 equationsource: Optional[Dict] = None,
                 **kwargs):
        '''Initialize the Equilibria class.'''
        # set
        components_ = [component.strip() for component in components]
        self.__components = components_
        self.__comp_num = len(components_)

        # model source
        self._datasource = datasource
        self._equationsource = equationsource

        # NOTE: init activity model
        self.Activity_ = Activity(
            datasource=self._datasource,
            equationsource=self._equationsource
        )

    def __repr__(self):
        des = "Phase Equilibria Calculations"
        return des

    @property
    def components(self):
        '''Get the components of the system.'''
        return self.__components

    @property
    def component_num(self):
        '''Get the number of components in the system.'''
        return self.__comp_num

    def __mole_fraction_comp(self,
                             z_i: List[float] | np.ndarray
                             ) -> Dict[str, float]:
        '''
        Convert mole fraction list to dictionary.

        Parameters
        ----------
        z_i : list | np.ndarray
            List of mole fractions of each component in the mixture.

        Returns
        -------
        dict
            Dictionary of mole fractions.
        '''
        # NOTE: check type
        if isinstance(z_i, np.ndarray):
            # convert to list
            z_i = [float(z) for z in z_i]
        # convert to dict
        return {str(self.components[i]): z_i[i]
                for i in range(self.component_num)}

    def __activity_coefficient(self,
                               activity_model: Literal['NRTL', 'UNIQUAC'],
                               activity: Activity,
                               x_i_comp: Dict[str, float],
                               T: float,
                               **kwargs
                               ) -> np.ndarray:
        '''
        Calculate activity coefficient using NRTL or UNIQUAC model.

        Parameters
        ----------
        activity_model : str
            Activity model to be used ('NRTL' or 'UNIQUAC').
        activity : Activity
            Activity object for calculating activity coefficients.
        x_i_comp : dict
            Dictionary of mole fractions of each component in the mixture.
        T : float
            Temperature [K].
        kwargs : dict
            Additional parameters for the calculation.

        Returns
        -------
        AcCo_i : array-like
            Activity coefficient of each component in the mixture.
        '''
        try:
            # NOTE: check model
            if activity_model == 'NRTL':
                # calculate activity
                res_ = activity.NRTL(
                    self.components,
                    x_i_comp,
                    T,
                    **kwargs)
                # extract
                AcCo_i = res_['value']
            elif activity_model == 'UNIQUAC':
                # calculate activity
                res_ = activity.UNIQUAC(
                    self.components,
                    x_i_comp,
                    T,
                    **kwargs)
                # extract
                AcCo_i = res_['value']
            else:
                # equals unity for ideal solution
                AcCo_i = np.ones(self.component_num)

            # res
            return AcCo_i
        except Exception as e:
            raise Exception(f'activity coefficient calculation failed! {e}')

    def __check_activity_coefficients(self,
                                      activity_model: Literal['NRTL', 'UNIQUAC'],
                                      activity: Activity,
                                      z_i_comp: Dict[str, float],
                                      T_value: float,
                                      **kwargs
                                      ) -> np.ndarray:
        '''
        Check if activity coefficients are provided by the user.

        Returns
        -------
        AcCo_i : np.ndarray
            Activity coefficients as a numpy array.
        '''
        try:
            # SECTION: activity coefficient
            # ! check kwargs
            activity_coefficients_ext = kwargs.get(
                'activity_coefficients', None)

            # check calculated activity coefficients
            if activity_coefficients_ext:
                # set
                AcCo_i_ = activity_coefficients_ext

                # ! set activity coefficients
                AcCo_i = self.set_calculated_activity_coefficient(
                    AcCo_i_)
            else:
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    z_i_comp,
                    T_value,
                    **kwargs
                )

            # res
            return AcCo_i
        except Exception as e:
            raise Exception(f'check activity coefficients failed! {e}')

    def set_calculated_activity_coefficient(self,
                                            AcCo_i_cal:
                                            Dict[str, float] | np.ndarray | list,
                                            ) -> np.ndarray:
        """
        Set calculated activity coefficients, provided by the user.

        Parameters
        ----------
        AcCo_i_cal : dict | np.ndarray | list
            Activity coefficients for each component in the mixture.
            If a dict is provided, it should contain component names as keys
            and their corresponding activity coefficients as values.
            If a numpy array or list is provided, it should contain activity
            coefficients in the same order as the components.

        Returns
        -------
        AcCo_i : np.ndarray
            Activity coefficients as a numpy array.
        """
        try:
            # NOTE: check type
            if isinstance(AcCo_i_cal, dict):
                # convert to array based on components
                AcCo_i = np.zeros(self.component_num)
                for i, component in enumerate(self.components):
                    if component in AcCo_i_cal:
                        AcCo_i[i] = AcCo_i_cal[component]
                    else:
                        raise Exception(
                            f"activity_coefficients for {component} not found!")
            elif isinstance(AcCo_i_cal, np.ndarray):
                AcCo_i = AcCo_i_cal
            elif isinstance(AcCo_i_cal, list):
                AcCo_i = np.array(AcCo_i_cal)
            else:
                raise Exception(
                    'activity_coefficients must be a dict, list or numpy array!')
            return AcCo_i
        except Exception as e:
            raise Exception(f'set activity coefficient failed! {e}')

    def _BP(self, params, **kwargs):
        '''
        The bubble-point pressure (BP) calculation determines the pressure at which the first bubble of vapor forms when a liquid mixture is heated at a constant temperature. It is used to find the pressure for a given temperature at which the liquid will begin to vaporize. This calculation is based on `Raoult's law`, which states that the vapor pressure of each component in the mixture is proportional to its mole fraction in the liquid phase.

        Parameters
        ----------
        params : dict
            Dictionary containing the following:
            - zi : liquid mole fraction
            - T: temperature [K]
        kwargs : dict
            additional parameters for the calculation
                - `activity_inputs`: information for the calculation

        Returns
        -------
        dict
            Dictionary containing the following:
            - bubble_pressure: bubble pressure [Pa]
            - temperature: temperature [K]
            - feed_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: `Temperature` (T) and `mole fraction of the components in the liquid phase` (xi).
        then zi = xi.
        - Computed Information: Pressure (P) and mole fraction in the vapor phase (yi).

        - The solution is obtained by the following steps:
            1. Calculate the vapor pressure of each component at the given temperature (T).
            2. Calculate the bubble pressure (P) using `Raoult's law` or other methods.
            3. Calculate the mole fraction in the vapor phase (yi) using `Raoult's law` or other methods.

        - inputs hold the information for the calculation:
            `alpha_ij`: binary interaction parameter
            `dg_ij`: dg_ij parameter
            `dU_ij`: dU_ij parameter
            `a_ij`: a_ij parameter
            `b_ij`: b_ij parameter
            `c_ij`: c_ij parameter
            `d_ij`: d_ij parameter
        '''
        try:
            # SECTION: data
            # NOTE: params
            # mole fraction [array]
            z_i = params['mole_fraction']
            # temperature [K]
            T = params['temperature']
            # equilibrium model
            eq_model = params['equilibrium_model']
            # fugacity model
            fugacity_model = params['fugacity_model']
            # activity model
            activity_model = params['activity_model']

            # NOTE: vapor pressure [Pa]
            VaPr_comp = params['vapor_pressure']

            # NOTE: mole fraction dict
            # convert to dict
            z_i_comp = self.__mole_fraction_comp(z_i)

            # temperature [K]
            T_value = T['value']

            # SECTION: activity coefficient
            # ! check kwargs
            activity_coefficients_ext = kwargs.get(
                'activity_coefficients', None)
            # check calculated activity coefficients
            if activity_coefficients_ext:
                # set
                AcCo_i_ = activity_coefficients_ext

                # ! set activity coefficients
                AcCo_i = self.set_calculated_activity_coefficient(
                    AcCo_i_)
            else:
                # NOTE: init model
                # init NRTL model
                # activity = Activity(
                #     datasource=self._datasource,
                #     equationsource=self._equationsource
                # )

                activity = self.Activity_

                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    z_i_comp,
                    T_value,
                    **kwargs
                )

            # SECTION: vapor pressure calculation
            # vapor pressure [Pa]
            VaPr_i = np.zeros(self.component_num)

            # looping over components
            for i, component in enumerate(self.components):
                # vapor pressure [Pa]
                VaPr_i[i] = VaPr_comp[component]['value']

            # bubble pressure [Pa]
            BuPr = np.sum(AcCo_i*z_i*VaPr_i)

            # vapor mole fraction
            y_i = np.zeros(self.component_num)

            # looping over components
            for i in range(self.component_num):
                y_i[i] = z_i[i]*VaPr_i[i]*AcCo_i[i]/BuPr

            # NOTE: Ki ratio
            K_i = np.multiply(y_i, 1/z_i)

            # NOTE: results
            res = {
                "bubble_pressure": {
                    "value": float(BuPr),
                    "unit": "Pa"
                },
                "temperature": {
                    "value": T_value,
                    "unit": "K"
                },
                "feed_mole_fraction": z_i,
                "vapor_mole_fraction": y_i,
                "liquid_mole_fraction": z_i,
                "mole_fraction_sum": {
                    "zi": float(np.sum(z_i)),
                    "xi": float(np.sum(z_i)),
                    "yi": float(np.sum(y_i))
                },
                "vapor_pressure": {
                    "value": VaPr_i,
                    "unit": "Pa"
                },
                "activity_coefficient": {
                    "value": AcCo_i,
                    "unit": "dimensionless"
                },
                "K_ratio": {
                    "value": K_i,
                    "unit": "dimensionless"
                }
            }

            # res
            return res
        except Exception as e:
            raise Exception(f'bubble pressure calculation failed! {e}')

    def _DP(self, params, **kwargs):
        '''
        The dew-point pressure (DP) calculation determines the pressure at which the first drop of liquid condenses when a vapor mixture is cooled at a constant temperature. It is used to find the pressure at which vapor will begin to condense.

        Parameters
        ----------
        params : dict
            Dictionary containing the following:
            - zi : vapor mole fraction (yi)
            - T: temperature [K]
        kwargs : dict
            additional parameters for the calculation
            - `activity_inputs`: information for the calculation
            - `max_iter`: maximum number of iterations for the calculation, default is 500
            - `tolerance`: tolerance for the calculation, default is 1e-6

        Returns
        -------
        res : dict
            Dictionary containing the following:
            - dew_pressure: dew pressure [Pa]
            - temperature: temperature [K]
            - feed_mole_fraction: vapor mole fraction (yi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Temperature (T) and mole fraction of the components in the vapor phase (yi).
        then zi = yi.
        - Computed Information: Pressure (P) and mole fraction in the liquid phase (xi).


        - The solution is obtained by the following steps:
            1. Calculate the vapor pressure of each component at the given temperature (T).
            2. Calculate the dew pressure (P) using Raoult's law.
            3. Calculate the mole fraction in the liquid phase (xi) using `Raoult's law`.
        '''
        try:
            # SECTION: data
            # NOTE: params
            # mole fraction [array]
            y_i = params['mole_fraction']
            # temperature [K]
            T = params['temperature']
            # equilibrium model
            eq_model = params['equilibrium_model']
            # fugacity model
            fugacity_model = params['fugacity_model']
            # activity model
            activity_model = params['activity_model']

            # SECTION: kwargs
            # set max_iter
            max_iter = kwargs.get('max_iter', 500)
            # set tolerance
            tolerance = kwargs.get('tolerance', 1e-6)

            # SECTION: activity coefficient
            # NOTE: mole fraction dict
            # convert to dict
            y_i_comp = self.__mole_fraction_comp(y_i)

            # temperature [K]
            T_value = T['value']

            # NOTE: vapor pressure [Pa]
            VaPr_comp = params['vapor_pressure']

            # vapor pressure [Pa]
            VaPr_i = np.zeros(self.component_num)

            # NOTE: activity coefficient [dimensionless]
            AcCo_i = np.ones(self.component_num)

            # iteration counter
            m = 0

            # SECTION: check vle model
            if eq_model == 'raoult':
                # ! Raoult's law
                # looping over components
                for i, component in enumerate(self.components):
                    # vapor pressure [Pa]
                    VaPr_i[i] = VaPr_comp[component]['value']

                # NOTE: dew pressure [Pa]
                DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

                # NOTE: liquid mole fraction
                x_i = np.zeros(self.component_num)

                # looping over components
                for i in range(self.component_num):
                    # mole fraction
                    x_i[i] = y_i[i]*DePr/VaPr_i[i]*AcCo_i[i]

            elif eq_model == 'modified-raoult':
                # ! Modified Raoult's law
                # vapor pressure [Pa]
                # looping over components
                for i, component in enumerate(self.components):
                    # vapor pressure [Pa]
                    eq_ = VaPr_comp[component]['equation']
                    args_ = VaPr_comp[component]['args']
                    # update args
                    args_['T'] = T_value
                    # cal
                    res_ = eq_.cal(**args_)
                    # extract
                    res_value_ = res_['value']
                    res_unit_ = res_['unit']
                    # convert to Pa
                    unit_block_ = f"{res_unit_} => Pa"
                    VaPr_ = pycuc.to(res_value_, unit_block_)
                    # set
                    VaPr_i[i] = VaPr_

                # NOTE: initial dew pressure [Pa]
                DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

                # NOTE: iterate to find dew pressure while checking both DePr and x_i

                # set x loop
                x_i_comp = y_i_comp.copy()
                x_i = y_i.copy()

                # init activity model
                # activity = Activity(
                #     datasource=self._datasource,
                #     equationsource=self._equationsource
                # )

                activity = self.Activity_

                # SECTION: iteration
                for m in range(max_iter):

                    # NOTE: check activity coefficients is provided by the user
                    AcCo_i = self.__check_activity_coefficients(
                        activity_model,
                        activity,
                        x_i_comp,
                        T_value,
                        **kwargs
                    )

                    # NOTE: calculate
                    # AcCo_i = self.__activity_coefficient(
                    #     activity_model,
                    #     activity,
                    #     x_i_comp,
                    #     T_value,
                    #     **kwargs
                    # )

                    # NOTE: dew pressure [Pa]
                    DePr_new = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

                    # NOTE: liquid mole fraction
                    x_i_new = (y_i*DePr_new)/(VaPr_i*AcCo_i)
                    x_i_new /= np.sum(x_i_new)

                    # NOTE: check convergence
                    if (np.all(np.abs(x_i - x_i_new) < tolerance) and
                            np.abs(DePr - DePr_new) < tolerance):
                        break

                    # update loop
                    x_i = x_i_new
                    x_i_comp = self.__mole_fraction_comp(x_i)
                    DePr = DePr_new

            # NOTE: k-ratio
            K_i = np.multiply(y_i, 1/x_i)

            # res
            res = {
                "dew_pressure": {
                    "value": float(DePr),
                    "unit": "Pa"
                },
                "temperature": {
                    "value": T,
                    "unit": "K"
                },
                "feed_mole_fraction": y_i,
                "vapor_mole_fraction": y_i,
                "liquid_mole_fraction": x_i,
                "mole_fraction_sum": {
                    "zi": float(np.sum(y_i)),
                    "xi": float(np.sum(x_i)),
                    "yi": float(np.sum(y_i))
                },
                "vapor_pressure": {
                    "value": VaPr_i,
                    "unit": "Pa"
                },
                "activity_coefficient": {
                    "value": AcCo_i,
                    "unit": "dimensionless"
                },
                "K_ratio": {
                    "value": K_i,
                    "unit": "dimensionless"
                },
                "max_iter": max_iter,
                "iteration": m,
                "tolerance": tolerance
            }

            # res
            return res
        except Exception as e:
            raise Exception(f'dew pressure calculation failed! {e}')

    def _BT(self, params, **kwargs):
        '''
        The `bubble-point temperature` (BT) calculation determines the temperature at which the first bubble of vapor forms when a liquid mixture is heated at a constant pressure. It helps identify the temperature at which the liquid will start vaporizing.

        Parameters
        ----------
        params : dict
            Dictionary containing the following:
            - zi : liquid mole fraction (xi)
            - P: pressure [Pa]
        kwargs : dict
            additional parameters for the calculation
            - `guess_temperature`: initial guess temperature [K], default is 295 K
            - `activity_inputs`: activity model inputs consists of:
                - activity_model: activity model name (NRTL, UNIQUAC)
                - alpha_ij: binary interaction parameter
                - dg_ij: dg_ij parameter
                - dU_ij: dU_ij parameter
                - a_ij: a_ij parameter
                - b_ij: b_ij parameter
                - c_ij: c_ij parameter
                - d_ij: d_ij parameter

        Returns
        -------
        res : dict
            Dictionary containing the following:
            - bubble_temperature: bubble temperature [K]
            - pressure: pressure [Pa]
            - feed_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: `Pressure` (P) and `mole fraction of the components in the liquid phase` (xi).
        then zi = xi.
        - Computed Information: `Temperature` (T) and `mole fraction in the vapor phase` (yi).

        The solution is obtained by the following steps:
            1. Choose an initial guess for the bubble-point temperature (Tg).
            2. Calculate the vapor pressure of each component at the guessed temperature (Tg).
            3. Calculate the bubble pressure (Pb) using `Raoult's law` or other methods.
            4. Check if the calculated bubble pressure (Pb) equals the given pressure (P).
            5. Calculate the vapor mole fraction (yi) using `Raoult's law` or other methods.
        '''
        try:
            # SECTION: data
            # NOTE: params
            # mole fraction [array]
            z_i = params['mole_fraction']
            # pressure [Pa]
            P = params['pressure']
            # equilibrium model
            eq_model = params['equilibrium_model']
            # fugacity model
            fugacity_model = params['fugacity_model']
            # activity model
            activity_model = params['activity_model']
            # solver method
            solver_method = params['solver_method']

            # SECTION: vapor pressure calculation
            # NOTE: vapor pressure equation [Pa]
            VaPr_comp = params['vapor_pressure']

            # SECTION: kwargs
            # temperature guess [K]
            T_g0 = kwargs.get('guess_temperature', 295)

            # SECTION: activity coefficient
            # NOTE: mole fraction dict
            # convert to dict
            z_i_comp = self.__mole_fraction_comp(z_i)

            # NOTE: init model
            # init NRTL model
            # activity = Activity(
            #     datasource=self._datasource,
            #     equationsource=self._equationsource
            # )

            activity = self.Activity_

            # activity inputs
            activity_inputs = kwargs.get('activity_inputs', {})

            # SECTION: optimization
            # pressure [Pa]
            P_value = P['value']

            # params
            _params = {
                'mole_fraction': z_i,
                'mole_fraction_comp': z_i_comp,
                'pressure': P_value,
                'vapor_pressure': VaPr_comp,
                'activity_model': activity_model,
                'activity': activity,
                'activity_inputs': activity_inputs
            }

            # NOTE: bubble temperature [K]
            T = None

            # check solver method
            if solver_method == 'fsolve':
                # ! fsolve
                _res0 = optimize.fsolve(
                    self.fBT, T_g0, args=(_params,), full_output=True)
                # extract
                T, infodict, ier, msg = _res0

                # check if root found
                if ier != 1:
                    raise Exception(f'root not found!, {msg}')

                # check
                if len(T) == 1:
                    T = float(T[0])
            elif solver_method == 'root':
                # ! root
                _res0 = optimize.root(
                    self.fBT, T_g0, args=(_params,))

                # check
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]

            elif solver_method == 'least-squares':
                # ! least-squares
                _res0 = optimize.least_squares(
                    self.fBT, T_g0, args=(_params,))

                # check
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]

            else:
                raise Exception('solver method not found!')

            # SECTION: vapor pressure [Pa]
            # at T (Tg)
            VaPr = np.zeros(self.component_num)

            # looping over components
            for i, component in enumerate(self.components):
                # vapor pressure [Pa]
                eq_ = VaPr_comp[component]['value']
                args_ = VaPr_comp[component]['args']
                # update args
                args_['T'] = T
                # cal
                res_ = eq_.cal(**args_)
                # extract
                res_value_ = res_['value']
                res_unit_ = res_['unit']
                # convert to Pa
                unit_block_ = f"{res_unit_} => Pa"
                VaPr_ = pycuc.to(res_value_, unit_block_)
                # set
                VaPr[i] = VaPr_

            # SECTION: calculate activity coefficient
            # NOTE: calculate
            AcCo_i = self.__activity_coefficient(
                activity_model,
                activity,
                z_i_comp,
                T,
                **kwargs
            )

            # vapor mole fraction
            yi = np.zeros(self.component_num)
            for i in range(self.component_num):
                yi[i] = z_i[i]*AcCo_i[i]*VaPr[i]/P_value

            # NOTE: k-ratio
            K_i = np.multiply(yi, 1/z_i)

            # NOTE: results
            res = {
                "bubble_temperature": {
                    "value": float(T),
                    "unit": "K"
                },
                "pressure": {
                    "value": P_value,
                    "unit": "Pa"
                },
                "feed_mole_fraction": z_i,
                "liquid_mole_fraction": z_i,
                "vapor_mole_fraction": yi,
                "mole_fraction_sum": {
                    "zi": float(np.sum(z_i)),
                    "xi": float(np.sum(z_i)),
                    "yi": float(np.sum(yi))
                },
                "vapor_pressure": {
                    "value": VaPr,
                    "unit": "Pa"
                },
                "K_ratio": {
                    "value": K_i,
                    "unit": "dimensionless"
                },
                "activity_coefficient": {
                    "value": AcCo_i,
                    "unit": "dimensionless"
                }
            }

            # res
            return res
        except Exception as e:
            raise Exception(f"bubble temperature calculation failed! {e}")

    def fBT(self, x, params) -> float:
        '''
        bubble temperature function

        Parameters
        ----------
        x : array-like
            Guess temperature [K]
        params : dict
            Dictionary containing the following:
            - mole_fraction : liquid mole fraction (xi)
            - pressure : pressure [Pa]
            - vapor_pressure : vapor pressure equation [Pa]
        '''
        # NOTE: temperature (loop)
        T = x[0]

        # NOTE: params
        # mole fraction [array]
        z_i = params['mole_fraction']
        z_i_comp = params['mole_fraction_comp']
        # pressure [Pa]
        P = params['pressure']
        # vapor pressure calculation
        VaPr_comp = params['vapor_pressure']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']

        # SECTION: activity coefficient
        # NOTE: calculate
        AcCo_i = self.__activity_coefficient(
            activity_model,
            activity,
            z_i_comp,
            T,
            activity_inputs=activity_inputs
        )

        # NOTE: calculate vapor-pressure
        # vapor pressure [Pa]
        VaPr_i = np.zeros(self.component_num)

        # looping over components
        for i, component in enumerate(self.components):
            # vapor pressure [?]
            eq_ = VaPr_comp[component]['value']
            args_ = VaPr_comp[component]['args']
            # update args
            args_['T'] = T
            # cal
            res_ = eq_.cal(**args_)
            # extract
            res_value_ = res_['value']
            res_unit_ = res_['unit']
            # convert to Pa
            unit_block_ = f"{res_unit_} => Pa"
            VaPr_ = pycuc.to(res_value_, unit_block_)

            # save
            VaPr_i[i] = VaPr_

        # NOTE: bubble pressure [Pa]
        BuPr = np.dot(z_i*AcCo_i, VaPr_i)

        # NOTE: loss
        loss = abs((P/BuPr) - 1)

        return loss

    def _DT(self, params, **kwargs):
        '''
        The `dew-point temperature` (DT) calculation determines the temperature at which the first drop of liquid condenses when a vapor mixture is cooled at a constant pressure. It identifies the temperature at which vapor will start to condense.

        Parameters
        ----------
        params : dict
            Dictionary containing the following:
            - zi : vapor mole fraction (yi)
            - P: pressure [Pa]
        kwargs : dict
            additional parameters for the calculation
            - `guess_temperature`: initial guess temperature [K], default is 295 K
            - `minimum_temperature`: minimum temperature [K], default is 200 K
            - `maximum_temperature`: maximum temperature [K], default is 1000 K
            - `max_iter`: maximum number of iterations for the calculation, default is 500
            - `tolerance`: tolerance for the calculation, default is 1e-6
            - `epsilon`: small value to avoid numerical issues, default is 1e-8

        Returns
        -------
        res : dict
            Dictionary containing the following:
            - dew_temperature: dew temperature [K]
            - pressure: pressure [Pa]
            - feed_mole_fraction: vapor mole fraction (yi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_pressure: vapor pressure [Pa]
            - activity_coefficient: activity coefficient [dimensionless]
            - K_ratio: K-ratio [dimensionless]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi).
        then zi = yi.
        - Computed Information: Temperature (T) and mole fraction in the liquid phase (xi).

        The solution is obtained by the following steps:
            1. Choose an initial guess for the dew-point temperature (Tg).
            2. Calculate the vapor pressure of each component at the guessed temperature (Tg).
            3. Calculate the dew pressure (Pb) using `Raoult's law` or other methods.
            4. Check if the calculated dew pressure (Pb) equals the given pressure (P).
            5. Calculate the liquid mole fraction (xi) using `Raoult's law` or other methods.
        '''
        try:
            # SECTION: data
            # NOTE: params
            # mole fraction [array]
            z_i = params['mole_fraction']
            # pressure [Pa]
            P = params['pressure']
            # equilibrium model
            eq_model = params['equilibrium_model']
            # fugacity model
            fugacity_model = params['fugacity_model']
            # activity model
            activity_model = params['activity_model']
            # solver method
            solver_method = params['solver_method']

            # SECTION: vapor pressure calculation
            # NOTE: vapor pressure equation [Pa]
            VaPr_comp = params['vapor_pressure']

            # NOTE: kwargs
            # temperature guess [K]
            T_g0 = kwargs.get('guess_temperature', 295)

            # SECTION: activity coefficient
            # NOTE: mole fraction dict
            # convert to dict
            z_i_comp = self.__mole_fraction_comp(z_i)

            # NOTE: init model
            # init NRTL model
            # activity = Activity(
            #     datasource=self._datasource,
            #     equationsource=self._equationsource
            # )

            activity = self.Activity_

            # activity inputs
            activity_inputs = kwargs.get('activity_inputs', {})

            # activity coefficient [dimensionless]
            AcCo_i = np.ones(self.component_num)

            # SECTION: optimization
            # pressure [Pa]
            P_value = P['value']

            # guess values
            # set max_iter
            max_iter = kwargs.get('max_iter', 500)
            # set tolerance
            tolerance = kwargs.get('tolerance', 1e-6)
            # eps
            eps = kwargs.get('eps', 1e-8)
            # guess temperature [K]
            T_g0 = kwargs.get('guess_temperature', 295)
            T_min = kwargs.get('minimum_temperature', 200)
            T_max = kwargs.get('maximum_temperature', 1000)

            # params
            _params = {
                'equilibrium_model': eq_model,
                'mole_fraction': z_i,
                'mole_fraction_comp': z_i_comp,
                'pressure': P_value,
                'vapor_pressure': VaPr_comp,
                'activity_model': activity_model,
                'activity': activity,
                'activity_inputs': activity_inputs
            }

            # NOTE: dew temperature [K]
            T = None

            # SECTION: solution
            # check solver method
            if solver_method == 'fsolve' and eq_model == 'raoult':
                # ! fsolve and raoult
                # bubble temperature calculation
                _res0 = optimize.fsolve(
                    self.fDT,
                    T_g0,
                    args=(_params),
                    full_output=True)

                # extract
                T, infodict, ier, msg = _res0

                # check if root found
                if ier != 1:
                    raise Exception(f'root not found!, {msg}')

                # check
                if len(T) == 1:
                    T = float(T[0])

            elif solver_method == 'root' and eq_model == 'raoult':
                # ! root and raoult
                _res0 = optimize.root(
                    self.fDT, T_g0, args=(_params,))

                # check
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]

            elif solver_method == 'least-squares' and eq_model == 'raoult':
                # ! least-squares and raoult
                # NOTE: Bounds
                lower_bounds = T_min
                upper_bounds = T_max - eps

                bounds = (lower_bounds, upper_bounds)

                _res0 = optimize.least_squares(
                    self.fDT,
                    T_g0,
                    args=(_params,),
                    bounds=bounds,)

                # check
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]

            elif solver_method == 'least-squares' and eq_model == 'modified-raoult':
                # ! least-squares and modified-raoult
                # NOTE: initial guess
                N = self.component_num
                # Initial guess: beta, x1, x2, ..., x_{N-1}
                x0 = [T_g0] + [1.0 / N] * (N - 1)

                # NOTE: Bounds
                lower_bounds = [T_min] + [eps] * \
                    (N - 1)         # β and each x_i ≥ eps
                upper_bounds = [T_max - eps] + [1.0] * \
                    (N - 1)   # β ≤ 1−eps, x_i ≤ 1

                bounds = (lower_bounds, upper_bounds)

                _res0 = optimize.least_squares(
                    self.fDT3,
                    x0,
                    args=(_params,),
                    bounds=bounds,)

                # check
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]
                # liquid mole fraction
                x_i = np.zeros_like(z_i)
                x_i[:-1] = _res0.x[1:]
                x_i[-1] = 1 - np.sum(x_i[:-1])
                # to array
                x_i = np.array(x_i)
                x_i_ = [float(i) for i in x_i]
                # comp
                x_i_comp = self.__mole_fraction_comp(x_i_)

                # SECTION: calculate activity coefficient
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T,
                    **kwargs
                )

            elif solver_method == 'fsolve' and eq_model == 'modified-raoult':
                # ! fsolve or root and modified-raoult
                # NOTE: initial guess
                N = self.component_num
                # Initial guess: beta, x1, x2, ..., x_{N-1}
                x0 = [T_g0] + [1.0 / N] * (N - 1)

                # NOTE: bubble temperature calculation
                _res0 = optimize.fsolve(
                    self.fDT3,
                    x0,
                    args=(_params),
                    full_output=True)

                # extract
                x, infodict, ier, msg = _res0

                # check if root found
                if ier != 1:
                    raise Exception(f'root not found!, {msg}')

                # extract
                T = x[0]
                # liquid mole fraction
                x_i = np.zeros_like(z_i)
                x_i[:-1] = x[1:]
                x_i[-1] = 1 - np.sum(x_i[:-1])
                # to array
                x_i = np.array(x_i)
                x_i_ = [float(i) for i in x_i]
                # comp
                x_i_comp = self.__mole_fraction_comp(x_i_)

                # SECTION: calculate activity coefficient
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T,
                    **kwargs
                )

            elif solver_method == 'root' and eq_model == 'modified-raoult':
                # ! root and modified-raoult
                # NOTE: initial guess
                N = self.component_num
                # Initial guess: beta, x1, x2, ..., x_{N-1}
                x0 = [T_g0] + [1.0 / N] * (N - 1)

                # NOTE: bubble temperature calculation
                _res0 = optimize.root(
                    self.fDT3,
                    x0,
                    args=(_params),)

                # check if root found
                if not _res0.success:
                    raise Exception(f'root not found!, {_res0.message}')

                # extract
                T = _res0.x[0]
                # liquid mole fraction
                x_i = np.zeros_like(z_i)
                x_i[:-1] = _res0.x[1:]
                x_i[-1] = 1 - np.sum(x_i[:-1])
                # to array
                x_i = np.array(x_i)
                x_i_ = [float(i) for i in x_i]
                # comp
                x_i_comp = self.__mole_fraction_comp(x_i_)

                # SECTION: calculate activity coefficient
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T,
                    **kwargs
                )

            else:
                raise Exception('solver method not found!')

            # NOTE: vapor pressure [Pa]
            # at T (Tg)
            VaPr_i = np.zeros(self.component_num)

            # looping over components
            for i, component in enumerate(self.components):
                # vapor pressure [?]
                eq_ = VaPr_comp[component]['value']
                args_ = VaPr_comp[component]['args']
                # update args
                args_['T'] = T
                # cal
                res_ = eq_.cal(**args_)
                # extract
                res_value_ = res_['value']
                res_unit_ = res_['unit']
                # convert to Pa
                unit_block_ = f"{res_unit_} => Pa"
                VaPr_ = pycuc.to(res_value_, unit_block_)

                # save
                VaPr_i[i] = VaPr_

            # NOTE: liquid mole fraction
            if eq_model == 'raoult':
                x_i = np.zeros(self.component_num)
                # looping over components
                for i in range(self.component_num):
                    x_i[i] = (z_i[i]*P_value)/(VaPr_i[i]*AcCo_i[i])

            # NOTE: k-ratio
            K_i = np.multiply(x_i, 1/z_i)

            # NOTE: results
            res = {
                "dew_temperature": {
                    "value": float(T),
                    "unit": "K"
                },
                "pressure": {
                    "value": P_value,
                    "unit": "Pa"
                },
                "feed_mole_fraction": z_i,
                "liquid_mole_fraction": x_i,
                "vapor_mole_fraction": z_i,
                "mole_fraction_sum": {
                    "xi": float(np.sum(x_i)),
                    "yi": float(np.sum(z_i)),
                    "zi": float(np.sum(z_i))
                },
                "vapor_pressure": {
                    "value": VaPr_i,
                    "unit": "Pa"
                },
                "K_ratio": {
                    "value": K_i,
                    "unit": "dimensionless"
                },
                "activity_coefficient": {
                    "value": AcCo_i,
                    "unit": "dimensionless"
                }
            }

            # res
            return res
        except Exception as e:
            raise Exception(f'dew temperature calculation failed! {e}')

    def fDT(self, x, params) -> float:
        '''
        dew temperature function, find `temperature at dew point` using raoult'law assumption.

        Parameters
        ----------
        x : array-like
            Guess temperature [K]
        params : dict
            Dictionary containing the following:
            - mole_fraction : liquid mole fraction (xi)
            - pressure : pressure [Pa]
            - vapor_pressure : vapor pressure equation [Pa]

        Returns
        -------
        loss : float
            Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P).

        Notes
        -----
        The dew temperature function calculates the dew pressure using the given temperature and compares it with the provided pressure (P). It returns the absolute difference as the loss function value.

        - Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi).
        then zi = yi.
        - Computed Information: Temperature (T)
        - The mole fraction in the liquid phase (xi) is calculated later.
        '''
        # NOTE: temperature (loop)
        T = x[0]

        # NOTE: params
        # mole fraction [array]
        y_i = params['mole_fraction']
        # pressure [Pa]
        P = params['pressure']
        # vapor pressure calculation
        VaPr_comp = params['vapor_pressure']

        # NOTE: calculate vapor-pressure
        # vapor pressure [Pa]
        VaPr_i = np.zeros(self.component_num)

        # looping over components
        for i, component in enumerate(self.components):
            # vapor pressure [?]
            eq_ = VaPr_comp[component]['value']
            args_ = VaPr_comp[component]['args']
            # update args
            args_['T'] = T
            # cal
            res_ = eq_.cal(**args_)
            # extract
            res_value_ = res_['value']
            res_unit_ = res_['unit']
            # convert to Pa
            unit_block_ = f"{res_unit_} => Pa"
            VaPr_ = pycuc.to(res_value_, unit_block_)

            # save
            VaPr_i[i] = VaPr_

        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # NOTE: dew pressure [Pa]
        DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

        # SECTION: loss function
        loss = abs(float(P/DePr) - 1)

        return loss

    def fDT1(self, x, params) -> float:
        '''
        dew temperature function using modified raoult'law assumption, but does not return the liquid mole fraction.

        Parameters
        ----------
        x : array-like
            Guess temperature [K]
        params : dict
            Dictionary containing the following:
            - mole_fraction : liquid mole fraction (xi)
            - pressure : pressure [Pa]
            - vapor_pressure : vapor pressure equation [Pa]
        '''
        # NOTE: temperature (loop)
        T = x[0]

        # NOTE: params
        # equilibrium model
        eq_model = params['equilibrium_model']
        # mole fraction [array]
        y_i = params['mole_fraction']
        y_i_comp = params['mole_fraction_comp']
        # pressure [Pa]
        P = params['pressure']
        # vapor pressure calculation
        VaPr_comp = params['vapor_pressure']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']
        # max iteration
        max_iter = params.get('max_iter', 500)
        # tolerance
        tolerance = params.get('tolerance', 1e-6)

        # NOTE: calculate vapor-pressure
        # vapor pressure [Pa]
        VaPr_i = np.zeros(self.component_num)

        # looping over components
        for i, component in enumerate(self.components):
            # vapor pressure [?]
            eq_ = VaPr_comp[component]['value']
            args_ = VaPr_comp[component]['args']
            # update args
            args_['T'] = T
            # cal
            res_ = eq_.cal(**args_)
            # extract
            res_value_ = res_['value']
            res_unit_ = res_['unit']
            # convert to Pa
            unit_block_ = f"{res_unit_} => Pa"
            VaPr_ = pycuc.to(res_value_, unit_block_)

            # save
            VaPr_i[i] = VaPr_

        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # NOTE: dew pressure [Pa]
        DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

        # SECTION: equilibrium model
        if eq_model == 'modified-raoult':
            # ! Modified Raoult's law
            # set x loop
            x_i_comp = y_i_comp.copy()
            x_i = y_i.copy()

            for _ in range(max_iter):
                # NOTE: calculate activity coefficient
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T,
                    activity_inputs=activity_inputs
                )

                # NOTE: liquid mole fraction
                x_i_new = (y_i*DePr)/(VaPr_i*AcCo_i)
                x_i_new /= np.sum(x_i_new)

                # NOTE: check convergence
                if np.all(np.abs(x_i - x_i_new) < tolerance):
                    break

                # NOTE: dew pressure [Pa]
                DePr_new = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

                # update loop
                x_i = x_i_new
                x_i_comp = self.__mole_fraction_comp(x_i)
                DePr = DePr_new

        # SECTION: loss function
        loss = abs(float(P/DePr) - 1)

        return loss

    def fDT2(self, x, params) -> List[float]:
        '''
        dew temperature function using modified raoult'law assumption, the loss function is the absolute difference between the calculated dew pressure and the given pressure (P).

        Parameters
        ----------
        x : array-like
            Guess temperature [K]
        params : dict
            Dictionary containing the following:
            - mole_fraction : liquid mole fraction (xi)
            - pressure : pressure [Pa]
            - vapor_pressure : vapor pressure equation [Pa]

        Returns
        -------
        loss : list[float]
            Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P).
            - Single scalar value.
        '''
        # NOTE: temperature (loop)
        T = x[0]

        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        x_i[-1] = 1 - np.sum(x_i[:-1])

        if np.any(x_i <= 1e-8) or T < 200 or T > 1000:
            return [1e6]

        # NOTE: params
        # equilibrium model
        eq_model = params['equilibrium_model']
        # mole fraction [array]
        y_i = params['mole_fraction']
        y_i_comp = params['mole_fraction_comp']
        # pressure [Pa]
        P = params['pressure']
        # vapor pressure calculation
        VaPr_comp = params['vapor_pressure']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']
        # max iteration
        max_iter = params.get('max_iter', 500)
        # tolerance
        tolerance = params.get('tolerance', 1e-6)

        # NOTE: calculate vapor-pressure
        # vapor pressure [Pa]
        VaPr_i = np.zeros(self.component_num)

        # looping over components
        for i, component in enumerate(self.components):
            # vapor pressure [?]
            eq_ = VaPr_comp[component]['value']
            args_ = VaPr_comp[component]['args']
            # update args
            args_['T'] = T
            # cal
            res_ = eq_.cal(**args_)
            # extract
            res_value_ = res_['value']
            res_unit_ = res_['unit']
            # convert to Pa
            unit_block_ = f"{res_unit_} => Pa"
            VaPr_ = pycuc.to(res_value_, unit_block_)

            # save
            VaPr_i[i] = VaPr_

        # SECTION: calculate activity coefficient
        # set liquid mole fraction
        x_i_comp = self.__mole_fraction_comp(x_i)

        # calculate activity coefficient
        AcCo_i = self.__activity_coefficient(
            activity_model,
            activity,
            x_i_comp,
            T,
            activity_inputs=activity_inputs
        )

        # NOTE: dew pressure [Pa]
        DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

        # SECTION: loss function
        loss = abs((float(P)/float(DePr)) - 1)

        return [loss]

    def fDT3(self, x, params) -> float:
        '''
        dew temperature function using modified raoult'law assumption, a system of equations is solved to find the dew temperature and liquid mole fraction.

        Parameters
        ----------
        x : array-like
            Guess temperature [K]
        params : dict
            Dictionary containing the following:
            - mole_fraction : liquid mole fraction (xi)
            - pressure : pressure [Pa]
            - vapor_pressure : vapor pressure equation [Pa]

        Returns
        -------
        residuals : array-like
            Residuals of the system of equations.
            - Vector of N residuals, each component balance enforced
        '''
        # NOTE: temperature (loop)
        T = x[0]

        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        x_i[-1] = 1 - np.sum(x_i[:-1])

        if np.any(x_i <= 1e-8) or T < 200 or T > 1000:
            return [1e6]

        # NOTE: params
        # equilibrium model
        eq_model = params['equilibrium_model']
        # mole fraction [array]
        y_i = params['mole_fraction']
        y_i_comp = params['mole_fraction_comp']
        # pressure [Pa]
        P = params['pressure']
        # vapor pressure calculation
        VaPr_comp = params['vapor_pressure']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']
        # max iteration
        max_iter = params.get('max_iter', 500)
        # tolerance
        tolerance = params.get('tolerance', 1e-6)

        # NOTE: calculate vapor-pressure
        # vapor pressure [Pa]
        VaPr_i = np.zeros(self.component_num)

        # looping over components
        for i, component in enumerate(self.components):
            # vapor pressure [?]
            eq_ = VaPr_comp[component]['value']
            args_ = VaPr_comp[component]['args']
            # update args
            args_['T'] = T
            # cal
            res_ = eq_.cal(**args_)
            # extract
            res_value_ = res_['value']
            res_unit_ = res_['unit']
            # convert to Pa
            unit_block_ = f"{res_unit_} => Pa"
            VaPr_ = pycuc.to(res_value_, unit_block_)

            # save
            VaPr_i[i] = VaPr_

        # SECTION: calculate activity coefficient
        # set liquid mole fraction
        x_i_comp = self.__mole_fraction_comp(x_i)

        # calculate activity coefficient
        AcCo_i = self.__activity_coefficient(
            activity_model,
            activity,
            x_i_comp,
            T,
            activity_inputs=activity_inputs
        )

        # SECTION: loss function
        residuals = (y_i * P) / (AcCo_i * VaPr_i) - x_i

        return residuals

    def __cal_bubble_pressure(self, xi: np.ndarray, VaPr: np.ndarray) -> float:
        '''
        Calculate bubble pressure according to Raoult's law.

        Parameters
        ----------
        xi : array-like
            Liquid mole fraction of each component in the mixture.
        VaPr : array-like
            Vapor pressure of each component in the mixture.

        Returns
        -------
        BuPr : float
            Bubble pressure of the mixture.
        '''
        # bubble pressure
        BuPr = np.dot(xi, VaPr)

        # res
        return BuPr

    def __cal_dew_pressure(self, yi: np.ndarray, VaPr: np.ndarray) -> float:
        '''
        Calculate dew pressure according to Raoult's law.

        Parameters
        ----------
        yi : array-like
            Vapor mole fraction of each component in the mixture.
        VaPr : array-like
            Vapor pressure of each component in the mixture.

        Returns
        -------
        DePr : float
            Dew pressure of the mixture.
        '''
        # dew-point pressure
        DePr = 1/np.dot(yi, 1/VaPr)

        # res
        return DePr

    def _flash_checker(self,
                       z_i: List[float],
                       Pf: float,
                       Tf: float,
                       VaPr_comp: Dict) -> bool:
        '''
        Check if the flash occurs at the given pressure and temperature according to the bubble and dew pressures of the mixture; according to Raoult's law.

        Parameters
        ----------
        z_i : list
            Feed mole fraction of each component in the mixture.
        Pf : float
            Flash pressure [Pa]
        Tf : float
            Flash temperature [K]
        VaPr_comp : dict
            Dictionary containing the vapor pressure equations for each component.

        Returns
        -------
        bool
            True if the pressure and temperature are valid, False otherwise.

        Notes
        -----
        Assumption are as:

        - The system is in equilibrium.
        - The vapor pressure of each component is calculated using the provided equations.
        - The bubble pressure is calculated using Raoult's law.
        - The dew pressure is calculated using Raoult's law.
        '''
        try:
            # NOTE: mole fraction [array]
            z_i_ = np.array(z_i)

            # NOTE: vapor pressure [Pa]
            VaPr_i = np.zeros(self.component_num)

            # looping over components
            for i, component in enumerate(self.components):
                # vapor pressure [?]
                eq_ = VaPr_comp[component]['value']
                args_ = VaPr_comp[component]['args']
                # update args
                args_['T'] = Tf
                # cal
                res_ = eq_.cal(**args_)
                # extract
                res_value_ = res_['value']
                res_unit_ = res_['unit']
                # convert to Pa
                unit_block_ = f"{res_unit_} => Pa"
                VaPr_ = pycuc.to(res_value_, unit_block_)
                # save
                VaPr_i[i] = VaPr_

            # NOTE: calculate bubble pressure
            BuPr = self.__cal_bubble_pressure(z_i_, VaPr_i)
            # NOTE: calculate dew pressure
            DePr = self.__cal_dew_pressure(z_i_, VaPr_i)

            # NOTE: check if the given pressure and temperature are within the valid range
            if Pf < BuPr and Pf > DePr:
                # if the pressure is between the bubble and dew pressures, return False
                return True
            else:
                # if the pressure is outside the valid range, return False
                return False
        except Exception as e:
            # if any error occurs, return False
            raise Exception(f"Error in flash checker: {e}")

    def _IFL(self, params, **kwargs):
        '''
        The `isothermal-flash` (IFL) calculation This calculation determines the vapor and liquid phase compositions and amounts at a specified temperature and pressure. The system is "flashed" isothermally, meaning the temperature is kept constant while the phase behavior is calculated for a mixture.

        Parameters
        ----------
        params : dict
            Dictionary containing the following:
            - zi : feed mole fraction (zi)
            - P: pressure [Pa]
            - T: temperature [K]
        kwargs : dict
            additional parameters for the calculation
            - `activity_inputs`: additional inputs for the activity model, default is {}
            - `eps`: small value to avoid numerical issues, default is 1e-8
            - `guess_vapor_fraction`: vapor fraction (VF) [dimensionless], default is 0.5

        Returns
        -------
        res : dict
            Dictionary containing the following:
            - vapor_to_liquid_ratio: V/F ratio [dimensionless]
            - liquid_to_vapor_ratio: L/F ratio [dimensionless]
            - feed_mole_fraction: liquid mole fraction (zi)
            - liquid_mole_fraction: liquid mole fraction (xi)
            - vapor_mole_fraction: vapor mole fraction (yi)
            - vapor_pressure: vapor pressure [Pa]
            - K_ratio: K-ratio [dimensionless]
            - temperature: temperature [K]
            - pressure: pressure [Pa]

        Notes
        -----
        The summary of the calculation is as follows:

        - Known Information: Temperature (T), pressure (P), and mole fraction of the components in the liquid phase (zi).
        - Computed Information: Vapor-to-liquid ratio (V/F), liquid mole fraction (xi), vapor mole fraction (yi), and liquid-to-vapor ratio (L/F).

        The solution is obtained by the following steps:
            1. Calculate the K ratio (Ki) using Raoult's law.
            2. Choose an initial guess for the vapor-to-liquid ratio (V/F).
            3. Solve a system of non-linear equations to find the V/F ratio that satisfies the mass balance equations.
            4. Calculate the liquid and vapor mole fractions (xi and yi) using the V/F ratio and K ratio.
            5. Calculate the liquid-to-vapor ratio (L/F) from the V/F ratio.
        '''
        try:
            # SECTION: data
            # NOTE: params
            # mole fraction [array]
            z_i = params['mole_fraction']
            # pressure [Pa]
            P = params['pressure']
            # temperature [K]
            T = params['temperature']
            # equilibrium model
            eq_model = params['equilibrium_model']
            # fugacity model
            fugacity_model = params['fugacity_model']
            # activity model
            activity_model = params['activity_model']
            # solver method
            solver_method = params['solver_method']

            # SECTION: vapor pressure calculation
            # NOTE: vapor pressure equation [Pa]
            VaPr_comp = params['vapor_pressure']

            # activity coefficient
            AcCo_i = np.ones(self.component_num)

            # NOTE: kwarg
            # activity inputs
            activity_inputs = kwargs.get('activity_inputs', {})

            # NOTE: set values
            # pressure [Pa]
            P_value = P['value']
            # temperature [K]
            T_value = T['value']

            # NOTE: ki ratio (Raoult's law)
            # pressure and temperature are constant (Raoutl's law)
            K_i = np.zeros(self.component_num)

            # vapor pressure [Pa]
            VaPr_i = np.zeros(self.component_num)

            # looping over components
            for i, component in enumerate(self.components):
                # vapor pressure [?]
                eq_ = VaPr_comp[component]['value']
                args_ = VaPr_comp[component]['args']
                # update args
                args_['T'] = T_value
                # cal
                res_ = eq_.cal(**args_)
                # extract
                res_value_ = res_['value']
                res_unit_ = res_['unit']
                # convert to Pa
                unit_block_ = f"{res_unit_} => Pa"
                VaPr_ = pycuc.to(res_value_, unit_block_)
                # set
                VaPr_i[i] = VaPr_
                K_i[i] = VaPr_/P_value

            # SECTION: activity model
            # init
            # activity = Activity(
            #     datasource=self._datasource,
            #     equationsource=self._equationsource
            # )

            activity = self.Activity_

            # SECTION: optimization
            # NOTE: guess values
            VF0_g = kwargs.get('guess_vapor_fraction', 0.5)
            # eps
            eps = kwargs.get('eps', 1e-8)

            # NOTE: mole fraction dict
            z_i_comp = self.__mole_fraction_comp(z_i)

            # NOTE: set params
            _params = {
                'equilibrium_model': eq_model,
                'mole_fraction': z_i,
                'mole_fraction_comp': z_i_comp,
                'temperature': T_value,
                'pressure': P_value,
                'K_ratio': K_i,
                'activity_model': activity_model,
                'activity': activity,
                'activity_inputs': activity_inputs,
            }

            # NOTE: solver message
            solver_message = None

            # NOTE: check solver method
            if solver_method == 'least_squares':
                # ! least-squares
                # NOTE: initial guess
                N = self.component_num
                # Initial guess: beta, x1, x2, ..., x_{N-1}
                x0 = [VF0_g] + [1.0 / N] * (N - 1)

                # NOTE: Bounds
                lower_bounds = [eps] + [eps] * \
                    (N - 1)         # β and each x_i ≥ eps
                upper_bounds = [1.0 - eps] + [1.0] * \
                    (N - 1)   # β ≤ 1−eps, x_i ≤ 1

                bounds = (lower_bounds, upper_bounds)

                # NOTE: set function
                if eq_model == 'raoult':
                    fn = self.fIFL
                elif eq_model == 'modified-raoult':
                    fn = self.fIFL1
                else:
                    raise Exception('equilibrium model not found!')

                # V/F
                _res = optimize.least_squares(
                    fn,
                    x0,
                    args=(_params,),
                    bounds=bounds)

                # check if root found
                if _res.success is False:
                    raise Exception(f'root not found!, {_res.message}')

                # solver message
                solver_message = _res.message

                # -> result analysis
                V_F_ratio = _res.x[0]
                x_i = np.zeros_like(z_i)
                x_i[:-1] = _res.x[1:]
                x_i[-1] = 1 - np.sum(x_i[:-1])
                # to array
                x_i = np.array(x_i)
                # comp
                x_i_comp = self.__mole_fraction_comp(x_i)

                # SECTION: calculate activity coefficient
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T_value,
                    **kwargs
                )

            elif solver_method == 'minimize':
                # ! minimize
                # NOTE: initial guess
                N = self.component_num
                # Initial guess: beta, x1, x2, ..., x_{N-1}
                x0 = [VF0_g] + [1.0 / N] * (N - 1)

                # NOTE: Bounds
                # individual bounds (min, max)
                bounds = [(eps, 1.0 - eps)] + [(eps, 1.0)
                                               for _ in range(N - 1)]

                # V/F
                _res = optimize.minimize(
                    self.fIFL2,
                    x0=x0,
                    args=(_params,),
                    constraints=self.flash_constraints(
                        z_i,
                        K_i,
                        T_value,
                        P_value,
                        VaPr_i,
                        equilibrium_model=eq_model,
                        activity_model=activity_model,
                        activity=activity,
                        **kwargs),
                    bounds=bounds,
                )

                # check if root found
                if _res.success is False:
                    raise Exception(f'root not found! {_res.message}')

                # set
                solver_message = _res.message

                # -> result analysis
                V_F_ratio = _res.x[0]
                x_i = np.zeros_like(z_i)
                x_i[:-1] = _res.x[1:]
                x_i[-1] = 1 - np.sum(x_i[:-1])
                # to array
                x_i = np.array(x_i)
                # comp
                x_i_comp = self.__mole_fraction_comp(x_i)

                # SECTION: calculate activity coefficient
                # NOTE: calculate
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_i_comp,
                    T_value,
                    **kwargs
                )

            else:
                raise Exception('solver method not found!')

            # SECTION: liquid/vapor mole fraction
            xy_ = self.xy_flash(V_F_ratio, z_i, K_i, AcCo_i)
            # set
            xi = xy_['liquid']
            yi = xy_['vapor']

            # NOTE: calculate L/F
            L_F_ratio = 1 - V_F_ratio

            # SECTION: calculate activity coefficient
            # liquid mole fraction
            xi_comp = self.__mole_fraction_comp(xi)

            # NOTE: check model
            AcCo_i = self.__activity_coefficient(
                activity_model,
                activity,
                xi_comp,
                T_value,
                **kwargs
            )

            # NOTE: results
            res = {
                "V_F_ratio": {
                    "value": float(V_F_ratio),
                    "unit": "dimensionless"
                },
                "L_F_ratio": {
                    "value": float(L_F_ratio),
                    "unit": "dimensionless"
                },
                "feed_mole_fraction": z_i,
                "liquid_mole_fraction": xi,
                "vapor_mole_fraction": yi,
                "mole_fraction_sum": {
                    "xi": float(np.sum(xi)),
                    "yi": float(np.sum(yi)),
                    "zi": float(np.sum(z_i))
                },
                "vapor_pressure": {
                    "value": VaPr_i,
                    "unit": "Pa"
                },
                "K_ratio": {
                    "value": K_i,
                    "unit": "dimensionless"
                },
                "temperature": {
                    "value": T_value,
                    "unit": "K"
                },
                "pressure": {
                    "value": P_value,
                    "unit": "Pa"
                },
                "activity_coefficient": {
                    "value": AcCo_i,
                    "unit": "dimensionless"
                },
                "solver_message": solver_message,
            }

            # res
            return res
        except Exception as e:
            raise Exception(f"flash isothermal failed! {e}")

    def fIFL(self, x, params):
        '''
        Flash isothermal function, according to `Raoult's law`.

        Parameters
        ----------
        x : array-like
            VF : vapor-feed ratio [dimensionless]
            x_i : liquid mole fraction [dimensionless]
        params : tuple
            Tuple containing the following:
            - zi : feed mole fraction [dimensionless]
            - Ki : K ratio [dimensionless]
        '''
        # NOTE: extract variables
        # vapor fraction
        VF = x[0]
        # liquid mole fraction
        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        # last x_i
        x_i[-1] = 1 - np.sum(x_i[:-1])

        # NOTE: params
        # feed mole fraction
        z_i = params['mole_fraction']
        # K ratio (P*/P) [dimensionless]
        K_i = params['K_ratio']

        # NOTE: activity coefficient
        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # NOTE: vapor mole fraction
        K_i = AcCo_i * K_i
        y_i = K_i * x_i

        # SECTION: check optimization region
        # Avoid division by zero or invalid regions
        if VF <= 0 or VF >= 1:
            return 1e6  # Return a large value to indicate invalid region

        # mole fraction
        # NOTE: check if x_i is valid
        if np.any(x_i <= 1e-8) or np.any(x_i >= 1 - 1e-8):
            return 1e6

        # NOTE: sum x_i = 1
        if np.abs(np.sum(x_i) - 1) > 1e-8:
            return 1e6

        # NOTE: function
        eqs = []

        # looping over components
        for i in range(self.component_num - 1):
            eq = z_i[i] - ((1 - VF) * x_i[i] + VF * y_i[i])
            eqs.append(eq)

        # Closure relations
        # eqs.append(np.sum(x_i) - 1)   # sum x_i = 1
        # eqs.append(np.sum(y_i) - 1)   # sum y_i = 1
        # x-y < eps
        eqs.append(np.sum(x_i) - np.sum(y_i) - 1e-8)  # sum x_i = sum y_i

        return eqs

    def fIFL1(self, x, params):
        '''
        Flash isothermal function, according to modified Raoult's law.

        Parameters
        ----------
        x : array-like
            VF : vapor-feed ratio [dimensionless]
            x_i : liquid mole fraction [dimensionless]
        params : tuple
            Tuple containing the following:
            - zi : feed mole fraction [dimensionless]
            - Ki : K ratio [dimensionless]
        '''
        # NOTE: extract variables
        # vapor fraction
        VF = x[0]
        # liquid mole fraction
        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        # last x_i
        x_i[-1] = 1 - np.sum(x_i[:-1])

        # NOTE: params
        # feed mole fraction
        z_i = params['mole_fraction']
        # equilibrium model
        eq_model = params['equilibrium_model']
        # temperature [K]
        T = params['temperature']
        # pressure [Pa]
        P = params['pressure']
        # K ratio (P*/P) [dimensionless]
        K_i = params['K_ratio']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']

        # NOTE: activity coefficient
        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # SECTION: equilibrium model
        # ! Modified Raoult's law
        # set x loop
        x_i_comp = self.__mole_fraction_comp(x_i)

        # NOTE: calculate activity coefficient
        # NOTE: check model
        if activity_model == 'NRTL':
            # calculate activity
            res_ = activity.NRTL(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        elif activity_model == 'UNIQUAC':
            # calculate activity
            res_ = activity.UNIQUAC(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        else:
            # equals unity for ideal solution
            AcCo_i = np.ones(self.component_num)

        # NOTE: vapor mole fraction
        K_i = AcCo_i * K_i
        y_i = K_i * x_i

        # SECTION: check optimization region
        # Avoid division by zero or invalid regions
        if VF <= 0 or VF >= 1:
            return 1e6  # Return a large value to indicate invalid region

        # mole fraction
        # NOTE: check if x_i is valid
        if np.any(x_i <= 1e-8) or np.any(x_i >= 1 - 1e-8):
            return 1e6

        # NOTE: sum x_i = 1
        if np.abs(np.sum(x_i) - 1) > 1e-8:
            return 1e6

        # NOTE: function
        eqs = []

        # looping over components
        for i in range(self.component_num - 1):
            eq = z_i[i] - ((1 - VF) * x_i[i] + VF * y_i[i])
            eqs.append(eq)

        # Closure relations
        # eqs.append(np.sum(x_i) - 1)   # sum x_i = 1
        # eqs.append(np.sum(y_i) - 1)   # sum y_i = 1
        # x-y < eps
        eqs.append(np.sum(x_i) - np.sum(y_i) - 1e-8)  # sum x_i = sum y_i

        return eqs

    def fIFL2(self, x, params):
        '''
        Flash isothermal function (nonlinear equations), according to Raoult's law.

        Parameters
        ----------
        x : array-like
            VF : vapor-feed ratio [dimensionless]
            x_i : liquid mole fraction [dimensionless]
        params : tuple
            Tuple containing the following:
            - zi: feed mole fraction (zi)
            - Ki: K ratio of each component in the mixture
        '''
        # NOTE: extract variables
        # vapor fraction
        VF = x[0]
        # liquid mole fraction
        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        # last x_i
        x_i[-1] = 1 - np.sum(x_i[:-1])

        # NOTE: params
        # equilibrium model
        eq_model = params['equilibrium_model']
        # feed mole fraction
        z_i = params['mole_fraction']
        # temperature [K]
        T = params['temperature']
        # pressure [Pa]
        P = params['pressure']
        # K ratio (P*/P) [dimensionless]
        K_i = params['K_ratio']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']

        # NOTE: activity coefficient
        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # SECTION: equilibrium model
        if eq_model == 'modified-raoult':
            # ! Modified Raoult's law
            # set x loop
            x_i_comp = self.__mole_fraction_comp(x_i)

            # NOTE: calculate activity coefficient
            # NOTE: check model
            if activity_model == 'NRTL':
                # calculate activity
                res_ = activity.NRTL(
                    self.components,
                    x_i_comp,
                    T,
                    activity_inputs=activity_inputs)
                # extract
                AcCo_i = res_['value']
            elif activity_model == 'UNIQUAC':
                # calculate activity
                res_ = activity.UNIQUAC(
                    self.components,
                    x_i_comp,
                    T,
                    activity_inputs=activity_inputs)
                # extract
                AcCo_i = res_['value']
            else:
                # equals unity for ideal solution
                AcCo_i = np.ones(self.component_num)

        # NOTE: update K_i
        K_i = AcCo_i * K_i
        # NOTE: vapor mole fraction
        y_i = K_i * x_i

        # SECTION: system of nonlinear equations (NLE)
        # Residuals of component balances
        residuals = z_i - ((1 - VF) * x_i + VF * y_i)

        # x_i sum = 1
        residuals = np.append(residuals, np.sum(x_i) - 1)

        # Objective: minimize sum of squared residuals
        return np.sum(residuals**2)

    def fIFL3(self, x, params):
        '''
        Flash isothermal function (nonlinear equations), according to Raoult's law.

        Parameters
        ----------
        x : array-like
            VF : vapor-feed ratio [dimensionless]
            x_i : liquid mole fraction [dimensionless]
        params : tuple
            Tuple containing the following:
            - zi: feed mole fraction (zi)
            - Ki: K ratio of each component in the mixture
        '''
        # NOTE: extract variables
        # vapor fraction
        VF = x[0]
        # liquid mole fraction
        x_i = np.zeros(self.component_num)
        x_i[:-1] = x[1:]
        # last x_i
        x_i[-1] = 1 - np.sum(x_i[:-1])

        # NOTE: params
        # equilibrium model
        eq_model = params['equilibrium_model']
        # feed mole fraction
        z_i = params['mole_fraction']
        # temperature [K]
        T = params['temperature']
        # pressure [Pa]
        P = params['pressure']
        # K ratio (P*/P) [dimensionless]
        K_i = params['K_ratio']
        # activity model
        activity_model = params['activity_model']
        # activity
        activity: Activity = params['activity']
        # activity inputs
        activity_inputs = params['activity_inputs']

        # NOTE: activity coefficient
        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

        # SECTION: equilibrium model
        if eq_model == 'modified-raoult':
            # ! Modified Raoult's law
            # set x loop
            x_i_comp = self.__mole_fraction_comp(x_i)

            # NOTE: calculate activity coefficient
            # NOTE: check model
            if activity_model == 'NRTL':
                # calculate activity
                res_ = activity.NRTL(
                    self.components,
                    x_i_comp,
                    T,
                    activity_inputs=activity_inputs)
                # extract
                AcCo_i = res_['value']
            elif activity_model == 'UNIQUAC':
                # calculate activity
                res_ = activity.UNIQUAC(
                    self.components,
                    x_i_comp,
                    T,
                    activity_inputs=activity_inputs)
                # extract
                AcCo_i = res_['value']
            else:
                # equals unity for ideal solution
                AcCo_i = np.ones(self.component_num)

        # NOTE: update K_i
        K_i = AcCo_i * K_i
        # NOTE: vapor mole fraction
        y_i = K_i * x_i

        # SECTION: system of nonlinear equations (NLE)
        # Residuals of component balances
        residuals = z_i - ((1 - VF) * x_i + VF * y_i)

        # x_i sum = 1
        residuals = np.append(residuals, np.sum(x_i) - 1)

        # Objective: minimize sum of squared residuals
        return np.sum(residuals**2)

    def flash_constraints(self,
                          z: np.ndarray,
                          K: np.ndarray,
                          T: float,
                          P: float,
                          P_sat: np.ndarray,
                          equilibrium_model: Literal[
                              'raoult', 'modified-raoult'
                          ] = 'raoult',
                          activity_model: Literal[
                              'NRTL', 'UNIQUAC'
                          ] = 'NRTL',
                          activity: Optional[Activity] = None,
                          **kwargs):
        '''
        Constraints for the flash calculation.

        Parameters
        ----------
        z : array-like
            Feed mole fraction of each component in the mixture.
        K : array-like
            K ratio of each component in the mixture.
        T : float
            Flash temperature [K].
        P : float
            Flash pressure [Pa].
        P_sat : array-like
            Saturation pressure of each component in the mixture.
        equilibrium_model : str, optional
            Equilibrium model to be used. The default is 'raoult'.
        activity_model : str, optional
            Activity model to be used. The default is 'NRTL'.
        activity : Activity, optional
            Activity object for calculating activity coefficients.
            The default is None.
        kwargs : dict, optional
            Additional parameters for the calculation.
            - `activity_inputs`: additional inputs for the activity model.
        '''
        N = len(z)

        def liquid_sum(vars):
            """Equality constraint for liquid phase."""
            x = np.zeros(N)
            x[:-1] = vars[1:]
            x[-1] = 1.0 - np.sum(x[:-1])

            # res
            return np.sum(x) - 1

        def vapor_sum(vars):
            """Equality constraint for vapor phase."""
            x = np.zeros(N)
            x[:-1] = vars[1:]
            x[-1] = 1.0 - np.sum(x[:-1])

            # set x comp
            x_comp = self.__mole_fraction_comp(x)

            # check model
            if equilibrium_model == 'raoult':
                # ! Raoult's law
                # calculate y
                y = K * x
            elif equilibrium_model == 'modified-raoult':
                # ! Modified Raoult's law
                # check
                if activity:
                    # calculate activity coefficient
                    AcCo_i = self.__activity_coefficient(
                        activity_model,
                        activity,
                        x_comp,
                        T,
                        **kwargs)
                else:
                    raise ValueError("Activity model not provided.")

                # calculate y
                y = AcCo_i * K * x
            else:
                raise ValueError("Invalid equilibrium model.")

            # res
            return np.sum(y) - 1

        return [
            # {'type': 'eq', 'fun': liquid_sum},
            {'type': 'eq', 'fun': vapor_sum},
        ]

    def xy_flash(self,
                 V_F_ratio: float,
                 z_i: np.ndarray,
                 K_i: np.ndarray,
                 AcCo_i: Optional[
                     np.ndarray
                 ] = None
                 ) -> Dict[str, np.ndarray]:
        '''
        Calculate liquid/vapor mole fraction (xi, yi) using V/F ratio and K ratio.

        Parameters
        ----------
        V_F_ratio : float
            Vapor-to-liquid ratio (V/F).
        zi : array-like
            Feed mole fraction of each component in the mixture.
        Ki : array-like
            K ratio of each component in the mixture.
        AcCo_i : array-like
            Activity coefficient of each component in the mixture.

        Returns
        -------
        dict
            Dictionary containing the following:
            - liquid: liquid mole fraction (xi) [dimensionless]
            - vapor: vapor mole fraction (yi) [dimensionless]
        '''
        try:
            # check
            if AcCo_i is None:
                # equals unity for ideal solution
                AcCo_i = np.ones(self.component_num)

            # init
            x_i = np.zeros(self.component_num)
            y_i = np.zeros(self.component_num)

            # liquid/vapor mole fraction
            for i in range(self.component_num):
                x_i[i] = (z_i[i])/(1+(V_F_ratio)*(K_i[i]*AcCo_i[i]-1))
                y_i[i] = K_i[i]*AcCo_i[i]*x_i[i]

            # res
            return {
                'liquid': x_i,
                'vapor': y_i
            }
        except Exception as e:
            raise Exception(f"Error in xy_flash calculation: {e}")

component_num property

Get the number of components in the system.

components property

Get the components of the system.

__activity_coefficient(activity_model, activity, x_i_comp, T, **kwargs)

Calculate activity coefficient using NRTL or UNIQUAC model.

Parameters

activity_model : str Activity model to be used ('NRTL' or 'UNIQUAC'). activity : Activity Activity object for calculating activity coefficients. x_i_comp : dict Dictionary of mole fractions of each component in the mixture. T : float Temperature [K]. kwargs : dict Additional parameters for the calculation.

Returns

AcCo_i : array-like Activity coefficient of each component in the mixture.

Source code in pyThermoFlash/docs/equilibria.py

def __activity_coefficient(self,
                           activity_model: Literal['NRTL', 'UNIQUAC'],
                           activity: Activity,
                           x_i_comp: Dict[str, float],
                           T: float,
                           **kwargs
                           ) -> np.ndarray:
    '''
    Calculate activity coefficient using NRTL or UNIQUAC model.

    Parameters
    ----------
    activity_model : str
        Activity model to be used ('NRTL' or 'UNIQUAC').
    activity : Activity
        Activity object for calculating activity coefficients.
    x_i_comp : dict
        Dictionary of mole fractions of each component in the mixture.
    T : float
        Temperature [K].
    kwargs : dict
        Additional parameters for the calculation.

    Returns
    -------
    AcCo_i : array-like
        Activity coefficient of each component in the mixture.
    '''
    try:
        # NOTE: check model
        if activity_model == 'NRTL':
            # calculate activity
            res_ = activity.NRTL(
                self.components,
                x_i_comp,
                T,
                **kwargs)
            # extract
            AcCo_i = res_['value']
        elif activity_model == 'UNIQUAC':
            # calculate activity
            res_ = activity.UNIQUAC(
                self.components,
                x_i_comp,
                T,
                **kwargs)
            # extract
            AcCo_i = res_['value']
        else:
            # equals unity for ideal solution
            AcCo_i = np.ones(self.component_num)

        # res
        return AcCo_i
    except Exception as e:
        raise Exception(f'activity coefficient calculation failed! {e}')

__cal_bubble_pressure(xi, VaPr)

Calculate bubble pressure according to Raoult's law.

Parameters

xi : array-like Liquid mole fraction of each component in the mixture. VaPr : array-like Vapor pressure of each component in the mixture.

Returns

BuPr : float Bubble pressure of the mixture.

Source code in pyThermoFlash/docs/equilibria.py

def __cal_bubble_pressure(self, xi: np.ndarray, VaPr: np.ndarray) -> float:
    '''
    Calculate bubble pressure according to Raoult's law.

    Parameters
    ----------
    xi : array-like
        Liquid mole fraction of each component in the mixture.
    VaPr : array-like
        Vapor pressure of each component in the mixture.

    Returns
    -------
    BuPr : float
        Bubble pressure of the mixture.
    '''
    # bubble pressure
    BuPr = np.dot(xi, VaPr)

    # res
    return BuPr

__cal_dew_pressure(yi, VaPr)

Calculate dew pressure according to Raoult's law.

Parameters

yi : array-like Vapor mole fraction of each component in the mixture. VaPr : array-like Vapor pressure of each component in the mixture.

Returns

DePr : float Dew pressure of the mixture.

Source code in pyThermoFlash/docs/equilibria.py

def __cal_dew_pressure(self, yi: np.ndarray, VaPr: np.ndarray) -> float:
    '''
    Calculate dew pressure according to Raoult's law.

    Parameters
    ----------
    yi : array-like
        Vapor mole fraction of each component in the mixture.
    VaPr : array-like
        Vapor pressure of each component in the mixture.

    Returns
    -------
    DePr : float
        Dew pressure of the mixture.
    '''
    # dew-point pressure
    DePr = 1/np.dot(yi, 1/VaPr)

    # res
    return DePr

__check_activity_coefficients(activity_model, activity, z_i_comp, T_value, **kwargs)

Check if activity coefficients are provided by the user.

Returns

AcCo_i : np.ndarray Activity coefficients as a numpy array.

Source code in pyThermoFlash/docs/equilibria.py

def __check_activity_coefficients(self,
                                  activity_model: Literal['NRTL', 'UNIQUAC'],
                                  activity: Activity,
                                  z_i_comp: Dict[str, float],
                                  T_value: float,
                                  **kwargs
                                  ) -> np.ndarray:
    '''
    Check if activity coefficients are provided by the user.

    Returns
    -------
    AcCo_i : np.ndarray
        Activity coefficients as a numpy array.
    '''
    try:
        # SECTION: activity coefficient
        # ! check kwargs
        activity_coefficients_ext = kwargs.get(
            'activity_coefficients', None)

        # check calculated activity coefficients
        if activity_coefficients_ext:
            # set
            AcCo_i_ = activity_coefficients_ext

            # ! set activity coefficients
            AcCo_i = self.set_calculated_activity_coefficient(
                AcCo_i_)
        else:
            # NOTE: calculate
            AcCo_i = self.__activity_coefficient(
                activity_model,
                activity,
                z_i_comp,
                T_value,
                **kwargs
            )

        # res
        return AcCo_i
    except Exception as e:
        raise Exception(f'check activity coefficients failed! {e}')

__init__(components, datasource=None, equationsource=None, **kwargs)

Initialize the Equilibria class.

Source code in pyThermoFlash/docs/equilibria.py

def __init__(self,
             components: List[str],
             datasource: Optional[Dict] = None,
             equationsource: Optional[Dict] = None,
             **kwargs):
    '''Initialize the Equilibria class.'''
    # set
    components_ = [component.strip() for component in components]
    self.__components = components_
    self.__comp_num = len(components_)

    # model source
    self._datasource = datasource
    self._equationsource = equationsource

    # NOTE: init activity model
    self.Activity_ = Activity(
        datasource=self._datasource,
        equationsource=self._equationsource
    )

__mole_fraction_comp(z_i)

Convert mole fraction list to dictionary.

Parameters

z_i : list | np.ndarray List of mole fractions of each component in the mixture.

Returns

dict Dictionary of mole fractions.

Source code in pyThermoFlash/docs/equilibria.py

def __mole_fraction_comp(self,
                         z_i: List[float] | np.ndarray
                         ) -> Dict[str, float]:
    '''
    Convert mole fraction list to dictionary.

    Parameters
    ----------
    z_i : list | np.ndarray
        List of mole fractions of each component in the mixture.

    Returns
    -------
    dict
        Dictionary of mole fractions.
    '''
    # NOTE: check type
    if isinstance(z_i, np.ndarray):
        # convert to list
        z_i = [float(z) for z in z_i]
    # convert to dict
    return {str(self.components[i]): z_i[i]
            for i in range(self.component_num)}

fBT(x, params)

bubble temperature function

Parameters

x : array-like Guess temperature [K] params : dict Dictionary containing the following: - mole_fraction : liquid mole fraction (xi) - pressure : pressure [Pa] - vapor_pressure : vapor pressure equation [Pa]

Source code in pyThermoFlash/docs/equilibria.py

def fBT(self, x, params) -> float:
    '''
    bubble temperature function

    Parameters
    ----------
    x : array-like
        Guess temperature [K]
    params : dict
        Dictionary containing the following:
        - mole_fraction : liquid mole fraction (xi)
        - pressure : pressure [Pa]
        - vapor_pressure : vapor pressure equation [Pa]
    '''
    # NOTE: temperature (loop)
    T = x[0]

    # NOTE: params
    # mole fraction [array]
    z_i = params['mole_fraction']
    z_i_comp = params['mole_fraction_comp']
    # pressure [Pa]
    P = params['pressure']
    # vapor pressure calculation
    VaPr_comp = params['vapor_pressure']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']

    # SECTION: activity coefficient
    # NOTE: calculate
    AcCo_i = self.__activity_coefficient(
        activity_model,
        activity,
        z_i_comp,
        T,
        activity_inputs=activity_inputs
    )

    # NOTE: calculate vapor-pressure
    # vapor pressure [Pa]
    VaPr_i = np.zeros(self.component_num)

    # looping over components
    for i, component in enumerate(self.components):
        # vapor pressure [?]
        eq_ = VaPr_comp[component]['value']
        args_ = VaPr_comp[component]['args']
        # update args
        args_['T'] = T
        # cal
        res_ = eq_.cal(**args_)
        # extract
        res_value_ = res_['value']
        res_unit_ = res_['unit']
        # convert to Pa
        unit_block_ = f"{res_unit_} => Pa"
        VaPr_ = pycuc.to(res_value_, unit_block_)

        # save
        VaPr_i[i] = VaPr_

    # NOTE: bubble pressure [Pa]
    BuPr = np.dot(z_i*AcCo_i, VaPr_i)

    # NOTE: loss
    loss = abs((P/BuPr) - 1)

    return loss

fDT(x, params)

dew temperature function, find temperature at dew point using raoult'law assumption.

Parameters

x : array-like Guess temperature [K] params : dict Dictionary containing the following: - mole_fraction : liquid mole fraction (xi) - pressure : pressure [Pa] - vapor_pressure : vapor pressure equation [Pa]

Returns

loss : float Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P).

Notes

The dew temperature function calculates the dew pressure using the given temperature and compares it with the provided pressure (P). It returns the absolute difference as the loss function value.

• Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi). then zi = yi.
• Computed Information: Temperature (T)
• The mole fraction in the liquid phase (xi) is calculated later.

Source code in pyThermoFlash/docs/equilibria.py

def fDT(self, x, params) -> float:
    '''
    dew temperature function, find `temperature at dew point` using raoult'law assumption.

    Parameters
    ----------
    x : array-like
        Guess temperature [K]
    params : dict
        Dictionary containing the following:
        - mole_fraction : liquid mole fraction (xi)
        - pressure : pressure [Pa]
        - vapor_pressure : vapor pressure equation [Pa]

    Returns
    -------
    loss : float
        Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P).

    Notes
    -----
    The dew temperature function calculates the dew pressure using the given temperature and compares it with the provided pressure (P). It returns the absolute difference as the loss function value.

    - Known Information: Pressure (P) and mole fraction of the components in the vapor phase (yi).
    then zi = yi.
    - Computed Information: Temperature (T)
    - The mole fraction in the liquid phase (xi) is calculated later.
    '''
    # NOTE: temperature (loop)
    T = x[0]

    # NOTE: params
    # mole fraction [array]
    y_i = params['mole_fraction']
    # pressure [Pa]
    P = params['pressure']
    # vapor pressure calculation
    VaPr_comp = params['vapor_pressure']

    # NOTE: calculate vapor-pressure
    # vapor pressure [Pa]
    VaPr_i = np.zeros(self.component_num)

    # looping over components
    for i, component in enumerate(self.components):
        # vapor pressure [?]
        eq_ = VaPr_comp[component]['value']
        args_ = VaPr_comp[component]['args']
        # update args
        args_['T'] = T
        # cal
        res_ = eq_.cal(**args_)
        # extract
        res_value_ = res_['value']
        res_unit_ = res_['unit']
        # convert to Pa
        unit_block_ = f"{res_unit_} => Pa"
        VaPr_ = pycuc.to(res_value_, unit_block_)

        # save
        VaPr_i[i] = VaPr_

    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # NOTE: dew pressure [Pa]
    DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

    # SECTION: loss function
    loss = abs(float(P/DePr) - 1)

    return loss

fDT1(x, params)

dew temperature function using modified raoult'law assumption, but does not return the liquid mole fraction.

Parameters

x : array-like Guess temperature [K] params : dict Dictionary containing the following: - mole_fraction : liquid mole fraction (xi) - pressure : pressure [Pa] - vapor_pressure : vapor pressure equation [Pa]

Source code in pyThermoFlash/docs/equilibria.py

def fDT1(self, x, params) -> float:
    '''
    dew temperature function using modified raoult'law assumption, but does not return the liquid mole fraction.

    Parameters
    ----------
    x : array-like
        Guess temperature [K]
    params : dict
        Dictionary containing the following:
        - mole_fraction : liquid mole fraction (xi)
        - pressure : pressure [Pa]
        - vapor_pressure : vapor pressure equation [Pa]
    '''
    # NOTE: temperature (loop)
    T = x[0]

    # NOTE: params
    # equilibrium model
    eq_model = params['equilibrium_model']
    # mole fraction [array]
    y_i = params['mole_fraction']
    y_i_comp = params['mole_fraction_comp']
    # pressure [Pa]
    P = params['pressure']
    # vapor pressure calculation
    VaPr_comp = params['vapor_pressure']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']
    # max iteration
    max_iter = params.get('max_iter', 500)
    # tolerance
    tolerance = params.get('tolerance', 1e-6)

    # NOTE: calculate vapor-pressure
    # vapor pressure [Pa]
    VaPr_i = np.zeros(self.component_num)

    # looping over components
    for i, component in enumerate(self.components):
        # vapor pressure [?]
        eq_ = VaPr_comp[component]['value']
        args_ = VaPr_comp[component]['args']
        # update args
        args_['T'] = T
        # cal
        res_ = eq_.cal(**args_)
        # extract
        res_value_ = res_['value']
        res_unit_ = res_['unit']
        # convert to Pa
        unit_block_ = f"{res_unit_} => Pa"
        VaPr_ = pycuc.to(res_value_, unit_block_)

        # save
        VaPr_i[i] = VaPr_

    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # NOTE: dew pressure [Pa]
    DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

    # SECTION: equilibrium model
    if eq_model == 'modified-raoult':
        # ! Modified Raoult's law
        # set x loop
        x_i_comp = y_i_comp.copy()
        x_i = y_i.copy()

        for _ in range(max_iter):
            # NOTE: calculate activity coefficient
            AcCo_i = self.__activity_coefficient(
                activity_model,
                activity,
                x_i_comp,
                T,
                activity_inputs=activity_inputs
            )

            # NOTE: liquid mole fraction
            x_i_new = (y_i*DePr)/(VaPr_i*AcCo_i)
            x_i_new /= np.sum(x_i_new)

            # NOTE: check convergence
            if np.all(np.abs(x_i - x_i_new) < tolerance):
                break

            # NOTE: dew pressure [Pa]
            DePr_new = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

            # update loop
            x_i = x_i_new
            x_i_comp = self.__mole_fraction_comp(x_i)
            DePr = DePr_new

    # SECTION: loss function
    loss = abs(float(P/DePr) - 1)

    return loss

fDT2(x, params)

dew temperature function using modified raoult'law assumption, the loss function is the absolute difference between the calculated dew pressure and the given pressure (P).

Parameters

x : array-like Guess temperature [K] params : dict Dictionary containing the following: - mole_fraction : liquid mole fraction (xi) - pressure : pressure [Pa] - vapor_pressure : vapor pressure equation [Pa]

Returns

loss : list[float] Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P). - Single scalar value.

Source code in pyThermoFlash/docs/equilibria.py

def fDT2(self, x, params) -> List[float]:
    '''
    dew temperature function using modified raoult'law assumption, the loss function is the absolute difference between the calculated dew pressure and the given pressure (P).

    Parameters
    ----------
    x : array-like
        Guess temperature [K]
    params : dict
        Dictionary containing the following:
        - mole_fraction : liquid mole fraction (xi)
        - pressure : pressure [Pa]
        - vapor_pressure : vapor pressure equation [Pa]

    Returns
    -------
    loss : list[float]
        Loss function value, which is the absolute difference between the calculated dew pressure and the given pressure (P).
        - Single scalar value.
    '''
    # NOTE: temperature (loop)
    T = x[0]

    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    x_i[-1] = 1 - np.sum(x_i[:-1])

    if np.any(x_i <= 1e-8) or T < 200 or T > 1000:
        return [1e6]

    # NOTE: params
    # equilibrium model
    eq_model = params['equilibrium_model']
    # mole fraction [array]
    y_i = params['mole_fraction']
    y_i_comp = params['mole_fraction_comp']
    # pressure [Pa]
    P = params['pressure']
    # vapor pressure calculation
    VaPr_comp = params['vapor_pressure']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']
    # max iteration
    max_iter = params.get('max_iter', 500)
    # tolerance
    tolerance = params.get('tolerance', 1e-6)

    # NOTE: calculate vapor-pressure
    # vapor pressure [Pa]
    VaPr_i = np.zeros(self.component_num)

    # looping over components
    for i, component in enumerate(self.components):
        # vapor pressure [?]
        eq_ = VaPr_comp[component]['value']
        args_ = VaPr_comp[component]['args']
        # update args
        args_['T'] = T
        # cal
        res_ = eq_.cal(**args_)
        # extract
        res_value_ = res_['value']
        res_unit_ = res_['unit']
        # convert to Pa
        unit_block_ = f"{res_unit_} => Pa"
        VaPr_ = pycuc.to(res_value_, unit_block_)

        # save
        VaPr_i[i] = VaPr_

    # SECTION: calculate activity coefficient
    # set liquid mole fraction
    x_i_comp = self.__mole_fraction_comp(x_i)

    # calculate activity coefficient
    AcCo_i = self.__activity_coefficient(
        activity_model,
        activity,
        x_i_comp,
        T,
        activity_inputs=activity_inputs
    )

    # NOTE: dew pressure [Pa]
    DePr = 1/np.dot(y_i, 1/(VaPr_i*AcCo_i))

    # SECTION: loss function
    loss = abs((float(P)/float(DePr)) - 1)

    return [loss]

fDT3(x, params)

dew temperature function using modified raoult'law assumption, a system of equations is solved to find the dew temperature and liquid mole fraction.

Parameters

x : array-like Guess temperature [K] params : dict Dictionary containing the following: - mole_fraction : liquid mole fraction (xi) - pressure : pressure [Pa] - vapor_pressure : vapor pressure equation [Pa]

Returns

residuals : array-like Residuals of the system of equations. - Vector of N residuals, each component balance enforced

Source code in pyThermoFlash/docs/equilibria.py

def fDT3(self, x, params) -> float:
    '''
    dew temperature function using modified raoult'law assumption, a system of equations is solved to find the dew temperature and liquid mole fraction.

    Parameters
    ----------
    x : array-like
        Guess temperature [K]
    params : dict
        Dictionary containing the following:
        - mole_fraction : liquid mole fraction (xi)
        - pressure : pressure [Pa]
        - vapor_pressure : vapor pressure equation [Pa]

    Returns
    -------
    residuals : array-like
        Residuals of the system of equations.
        - Vector of N residuals, each component balance enforced
    '''
    # NOTE: temperature (loop)
    T = x[0]

    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    x_i[-1] = 1 - np.sum(x_i[:-1])

    if np.any(x_i <= 1e-8) or T < 200 or T > 1000:
        return [1e6]

    # NOTE: params
    # equilibrium model
    eq_model = params['equilibrium_model']
    # mole fraction [array]
    y_i = params['mole_fraction']
    y_i_comp = params['mole_fraction_comp']
    # pressure [Pa]
    P = params['pressure']
    # vapor pressure calculation
    VaPr_comp = params['vapor_pressure']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']
    # max iteration
    max_iter = params.get('max_iter', 500)
    # tolerance
    tolerance = params.get('tolerance', 1e-6)

    # NOTE: calculate vapor-pressure
    # vapor pressure [Pa]
    VaPr_i = np.zeros(self.component_num)

    # looping over components
    for i, component in enumerate(self.components):
        # vapor pressure [?]
        eq_ = VaPr_comp[component]['value']
        args_ = VaPr_comp[component]['args']
        # update args
        args_['T'] = T
        # cal
        res_ = eq_.cal(**args_)
        # extract
        res_value_ = res_['value']
        res_unit_ = res_['unit']
        # convert to Pa
        unit_block_ = f"{res_unit_} => Pa"
        VaPr_ = pycuc.to(res_value_, unit_block_)

        # save
        VaPr_i[i] = VaPr_

    # SECTION: calculate activity coefficient
    # set liquid mole fraction
    x_i_comp = self.__mole_fraction_comp(x_i)

    # calculate activity coefficient
    AcCo_i = self.__activity_coefficient(
        activity_model,
        activity,
        x_i_comp,
        T,
        activity_inputs=activity_inputs
    )

    # SECTION: loss function
    residuals = (y_i * P) / (AcCo_i * VaPr_i) - x_i

    return residuals

fIFL(x, params)

Flash isothermal function, according to Raoult's law.

Parameters

x : array-like VF : vapor-feed ratio [dimensionless] x_i : liquid mole fraction [dimensionless] params : tuple Tuple containing the following: - zi : feed mole fraction [dimensionless] - Ki : K ratio [dimensionless]

Source code in pyThermoFlash/docs/equilibria.py

def fIFL(self, x, params):
    '''
    Flash isothermal function, according to `Raoult's law`.

    Parameters
    ----------
    x : array-like
        VF : vapor-feed ratio [dimensionless]
        x_i : liquid mole fraction [dimensionless]
    params : tuple
        Tuple containing the following:
        - zi : feed mole fraction [dimensionless]
        - Ki : K ratio [dimensionless]
    '''
    # NOTE: extract variables
    # vapor fraction
    VF = x[0]
    # liquid mole fraction
    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    # last x_i
    x_i[-1] = 1 - np.sum(x_i[:-1])

    # NOTE: params
    # feed mole fraction
    z_i = params['mole_fraction']
    # K ratio (P*/P) [dimensionless]
    K_i = params['K_ratio']

    # NOTE: activity coefficient
    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # NOTE: vapor mole fraction
    K_i = AcCo_i * K_i
    y_i = K_i * x_i

    # SECTION: check optimization region
    # Avoid division by zero or invalid regions
    if VF <= 0 or VF >= 1:
        return 1e6  # Return a large value to indicate invalid region

    # mole fraction
    # NOTE: check if x_i is valid
    if np.any(x_i <= 1e-8) or np.any(x_i >= 1 - 1e-8):
        return 1e6

    # NOTE: sum x_i = 1
    if np.abs(np.sum(x_i) - 1) > 1e-8:
        return 1e6

    # NOTE: function
    eqs = []

    # looping over components
    for i in range(self.component_num - 1):
        eq = z_i[i] - ((1 - VF) * x_i[i] + VF * y_i[i])
        eqs.append(eq)

    # Closure relations
    # eqs.append(np.sum(x_i) - 1)   # sum x_i = 1
    # eqs.append(np.sum(y_i) - 1)   # sum y_i = 1
    # x-y < eps
    eqs.append(np.sum(x_i) - np.sum(y_i) - 1e-8)  # sum x_i = sum y_i

    return eqs

fIFL1(x, params)

Flash isothermal function, according to modified Raoult's law.

Parameters

x : array-like VF : vapor-feed ratio [dimensionless] x_i : liquid mole fraction [dimensionless] params : tuple Tuple containing the following: - zi : feed mole fraction [dimensionless] - Ki : K ratio [dimensionless]

Source code in pyThermoFlash/docs/equilibria.py

def fIFL1(self, x, params):
    '''
    Flash isothermal function, according to modified Raoult's law.

    Parameters
    ----------
    x : array-like
        VF : vapor-feed ratio [dimensionless]
        x_i : liquid mole fraction [dimensionless]
    params : tuple
        Tuple containing the following:
        - zi : feed mole fraction [dimensionless]
        - Ki : K ratio [dimensionless]
    '''
    # NOTE: extract variables
    # vapor fraction
    VF = x[0]
    # liquid mole fraction
    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    # last x_i
    x_i[-1] = 1 - np.sum(x_i[:-1])

    # NOTE: params
    # feed mole fraction
    z_i = params['mole_fraction']
    # equilibrium model
    eq_model = params['equilibrium_model']
    # temperature [K]
    T = params['temperature']
    # pressure [Pa]
    P = params['pressure']
    # K ratio (P*/P) [dimensionless]
    K_i = params['K_ratio']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']

    # NOTE: activity coefficient
    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # SECTION: equilibrium model
    # ! Modified Raoult's law
    # set x loop
    x_i_comp = self.__mole_fraction_comp(x_i)

    # NOTE: calculate activity coefficient
    # NOTE: check model
    if activity_model == 'NRTL':
        # calculate activity
        res_ = activity.NRTL(
            self.components,
            x_i_comp,
            T,
            activity_inputs=activity_inputs)
        # extract
        AcCo_i = res_['value']
    elif activity_model == 'UNIQUAC':
        # calculate activity
        res_ = activity.UNIQUAC(
            self.components,
            x_i_comp,
            T,
            activity_inputs=activity_inputs)
        # extract
        AcCo_i = res_['value']
    else:
        # equals unity for ideal solution
        AcCo_i = np.ones(self.component_num)

    # NOTE: vapor mole fraction
    K_i = AcCo_i * K_i
    y_i = K_i * x_i

    # SECTION: check optimization region
    # Avoid division by zero or invalid regions
    if VF <= 0 or VF >= 1:
        return 1e6  # Return a large value to indicate invalid region

    # mole fraction
    # NOTE: check if x_i is valid
    if np.any(x_i <= 1e-8) or np.any(x_i >= 1 - 1e-8):
        return 1e6

    # NOTE: sum x_i = 1
    if np.abs(np.sum(x_i) - 1) > 1e-8:
        return 1e6

    # NOTE: function
    eqs = []

    # looping over components
    for i in range(self.component_num - 1):
        eq = z_i[i] - ((1 - VF) * x_i[i] + VF * y_i[i])
        eqs.append(eq)

    # Closure relations
    # eqs.append(np.sum(x_i) - 1)   # sum x_i = 1
    # eqs.append(np.sum(y_i) - 1)   # sum y_i = 1
    # x-y < eps
    eqs.append(np.sum(x_i) - np.sum(y_i) - 1e-8)  # sum x_i = sum y_i

    return eqs

fIFL2(x, params)

Flash isothermal function (nonlinear equations), according to Raoult's law.

Parameters

x : array-like VF : vapor-feed ratio [dimensionless] x_i : liquid mole fraction [dimensionless] params : tuple Tuple containing the following: - zi: feed mole fraction (zi) - Ki: K ratio of each component in the mixture

Source code in pyThermoFlash/docs/equilibria.py

def fIFL2(self, x, params):
    '''
    Flash isothermal function (nonlinear equations), according to Raoult's law.

    Parameters
    ----------
    x : array-like
        VF : vapor-feed ratio [dimensionless]
        x_i : liquid mole fraction [dimensionless]
    params : tuple
        Tuple containing the following:
        - zi: feed mole fraction (zi)
        - Ki: K ratio of each component in the mixture
    '''
    # NOTE: extract variables
    # vapor fraction
    VF = x[0]
    # liquid mole fraction
    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    # last x_i
    x_i[-1] = 1 - np.sum(x_i[:-1])

    # NOTE: params
    # equilibrium model
    eq_model = params['equilibrium_model']
    # feed mole fraction
    z_i = params['mole_fraction']
    # temperature [K]
    T = params['temperature']
    # pressure [Pa]
    P = params['pressure']
    # K ratio (P*/P) [dimensionless]
    K_i = params['K_ratio']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']

    # NOTE: activity coefficient
    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # SECTION: equilibrium model
    if eq_model == 'modified-raoult':
        # ! Modified Raoult's law
        # set x loop
        x_i_comp = self.__mole_fraction_comp(x_i)

        # NOTE: calculate activity coefficient
        # NOTE: check model
        if activity_model == 'NRTL':
            # calculate activity
            res_ = activity.NRTL(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        elif activity_model == 'UNIQUAC':
            # calculate activity
            res_ = activity.UNIQUAC(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        else:
            # equals unity for ideal solution
            AcCo_i = np.ones(self.component_num)

    # NOTE: update K_i
    K_i = AcCo_i * K_i
    # NOTE: vapor mole fraction
    y_i = K_i * x_i

    # SECTION: system of nonlinear equations (NLE)
    # Residuals of component balances
    residuals = z_i - ((1 - VF) * x_i + VF * y_i)

    # x_i sum = 1
    residuals = np.append(residuals, np.sum(x_i) - 1)

    # Objective: minimize sum of squared residuals
    return np.sum(residuals**2)

fIFL3(x, params)

Flash isothermal function (nonlinear equations), according to Raoult's law.

Parameters

x : array-like VF : vapor-feed ratio [dimensionless] x_i : liquid mole fraction [dimensionless] params : tuple Tuple containing the following: - zi: feed mole fraction (zi) - Ki: K ratio of each component in the mixture

Source code in pyThermoFlash/docs/equilibria.py

def fIFL3(self, x, params):
    '''
    Flash isothermal function (nonlinear equations), according to Raoult's law.

    Parameters
    ----------
    x : array-like
        VF : vapor-feed ratio [dimensionless]
        x_i : liquid mole fraction [dimensionless]
    params : tuple
        Tuple containing the following:
        - zi: feed mole fraction (zi)
        - Ki: K ratio of each component in the mixture
    '''
    # NOTE: extract variables
    # vapor fraction
    VF = x[0]
    # liquid mole fraction
    x_i = np.zeros(self.component_num)
    x_i[:-1] = x[1:]
    # last x_i
    x_i[-1] = 1 - np.sum(x_i[:-1])

    # NOTE: params
    # equilibrium model
    eq_model = params['equilibrium_model']
    # feed mole fraction
    z_i = params['mole_fraction']
    # temperature [K]
    T = params['temperature']
    # pressure [Pa]
    P = params['pressure']
    # K ratio (P*/P) [dimensionless]
    K_i = params['K_ratio']
    # activity model
    activity_model = params['activity_model']
    # activity
    activity: Activity = params['activity']
    # activity inputs
    activity_inputs = params['activity_inputs']

    # NOTE: activity coefficient
    # equals unity for ideal solution
    AcCo_i = np.ones(self.component_num)

    # SECTION: equilibrium model
    if eq_model == 'modified-raoult':
        # ! Modified Raoult's law
        # set x loop
        x_i_comp = self.__mole_fraction_comp(x_i)

        # NOTE: calculate activity coefficient
        # NOTE: check model
        if activity_model == 'NRTL':
            # calculate activity
            res_ = activity.NRTL(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        elif activity_model == 'UNIQUAC':
            # calculate activity
            res_ = activity.UNIQUAC(
                self.components,
                x_i_comp,
                T,
                activity_inputs=activity_inputs)
            # extract
            AcCo_i = res_['value']
        else:
            # equals unity for ideal solution
            AcCo_i = np.ones(self.component_num)

    # NOTE: update K_i
    K_i = AcCo_i * K_i
    # NOTE: vapor mole fraction
    y_i = K_i * x_i

    # SECTION: system of nonlinear equations (NLE)
    # Residuals of component balances
    residuals = z_i - ((1 - VF) * x_i + VF * y_i)

    # x_i sum = 1
    residuals = np.append(residuals, np.sum(x_i) - 1)

    # Objective: minimize sum of squared residuals
    return np.sum(residuals**2)

flash_constraints(z, K, T, P, P_sat, equilibrium_model='raoult', activity_model='NRTL', activity=None, **kwargs)

Constraints for the flash calculation.

Parameters

z : array-like Feed mole fraction of each component in the mixture. K : array-like K ratio of each component in the mixture. T : float Flash temperature [K]. P : float Flash pressure [Pa]. P_sat : array-like Saturation pressure of each component in the mixture. equilibrium_model : str, optional Equilibrium model to be used. The default is 'raoult'. activity_model : str, optional Activity model to be used. The default is 'NRTL'. activity : Activity, optional Activity object for calculating activity coefficients. The default is None. kwargs : dict, optional Additional parameters for the calculation. - activity_inputs: additional inputs for the activity model.

Source code in pyThermoFlash/docs/equilibria.py

def flash_constraints(self,
                      z: np.ndarray,
                      K: np.ndarray,
                      T: float,
                      P: float,
                      P_sat: np.ndarray,
                      equilibrium_model: Literal[
                          'raoult', 'modified-raoult'
                      ] = 'raoult',
                      activity_model: Literal[
                          'NRTL', 'UNIQUAC'
                      ] = 'NRTL',
                      activity: Optional[Activity] = None,
                      **kwargs):
    '''
    Constraints for the flash calculation.

    Parameters
    ----------
    z : array-like
        Feed mole fraction of each component in the mixture.
    K : array-like
        K ratio of each component in the mixture.
    T : float
        Flash temperature [K].
    P : float
        Flash pressure [Pa].
    P_sat : array-like
        Saturation pressure of each component in the mixture.
    equilibrium_model : str, optional
        Equilibrium model to be used. The default is 'raoult'.
    activity_model : str, optional
        Activity model to be used. The default is 'NRTL'.
    activity : Activity, optional
        Activity object for calculating activity coefficients.
        The default is None.
    kwargs : dict, optional
        Additional parameters for the calculation.
        - `activity_inputs`: additional inputs for the activity model.
    '''
    N = len(z)

    def liquid_sum(vars):
        """Equality constraint for liquid phase."""
        x = np.zeros(N)
        x[:-1] = vars[1:]
        x[-1] = 1.0 - np.sum(x[:-1])

        # res
        return np.sum(x) - 1

    def vapor_sum(vars):
        """Equality constraint for vapor phase."""
        x = np.zeros(N)
        x[:-1] = vars[1:]
        x[-1] = 1.0 - np.sum(x[:-1])

        # set x comp
        x_comp = self.__mole_fraction_comp(x)

        # check model
        if equilibrium_model == 'raoult':
            # ! Raoult's law
            # calculate y
            y = K * x
        elif equilibrium_model == 'modified-raoult':
            # ! Modified Raoult's law
            # check
            if activity:
                # calculate activity coefficient
                AcCo_i = self.__activity_coefficient(
                    activity_model,
                    activity,
                    x_comp,
                    T,
                    **kwargs)
            else:
                raise ValueError("Activity model not provided.")

            # calculate y
            y = AcCo_i * K * x
        else:
            raise ValueError("Invalid equilibrium model.")

        # res
        return np.sum(y) - 1

    return [
        # {'type': 'eq', 'fun': liquid_sum},
        {'type': 'eq', 'fun': vapor_sum},
    ]

set_calculated_activity_coefficient(AcCo_i_cal)

Set calculated activity coefficients, provided by the user.

Parameters

AcCo_i_cal : dict | np.ndarray | list Activity coefficients for each component in the mixture. If a dict is provided, it should contain component names as keys and their corresponding activity coefficients as values. If a numpy array or list is provided, it should contain activity coefficients in the same order as the components.

Returns

AcCo_i : np.ndarray Activity coefficients as a numpy array.

Source code in pyThermoFlash/docs/equilibria.py

def set_calculated_activity_coefficient(self,
                                        AcCo_i_cal:
                                        Dict[str, float] | np.ndarray | list,
                                        ) -> np.ndarray:
    """
    Set calculated activity coefficients, provided by the user.

    Parameters
    ----------
    AcCo_i_cal : dict | np.ndarray | list
        Activity coefficients for each component in the mixture.
        If a dict is provided, it should contain component names as keys
        and their corresponding activity coefficients as values.
        If a numpy array or list is provided, it should contain activity
        coefficients in the same order as the components.

    Returns
    -------
    AcCo_i : np.ndarray
        Activity coefficients as a numpy array.
    """
    try:
        # NOTE: check type
        if isinstance(AcCo_i_cal, dict):
            # convert to array based on components
            AcCo_i = np.zeros(self.component_num)
            for i, component in enumerate(self.components):
                if component in AcCo_i_cal:
                    AcCo_i[i] = AcCo_i_cal[component]
                else:
                    raise Exception(
                        f"activity_coefficients for {component} not found!")
        elif isinstance(AcCo_i_cal, np.ndarray):
            AcCo_i = AcCo_i_cal
        elif isinstance(AcCo_i_cal, list):
            AcCo_i = np.array(AcCo_i_cal)
        else:
            raise Exception(
                'activity_coefficients must be a dict, list or numpy array!')
        return AcCo_i
    except Exception as e:
        raise Exception(f'set activity coefficient failed! {e}')

xy_flash(V_F_ratio, z_i, K_i, AcCo_i=None)

Calculate liquid/vapor mole fraction (xi, yi) using V/F ratio and K ratio.

Parameters

V_F_ratio : float Vapor-to-liquid ratio (V/F). zi : array-like Feed mole fraction of each component in the mixture. Ki : array-like K ratio of each component in the mixture. AcCo_i : array-like Activity coefficient of each component in the mixture.

Returns

dict Dictionary containing the following: - liquid: liquid mole fraction (xi) [dimensionless] - vapor: vapor mole fraction (yi) [dimensionless]

Source code in pyThermoFlash/docs/equilibria.py

def xy_flash(self,
             V_F_ratio: float,
             z_i: np.ndarray,
             K_i: np.ndarray,
             AcCo_i: Optional[
                 np.ndarray
             ] = None
             ) -> Dict[str, np.ndarray]:
    '''
    Calculate liquid/vapor mole fraction (xi, yi) using V/F ratio and K ratio.

    Parameters
    ----------
    V_F_ratio : float
        Vapor-to-liquid ratio (V/F).
    zi : array-like
        Feed mole fraction of each component in the mixture.
    Ki : array-like
        K ratio of each component in the mixture.
    AcCo_i : array-like
        Activity coefficient of each component in the mixture.

    Returns
    -------
    dict
        Dictionary containing the following:
        - liquid: liquid mole fraction (xi) [dimensionless]
        - vapor: vapor mole fraction (yi) [dimensionless]
    '''
    try:
        # check
        if AcCo_i is None:
            # equals unity for ideal solution
            AcCo_i = np.ones(self.component_num)

        # init
        x_i = np.zeros(self.component_num)
        y_i = np.zeros(self.component_num)

        # liquid/vapor mole fraction
        for i in range(self.component_num):
            x_i[i] = (z_i[i])/(1+(V_F_ratio)*(K_i[i]*AcCo_i[i]-1))
            y_i[i] = K_i[i]*AcCo_i[i]*x_i[i]

        # res
        return {
            'liquid': x_i,
            'vapor': y_i
        }
    except Exception as e:
        raise Exception(f"Error in xy_flash calculation: {e}")

Activity 🔄

Activity

Source code in pyThermoFlash/docs/activity.py

class Activity:
    def __init__(self,
                 datasource: Optional[Dict[str, Any]] = None,
                 equationsource: Optional[Dict[str, Any]] = None
                 ):
        '''
        Activity model class for pyThermoFlash.

        Parameters
        ----------
        datasource : dict, optional
            Data source for the model, containing data for components and equations.
        equationsource : dict, optional
            Equation source for the model, containing equations for components.
        '''
        # set datasource
        self.datasource = datasource
        # set equationsource
        self.equationsource = equationsource

    def __repr__(self):
        return "Activity model class for pyThermoFlash"

    def NRTL(self,
             components: List[str],
             z_i_comp: Dict[str, float],
             temperature: float,
             **kwargs):
        '''
        NRTL activity model for calculating activity coefficients.

        Parameters
        ----------
        components : list
            List of component names.
        z_i_comp : dict
            Dictionary of component names and their respective mole fractions.
        temperature : float
            Temperature in Kelvin.
        kwargs : dict
            Additional parameters for the model.
            - interaction-energy-parameter : list, optional
                Interaction energy parameters for the components.
        '''
        try:
            # SECTION: check src
            # extract activity model inputs
            activity_inputs = kwargs.get('activity_inputs', None)
            nrtl_datasource = kwargs.get('NRTL', None)

            # check
            if activity_inputs is None:
                # ! check nrtl inputs in datasource
                activity_inputs = nrtl_datasource

            # check if activity_inputs is None
            if activity_inputs is None:
                raise ValueError(
                    "No valid source provided for activity model (NRTL) inputs.")

            # SECTION: check if activity_inputs is a dictionary
            if activity_inputs is not None:
                # check if activity_inputs is a dictionary
                if not isinstance(activity_inputs, dict):
                    raise ValueError(
                        "activity_inputs must be a dictionary.")
                # check if activity_inputs is empty
                if len(activity_inputs) == 0:
                    raise ValueError(
                        "activity_inputs cannot be empty.")

                # NOTE: update activity_inputs with nrtl_datasource
                if nrtl_datasource is not None:
                    activity_inputs.update(nrtl_datasource)

            # SECTION: extract activity model inputs
            # NOTE: method 1
            # ! Δg_ij, interaction energy parameter
            dg_ij_src = activity_inputs.get('dg_ij', None)
            if dg_ij_src is None:
                dg_ij_src = activity_inputs.get('dg', None)

            # NOTE: method 2
            # ! constants a, b, c, and d
            a_ij_src = activity_inputs.get('a_ij', None)
            if a_ij_src is None:
                a_ij_src = activity_inputs.get('a', None)
            b_ij_src = activity_inputs.get('b_ij', None)
            if b_ij_src is None:
                b_ij_src = activity_inputs.get('b', None)
            c_ij_src = activity_inputs.get('c_ij', None)
            if c_ij_src is None:
                c_ij_src = activity_inputs.get('c', None)
            d_ij_src = activity_inputs.get('d_ij', None)
            if d_ij_src is None:
                d_ij_src = activity_inputs.get('d', None)

            # NOTE: α_ij, non-randomness parameter
            alpha_ij_src = activity_inputs.get('alpha_ij', None)
            if alpha_ij_src is None:
                alpha_ij_src = activity_inputs.get('alpha', None)

            # SECTION: init NRTL model
            # activity model
            activity = ptm.activity(
                components=components, model_name='NRTL')
            # set
            activity_nrtl = activity.nrtl
            # check
            if activity_nrtl is None:
                raise ValueError(
                    "Failed to initialize NRTL activity model. Please check the components and model name.")

            # NOTE: check method
            tau_ij_cal_method = 0
            if dg_ij_src is None:
                # check if a_ij, b_ij, c_ij are provided
                if a_ij_src is None or b_ij_src is None or c_ij_src is None or d_ij_src is None:
                    raise ValueError(
                        "No valid source provided for interaction energy parameter (Δg_ij) or constants a, b, c, and d.")
                # set method
                tau_ij_cal_method = 2

                # ! a_ij
                if isinstance(a_ij_src, TableMatrixData):
                    a_ij = a_ij_src.mat('a', components)
                elif isinstance(a_ij_src, list):
                    a_ij = np.array(a_ij_src)
                elif isinstance(a_ij_src, np.ndarray):
                    a_ij = a_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (a_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! b_ij
                if isinstance(b_ij_src, TableMatrixData):
                    b_ij = b_ij_src.mat('b', components)
                elif isinstance(b_ij_src, list):
                    b_ij = np.array(b_ij_src)
                elif isinstance(b_ij_src, np.ndarray):
                    b_ij = b_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (b_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! c_ij
                if isinstance(c_ij_src, TableMatrixData):
                    c_ij = c_ij_src.mat('c', components)
                elif isinstance(c_ij_src, list):
                    c_ij = np.array(c_ij_src)
                elif isinstance(c_ij_src, np.ndarray):
                    c_ij = c_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (c_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! d_ij
                if isinstance(d_ij_src, TableMatrixData):
                    d_ij = d_ij_src.mat('d', components)
                elif isinstance(d_ij_src, list):
                    d_ij = np.array(d_ij_src)
                elif isinstance(d_ij_src, np.ndarray):
                    d_ij = d_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (d_ij). Must be TableMatrixData, list of lists, or numpy array.")
            elif dg_ij_src is not None:
                # use dg_ij
                if isinstance(dg_ij_src, TableMatrixData):
                    dg_ij = dg_ij_src.mat('dg', components)
                elif isinstance(dg_ij_src, list):
                    dg_ij = np.array(dg_ij_src)
                elif isinstance(dg_ij_src, np.ndarray):
                    dg_ij = dg_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (Δg_ij). Must be TableMatrixData, list of lists, or numpy array.")
                # set method
                tau_ij_cal_method = 1
            else:
                raise ValueError(
                    "No valid source provided for interaction energy parameter (Δg_ij) or constants A, B, C.")

            # SECTION: extract data
            # α_ij, non-randomness parameter
            if isinstance(alpha_ij_src, TableMatrixData):
                alpha_ij = alpha_ij_src.mat('alpha', components)
            elif isinstance(alpha_ij_src, list):
                alpha_ij = np.array(alpha_ij_src)
            elif isinstance(alpha_ij_src, np.ndarray):
                alpha_ij = alpha_ij_src
            else:
                raise ValueError(
                    "Invalid source for non-randomness parameter (α_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # NOTE: calculate the binary interaction parameter matrix (tau_ij)
            if tau_ij_cal_method == 1:
                # check
                tau_ij, _ = activity_nrtl.cal_tau_ij_M1(
                    temperature=temperature,
                    dg_ij=dg_ij
                )
            elif tau_ij_cal_method == 2:
                tau_ij, _ = activity_nrtl.cal_tau_ij_M2(
                    temperature=temperature,
                    a_ij=a_ij,
                    b_ij=b_ij,
                    c_ij=c_ij,
                    d_ij=d_ij
                )
            else:
                raise ValueError(
                    "Invalid method for calculating tau_ij. Must be 1 or 2.")

            # NOTE: nrtl inputs
            inputs_ = {
                'mole_fraction': z_i_comp,
                "tau_ij": tau_ij,
                "alpha_ij": alpha_ij
            }

            # NOTE: calculate activity
            res_, _ = activity_nrtl.cal(model_input=inputs_)

            # res format
            # res = {
            #     'property_name': 'activity coefficients',
            #     'components': components,
            #     'mole_fraction': xi,
            #     'value': AcCo_i,
            #     'unit': 1,
            #     'symbol': "AcCo_i",
            #     'message': message,
            # }

            # res
            return res_
        except Exception as e:
            raise Exception(f"Failed to calculate NRTL activity: {e}") from e

    def UNIQUAC(self,
                components: List[str],
                z_i_comp: Dict[str, float],
                temperature: float,
                **kwargs):
        '''
        UNIQUAC activity model for calculating activity coefficients.

        Parameters
        ----------
        components : list
            List of component names.
        z_i_comp : dict
            Dictionary of component names and their respective mole fractions.
        temperature : float
            Temperature in Kelvin.
        kwargs : dict
            Additional parameters for the model.
            - interaction-energy-parameter : list, optional
                Interaction energy parameters for the components.
        '''
        try:
            # SECTION: check src
            # extract activity model inputs
            activity_inputs = kwargs.get('activity_inputs', None)
            uniquac_datasource = kwargs.get('UNIQUAC', None)

            # check
            if activity_inputs is None:
                # ! check nrtl inputs in datasource
                activity_inputs = uniquac_datasource

            # check if activity_inputs is None
            if activity_inputs is None:
                raise ValueError(
                    "No valid source provided for activity model (UNIQUAC) inputs.")

            # NOTE: check if activity_inputs is a dictionary
            if activity_inputs is not None:
                # check if activity_inputs is a dictionary
                if not isinstance(activity_inputs, dict):
                    raise ValueError(
                        "activity_inputs must be a dictionary.")
                # check if activity_inputs is empty
                if len(activity_inputs) == 0:
                    raise ValueError(
                        "activity_inputs cannot be empty.")

            # NOTE: update activity_inputs with uniquac_datasource
            if uniquac_datasource is not None:
                activity_inputs.update(uniquac_datasource)

            # NOTE: method 1
            # Δg_ij, interaction energy parameter
            dU_ij_src = activity_inputs.get('dU_ij', None)
            if dU_ij_src is None:
                dU_ij_src = activity_inputs.get('dU', None)

            # NOTE: method 2
            # constants a, b, c, and d
            a_ij_src = activity_inputs.get('a_ij', None)
            if a_ij_src is None:
                a_ij_src = activity_inputs.get('a', None)
            b_ij_src = activity_inputs.get('b_ij', None)
            if b_ij_src is None:
                b_ij_src = activity_inputs.get('b', None)
            c_ij_src = activity_inputs.get('c_ij', None)
            if c_ij_src is None:
                c_ij_src = activity_inputs.get('c', None)
            d_ij_src = activity_inputs.get('d_ij', None)
            if d_ij_src is None:
                d_ij_src = activity_inputs.get('d', None)

            # NOTE: r_i, relative van der Waals volume of component i
            r_i_src = activity_inputs.get('r_i', None)
            if r_i_src is None:
                r_i_src = activity_inputs.get('r', None)

            # final check if r_i is provided
            if r_i_src is None:
                raise ValueError("No valid source provided for r_i.")

            # check if r_i is a list or numpy array
            if isinstance(r_i_src, list):
                r_i = np.array(r_i_src)
            elif isinstance(r_i_src, np.ndarray):
                r_i = r_i_src
            else:
                raise ValueError(
                    "Invalid source for r_i. Must be a list or numpy array.")

            # NOTE: q_i, relative van der Waals area of component i
            q_i_src = activity_inputs.get('q_i', None)
            if q_i_src is None:
                q_i_src = activity_inputs.get('q', None)

            # final check if q_i is provided
            if q_i_src is None:
                raise ValueError("No valid source provided for q_i.")

            # check if q_i is a list or numpy array
            if isinstance(q_i_src, list):
                q_i = np.array(q_i_src)
            elif isinstance(q_i_src, np.ndarray):
                q_i = q_i_src
            else:
                raise ValueError(
                    "Invalid source for q_i. Must be a list or numpy array.")

            # SECTION: init NRTL model
            # activity model
            activity = ptm.activity(
                components=components, model_name='NRTL')
            # set
            activity_uniquac = activity.uniquac

            # check
            if activity_uniquac is None:
                raise ValueError(
                    "Failed to initialize UNIQUAC activity model. Please check the components and model name.")

            # NOTE: check method
            tau_ij_cal_method = 0
            if dU_ij_src is None:
                # check if a_ij, b_ij, c_ij are provided
                if (a_ij_src is None or
                    b_ij_src is None or
                    c_ij_src is None or
                        d_ij_src is None):
                    raise ValueError(
                        "No valid source provided for interaction energy parameter (Δg_ij) or constants a, b, c, and d.")
                # set method
                tau_ij_cal_method = 2

                # ! a_ij
                if isinstance(a_ij_src, TableMatrixData):
                    a_ij = a_ij_src.mat('a', components)
                elif isinstance(a_ij_src, list):
                    a_ij = np.array(a_ij_src)
                elif isinstance(a_ij_src, np.ndarray):
                    a_ij = a_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (a_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! b_ij
                if isinstance(b_ij_src, TableMatrixData):
                    b_ij = b_ij_src.mat('b', components)
                elif isinstance(b_ij_src, list):
                    b_ij = np.array(b_ij_src)
                elif isinstance(b_ij_src, np.ndarray):
                    b_ij = b_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (b_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! c_ij
                if isinstance(c_ij_src, TableMatrixData):
                    c_ij = c_ij_src.mat('c', components)
                elif isinstance(c_ij_src, list):
                    c_ij = np.array(c_ij_src)
                elif isinstance(c_ij_src, np.ndarray):
                    c_ij = c_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (c_ij). Must be TableMatrixData, list of lists, or numpy array.")

                # ! d_ij
                if isinstance(d_ij_src, TableMatrixData):
                    d_ij = d_ij_src.mat('d', components)
                elif isinstance(d_ij_src, list):
                    d_ij = np.array(d_ij_src)
                elif isinstance(d_ij_src, np.ndarray):
                    d_ij = d_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (d_ij). Must be TableMatrixData, list of lists, or numpy array.")
            elif dU_ij_src is not None:
                # use dU_ij
                if isinstance(dU_ij_src, TableMatrixData):
                    dU_ij = dU_ij_src.mat('dU', components)
                elif isinstance(dU_ij_src, list):
                    dU_ij = np.array(dU_ij_src)
                elif isinstance(dU_ij_src, np.ndarray):
                    dU_ij = dU_ij_src
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (Δg_ij). Must be TableMatrixData, list of lists, or numpy array.")
                # set method
                tau_ij_cal_method = 1
            else:
                raise ValueError(
                    "No valid source provided for interaction energy parameter (Δg_ij) or constants A, B, C.")

            # SECTION: calculate tau_ij
            # NOTE: calculate the binary interaction parameter matrix (tau_ij)
            if tau_ij_cal_method == 1:
                # check
                if isinstance(dU_ij, np.ndarray):
                    tau_ij, _ = activity_uniquac.cal_tau_ij_M1(
                        temperature=temperature,
                        dU_ij=dU_ij)
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (Δg_ij). Must be numpy array.")
            elif tau_ij_cal_method == 2:
                # check
                if (isinstance(a_ij, np.ndarray) and
                    isinstance(b_ij, np.ndarray) and
                    isinstance(c_ij, np.ndarray) and
                        isinstance(d_ij, np.ndarray)):
                    # calculate tau_ij
                    tau_ij, _ = activity_uniquac.cal_tau_ij_M2(
                        temperature=temperature,
                        a_ij=a_ij,
                        b_ij=b_ij,
                        c_ij=c_ij,
                        d_ij=d_ij)
                else:
                    raise ValueError(
                        "Invalid source for interaction energy parameter (a_ij, b_ij, c_ij, d_ij). Must be numpy array.")
            else:
                raise ValueError(
                    "Invalid method for calculating tau_ij. Must be 1 or 2.")

            # NOTE: nrtl inputs
            inputs_ = {
                'mole_fraction': z_i_comp,
                "tau_ij": tau_ij,
                "r_i": r_i,
                "q_i": q_i
            }

            # NOTE: calculate activity
            res_, _ = activity_uniquac.cal(model_input=inputs_)

            # res format
            # res = {
            #     'property_name': 'activity coefficients',
            #     'components': components,
            #     'mole_fraction': xi,
            #     'value': AcCo_i,
            #     'unit': 1,
            #     'symbol': "AcCo_i",
            #     'message': message,
            # }

            # res
            return res_
        except Exception as e:
            raise Exception(
                f"Failed to calculate UNIQUAC activity: {e}") from e

NRTL(components, z_i_comp, temperature, **kwargs)

NRTL activity model for calculating activity coefficients.

Parameters

components : list List of component names. z_i_comp : dict Dictionary of component names and their respective mole fractions. temperature : float Temperature in Kelvin. kwargs : dict Additional parameters for the model. - interaction-energy-parameter : list, optional Interaction energy parameters for the components.

Source code in pyThermoFlash/docs/activity.py

def NRTL(self,
         components: List[str],
         z_i_comp: Dict[str, float],
         temperature: float,
         **kwargs):
    '''
    NRTL activity model for calculating activity coefficients.

    Parameters
    ----------
    components : list
        List of component names.
    z_i_comp : dict
        Dictionary of component names and their respective mole fractions.
    temperature : float
        Temperature in Kelvin.
    kwargs : dict
        Additional parameters for the model.
        - interaction-energy-parameter : list, optional
            Interaction energy parameters for the components.
    '''
    try:
        # SECTION: check src
        # extract activity model inputs
        activity_inputs = kwargs.get('activity_inputs', None)
        nrtl_datasource = kwargs.get('NRTL', None)

        # check
        if activity_inputs is None:
            # ! check nrtl inputs in datasource
            activity_inputs = nrtl_datasource

        # check if activity_inputs is None
        if activity_inputs is None:
            raise ValueError(
                "No valid source provided for activity model (NRTL) inputs.")

        # SECTION: check if activity_inputs is a dictionary
        if activity_inputs is not None:
            # check if activity_inputs is a dictionary
            if not isinstance(activity_inputs, dict):
                raise ValueError(
                    "activity_inputs must be a dictionary.")
            # check if activity_inputs is empty
            if len(activity_inputs) == 0:
                raise ValueError(
                    "activity_inputs cannot be empty.")

            # NOTE: update activity_inputs with nrtl_datasource
            if nrtl_datasource is not None:
                activity_inputs.update(nrtl_datasource)

        # SECTION: extract activity model inputs
        # NOTE: method 1
        # ! Δg_ij, interaction energy parameter
        dg_ij_src = activity_inputs.get('dg_ij', None)
        if dg_ij_src is None:
            dg_ij_src = activity_inputs.get('dg', None)

        # NOTE: method 2
        # ! constants a, b, c, and d
        a_ij_src = activity_inputs.get('a_ij', None)
        if a_ij_src is None:
            a_ij_src = activity_inputs.get('a', None)
        b_ij_src = activity_inputs.get('b_ij', None)
        if b_ij_src is None:
            b_ij_src = activity_inputs.get('b', None)
        c_ij_src = activity_inputs.get('c_ij', None)
        if c_ij_src is None:
            c_ij_src = activity_inputs.get('c', None)
        d_ij_src = activity_inputs.get('d_ij', None)
        if d_ij_src is None:
            d_ij_src = activity_inputs.get('d', None)

        # NOTE: α_ij, non-randomness parameter
        alpha_ij_src = activity_inputs.get('alpha_ij', None)
        if alpha_ij_src is None:
            alpha_ij_src = activity_inputs.get('alpha', None)

        # SECTION: init NRTL model
        # activity model
        activity = ptm.activity(
            components=components, model_name='NRTL')
        # set
        activity_nrtl = activity.nrtl
        # check
        if activity_nrtl is None:
            raise ValueError(
                "Failed to initialize NRTL activity model. Please check the components and model name.")

        # NOTE: check method
        tau_ij_cal_method = 0
        if dg_ij_src is None:
            # check if a_ij, b_ij, c_ij are provided
            if a_ij_src is None or b_ij_src is None or c_ij_src is None or d_ij_src is None:
                raise ValueError(
                    "No valid source provided for interaction energy parameter (Δg_ij) or constants a, b, c, and d.")
            # set method
            tau_ij_cal_method = 2

            # ! a_ij
            if isinstance(a_ij_src, TableMatrixData):
                a_ij = a_ij_src.mat('a', components)
            elif isinstance(a_ij_src, list):
                a_ij = np.array(a_ij_src)
            elif isinstance(a_ij_src, np.ndarray):
                a_ij = a_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (a_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! b_ij
            if isinstance(b_ij_src, TableMatrixData):
                b_ij = b_ij_src.mat('b', components)
            elif isinstance(b_ij_src, list):
                b_ij = np.array(b_ij_src)
            elif isinstance(b_ij_src, np.ndarray):
                b_ij = b_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (b_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! c_ij
            if isinstance(c_ij_src, TableMatrixData):
                c_ij = c_ij_src.mat('c', components)
            elif isinstance(c_ij_src, list):
                c_ij = np.array(c_ij_src)
            elif isinstance(c_ij_src, np.ndarray):
                c_ij = c_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (c_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! d_ij
            if isinstance(d_ij_src, TableMatrixData):
                d_ij = d_ij_src.mat('d', components)
            elif isinstance(d_ij_src, list):
                d_ij = np.array(d_ij_src)
            elif isinstance(d_ij_src, np.ndarray):
                d_ij = d_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (d_ij). Must be TableMatrixData, list of lists, or numpy array.")
        elif dg_ij_src is not None:
            # use dg_ij
            if isinstance(dg_ij_src, TableMatrixData):
                dg_ij = dg_ij_src.mat('dg', components)
            elif isinstance(dg_ij_src, list):
                dg_ij = np.array(dg_ij_src)
            elif isinstance(dg_ij_src, np.ndarray):
                dg_ij = dg_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (Δg_ij). Must be TableMatrixData, list of lists, or numpy array.")
            # set method
            tau_ij_cal_method = 1
        else:
            raise ValueError(
                "No valid source provided for interaction energy parameter (Δg_ij) or constants A, B, C.")

        # SECTION: extract data
        # α_ij, non-randomness parameter
        if isinstance(alpha_ij_src, TableMatrixData):
            alpha_ij = alpha_ij_src.mat('alpha', components)
        elif isinstance(alpha_ij_src, list):
            alpha_ij = np.array(alpha_ij_src)
        elif isinstance(alpha_ij_src, np.ndarray):
            alpha_ij = alpha_ij_src
        else:
            raise ValueError(
                "Invalid source for non-randomness parameter (α_ij). Must be TableMatrixData, list of lists, or numpy array.")

        # NOTE: calculate the binary interaction parameter matrix (tau_ij)
        if tau_ij_cal_method == 1:
            # check
            tau_ij, _ = activity_nrtl.cal_tau_ij_M1(
                temperature=temperature,
                dg_ij=dg_ij
            )
        elif tau_ij_cal_method == 2:
            tau_ij, _ = activity_nrtl.cal_tau_ij_M2(
                temperature=temperature,
                a_ij=a_ij,
                b_ij=b_ij,
                c_ij=c_ij,
                d_ij=d_ij
            )
        else:
            raise ValueError(
                "Invalid method for calculating tau_ij. Must be 1 or 2.")

        # NOTE: nrtl inputs
        inputs_ = {
            'mole_fraction': z_i_comp,
            "tau_ij": tau_ij,
            "alpha_ij": alpha_ij
        }

        # NOTE: calculate activity
        res_, _ = activity_nrtl.cal(model_input=inputs_)

        # res format
        # res = {
        #     'property_name': 'activity coefficients',
        #     'components': components,
        #     'mole_fraction': xi,
        #     'value': AcCo_i,
        #     'unit': 1,
        #     'symbol': "AcCo_i",
        #     'message': message,
        # }

        # res
        return res_
    except Exception as e:
        raise Exception(f"Failed to calculate NRTL activity: {e}") from e

UNIQUAC(components, z_i_comp, temperature, **kwargs)

UNIQUAC activity model for calculating activity coefficients.

Parameters

components : list List of component names. z_i_comp : dict Dictionary of component names and their respective mole fractions. temperature : float Temperature in Kelvin. kwargs : dict Additional parameters for the model. - interaction-energy-parameter : list, optional Interaction energy parameters for the components.

Source code in pyThermoFlash/docs/activity.py

def UNIQUAC(self,
            components: List[str],
            z_i_comp: Dict[str, float],
            temperature: float,
            **kwargs):
    '''
    UNIQUAC activity model for calculating activity coefficients.

    Parameters
    ----------
    components : list
        List of component names.
    z_i_comp : dict
        Dictionary of component names and their respective mole fractions.
    temperature : float
        Temperature in Kelvin.
    kwargs : dict
        Additional parameters for the model.
        - interaction-energy-parameter : list, optional
            Interaction energy parameters for the components.
    '''
    try:
        # SECTION: check src
        # extract activity model inputs
        activity_inputs = kwargs.get('activity_inputs', None)
        uniquac_datasource = kwargs.get('UNIQUAC', None)

        # check
        if activity_inputs is None:
            # ! check nrtl inputs in datasource
            activity_inputs = uniquac_datasource

        # check if activity_inputs is None
        if activity_inputs is None:
            raise ValueError(
                "No valid source provided for activity model (UNIQUAC) inputs.")

        # NOTE: check if activity_inputs is a dictionary
        if activity_inputs is not None:
            # check if activity_inputs is a dictionary
            if not isinstance(activity_inputs, dict):
                raise ValueError(
                    "activity_inputs must be a dictionary.")
            # check if activity_inputs is empty
            if len(activity_inputs) == 0:
                raise ValueError(
                    "activity_inputs cannot be empty.")

        # NOTE: update activity_inputs with uniquac_datasource
        if uniquac_datasource is not None:
            activity_inputs.update(uniquac_datasource)

        # NOTE: method 1
        # Δg_ij, interaction energy parameter
        dU_ij_src = activity_inputs.get('dU_ij', None)
        if dU_ij_src is None:
            dU_ij_src = activity_inputs.get('dU', None)

        # NOTE: method 2
        # constants a, b, c, and d
        a_ij_src = activity_inputs.get('a_ij', None)
        if a_ij_src is None:
            a_ij_src = activity_inputs.get('a', None)
        b_ij_src = activity_inputs.get('b_ij', None)
        if b_ij_src is None:
            b_ij_src = activity_inputs.get('b', None)
        c_ij_src = activity_inputs.get('c_ij', None)
        if c_ij_src is None:
            c_ij_src = activity_inputs.get('c', None)
        d_ij_src = activity_inputs.get('d_ij', None)
        if d_ij_src is None:
            d_ij_src = activity_inputs.get('d', None)

        # NOTE: r_i, relative van der Waals volume of component i
        r_i_src = activity_inputs.get('r_i', None)
        if r_i_src is None:
            r_i_src = activity_inputs.get('r', None)

        # final check if r_i is provided
        if r_i_src is None:
            raise ValueError("No valid source provided for r_i.")

        # check if r_i is a list or numpy array
        if isinstance(r_i_src, list):
            r_i = np.array(r_i_src)
        elif isinstance(r_i_src, np.ndarray):
            r_i = r_i_src
        else:
            raise ValueError(
                "Invalid source for r_i. Must be a list or numpy array.")

        # NOTE: q_i, relative van der Waals area of component i
        q_i_src = activity_inputs.get('q_i', None)
        if q_i_src is None:
            q_i_src = activity_inputs.get('q', None)

        # final check if q_i is provided
        if q_i_src is None:
            raise ValueError("No valid source provided for q_i.")

        # check if q_i is a list or numpy array
        if isinstance(q_i_src, list):
            q_i = np.array(q_i_src)
        elif isinstance(q_i_src, np.ndarray):
            q_i = q_i_src
        else:
            raise ValueError(
                "Invalid source for q_i. Must be a list or numpy array.")

        # SECTION: init NRTL model
        # activity model
        activity = ptm.activity(
            components=components, model_name='NRTL')
        # set
        activity_uniquac = activity.uniquac

        # check
        if activity_uniquac is None:
            raise ValueError(
                "Failed to initialize UNIQUAC activity model. Please check the components and model name.")

        # NOTE: check method
        tau_ij_cal_method = 0
        if dU_ij_src is None:
            # check if a_ij, b_ij, c_ij are provided
            if (a_ij_src is None or
                b_ij_src is None or
                c_ij_src is None or
                    d_ij_src is None):
                raise ValueError(
                    "No valid source provided for interaction energy parameter (Δg_ij) or constants a, b, c, and d.")
            # set method
            tau_ij_cal_method = 2

            # ! a_ij
            if isinstance(a_ij_src, TableMatrixData):
                a_ij = a_ij_src.mat('a', components)
            elif isinstance(a_ij_src, list):
                a_ij = np.array(a_ij_src)
            elif isinstance(a_ij_src, np.ndarray):
                a_ij = a_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (a_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! b_ij
            if isinstance(b_ij_src, TableMatrixData):
                b_ij = b_ij_src.mat('b', components)
            elif isinstance(b_ij_src, list):
                b_ij = np.array(b_ij_src)
            elif isinstance(b_ij_src, np.ndarray):
                b_ij = b_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (b_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! c_ij
            if isinstance(c_ij_src, TableMatrixData):
                c_ij = c_ij_src.mat('c', components)
            elif isinstance(c_ij_src, list):
                c_ij = np.array(c_ij_src)
            elif isinstance(c_ij_src, np.ndarray):
                c_ij = c_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (c_ij). Must be TableMatrixData, list of lists, or numpy array.")

            # ! d_ij
            if isinstance(d_ij_src, TableMatrixData):
                d_ij = d_ij_src.mat('d', components)
            elif isinstance(d_ij_src, list):
                d_ij = np.array(d_ij_src)
            elif isinstance(d_ij_src, np.ndarray):
                d_ij = d_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (d_ij). Must be TableMatrixData, list of lists, or numpy array.")
        elif dU_ij_src is not None:
            # use dU_ij
            if isinstance(dU_ij_src, TableMatrixData):
                dU_ij = dU_ij_src.mat('dU', components)
            elif isinstance(dU_ij_src, list):
                dU_ij = np.array(dU_ij_src)
            elif isinstance(dU_ij_src, np.ndarray):
                dU_ij = dU_ij_src
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (Δg_ij). Must be TableMatrixData, list of lists, or numpy array.")
            # set method
            tau_ij_cal_method = 1
        else:
            raise ValueError(
                "No valid source provided for interaction energy parameter (Δg_ij) or constants A, B, C.")

        # SECTION: calculate tau_ij
        # NOTE: calculate the binary interaction parameter matrix (tau_ij)
        if tau_ij_cal_method == 1:
            # check
            if isinstance(dU_ij, np.ndarray):
                tau_ij, _ = activity_uniquac.cal_tau_ij_M1(
                    temperature=temperature,
                    dU_ij=dU_ij)
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (Δg_ij). Must be numpy array.")
        elif tau_ij_cal_method == 2:
            # check
            if (isinstance(a_ij, np.ndarray) and
                isinstance(b_ij, np.ndarray) and
                isinstance(c_ij, np.ndarray) and
                    isinstance(d_ij, np.ndarray)):
                # calculate tau_ij
                tau_ij, _ = activity_uniquac.cal_tau_ij_M2(
                    temperature=temperature,
                    a_ij=a_ij,
                    b_ij=b_ij,
                    c_ij=c_ij,
                    d_ij=d_ij)
            else:
                raise ValueError(
                    "Invalid source for interaction energy parameter (a_ij, b_ij, c_ij, d_ij). Must be numpy array.")
        else:
            raise ValueError(
                "Invalid method for calculating tau_ij. Must be 1 or 2.")

        # NOTE: nrtl inputs
        inputs_ = {
            'mole_fraction': z_i_comp,
            "tau_ij": tau_ij,
            "r_i": r_i,
            "q_i": q_i
        }

        # NOTE: calculate activity
        res_, _ = activity_uniquac.cal(model_input=inputs_)

        # res format
        # res = {
        #     'property_name': 'activity coefficients',
        #     'components': components,
        #     'mole_fraction': xi,
        #     'value': AcCo_i,
        #     'unit': 1,
        #     'symbol': "AcCo_i",
        #     'message': message,
        # }

        # res
        return res_
    except Exception as e:
        raise Exception(
            f"Failed to calculate UNIQUAC activity: {e}") from e

__init__(datasource=None, equationsource=None)

Activity model class for pyThermoFlash.

Parameters

datasource : dict, optional Data source for the model, containing data for components and equations. equationsource : dict, optional Equation source for the model, containing equations for components.

Source code in pyThermoFlash/docs/activity.py

def __init__(self,
             datasource: Optional[Dict[str, Any]] = None,
             equationsource: Optional[Dict[str, Any]] = None
             ):
    '''
    Activity model class for pyThermoFlash.

    Parameters
    ----------
    datasource : dict, optional
        Data source for the model, containing data for components and equations.
    equationsource : dict, optional
        Equation source for the model, containing equations for components.
    '''
    # set datasource
    self.datasource = datasource
    # set equationsource
    self.equationsource = equationsource