polykin.thermo.flash

flash2_PV

flash2_PV(
    F: float,
    z: FloatVector,
    P: float,
    beta: float,
    Kcalc: Callable[
        [float, float, FloatVector, FloatVector],
        FloatVector,
    ],
    *,
    T0: float = 300.0,
    maxiter: int = 50,
    atol_inner: float = 1e-09,
    rtol_outer: float = 1e-06
) -> FlashResult

Solve a 2-phase flash problem at given pressure and vapor fraction.

References

• J.F. Boston, H.I. Britt, "A radically different formulation and solution of the single-stage flash problem", Computers & Chemical Engineering, Volume 2, Issues 2-3, 1978, p. 109-122,

PARAMETER	DESCRIPTION
F	Feed mole flowrate (mol/s). TYPE: float
z	Feed mole fractions (mol/mol). TYPE: FloatVector
P	Pressure (Pa). TYPE: float
beta	Vapor phase fraction (mol/mol). TYPE: float
Kcalc	Function to calculate K-values, with signature Kcalc(T, P, x, y). TYPE: Callable[[float, float, FloatVector, FloatVector], FloatVector]
T0	Initial guess for temperature (K). TYPE: float DEFAULT: 300.0
maxiter	Maximum number of iterations. TYPE: int DEFAULT: 50
atol_inner	Absolute tolerance for the inner R-loop. TYPE: float DEFAULT: 1e-09
rtol_outer	Relative tolerance for the outer volatility parameters loop. TYPE: float DEFAULT: 1e-06

RETURNS	DESCRIPTION
FlashResult	Flash result.

Source code in src/polykin/thermo/flash/vle.py

def flash2_PV(
    F: float,
    z: FloatVector,
    P: float,
    beta: float,
    Kcalc: Callable[[float, float, FloatVector, FloatVector], FloatVector],
    *,
    T0: float = 300.0,
    maxiter: int = 50,
    atol_inner: float = 1e-9,
    rtol_outer: float = 1e-6,
) -> FlashResult:
    r"""Solve a 2-phase flash problem at given pressure and vapor fraction.

    **References**

    *  J.F. Boston, H.I. Britt, "A radically different formulation and solution
       of the single-stage flash problem", Computers & Chemical Engineering,
       Volume 2, Issues 2-3, 1978, p. 109-122,

    Parameters
    ----------
    F : float
        Feed mole flowrate (mol/s).
    z : FloatVector
        Feed mole fractions (mol/mol).
    P : float
        Pressure (Pa).
    beta : float
        Vapor phase fraction (mol/mol).
    Kcalc : Callable[[float, float, FloatVector, FloatVector], FloatVector]
        Function to calculate K-values, with signature `Kcalc(T, P, x, y)`.
    T0 : float
        Initial guess for temperature (K).
    maxiter : int
        Maximum number of iterations.   
    atol_inner : float
        Absolute tolerance for the inner R-loop.
    rtol_outer : float
        Relative tolerance for the outer volatility parameters loop.

    Returns
    -------
    FlashResult
        Flash result.
    """
    # Initial guesses
    z = z/z.sum()
    T = T0
    K = Kcalc(T, P, z, z)
    x = z/(1 + beta*(K - 1))
    y = K*x
    x /= x.sum()
    y /= y.sum()

    # Inner loop objective function
    def fobj(R, u, Kb0, beta, z) -> float:
        p = z/(1 - R + Kb0*R*exp(u))
        return 1 - beta - (1 - R)*p.sum()

    # Outer loop
    Tref = 300.0
    u, Kb, A, B = _parameters_PV(T, P, x, y, beta, Kcalc, Tref, all=True)
    Kb0 = Kb
    success = False
    for _ in range(maxiter):

        v_old = np.concatenate((u, [A]))

        # Inner R-loop
        if abs(beta - 0) <= eps:
            R = 0.0
        elif abs(beta - 1) < eps:
            R = 1.0
        else:
            sol = fzero_brent(lambda R: fobj(R, u, Kb0, beta, z),
                              0.0, 1.0,
                              maxiter=maxiter,
                              xtol=atol_inner,
                              ftol=atol_inner)
            R = sol.x
            if not sol.success:
                warnings.warn(
                    f"Inner R-loop did not converge after {maxiter} iterations.",
                    ConvergenceWarning)

        # Compute x, y
        p = z/(1 - R + Kb0*R*exp(u))
        eup = exp(u)*p
        sum_p = p.sum()
        sum_eup = eup.sum()
        Kb = sum_p/sum_eup
        T = 1/(1/Tref + (log(Kb) - A)/B)
        x = p/sum_p
        y = eup/sum_eup

        # Update u, A
        u, A = _parameters_PV(T, P, x, y, beta, Kcalc, Tref, all=False, B=B)
        v_new = np.concatenate((u, [A]))

        # Check convergence
        v0 = min(k for k in v_new if k > 0)
        if np.allclose(v_new, v_old, atol=v0*rtol_outer):
            success = True
            break

    else:
        warnings.warn(
            f"Outer loop did not converge after {maxiter} iterations.",
            ConvergenceWarning)

    # Overall balances
    V = F*beta
    L = F - V

    return FlashResult(success, T, P, F, L, V, beta, z, x, y, K)