polykin.thermo.flash

solve_Rachford_Rice

solve_Rachford_Rice(
    K: FloatVector,
    z: FloatVector,
    beta0: float,
    *,
    maxiter: int = 50,
    atol_res: float = 1e-09
) -> RachfordRiceResult

Solve the Rachford-Rice flash residual equation.

The numerical solution of the Rachford-Rice equation is carried out using a combination of the Newton and bisection methods, ensuring efficient and robust convergence.

References

• M.L. Michelsen and J. Mollerup, "Thermodynamic models: Fundamentals and computational aspects", Tie-Line publications, 2nd edition, 2007.

PARAMETER	DESCRIPTION
K	K-values. TYPE: FloatVector(N)
z	Feed mole fractions [mol/mol]. TYPE: FloatVector(N)
beta0	Initial guess for vapor phase fraction [mol/mol]. If NaN, an initial guess is automatically computed. TYPE: float
maxiter	Maximum number of iterations. TYPE: int DEFAULT: 50
atol_res	Absolute tolerance for residual. TYPE: float DEFAULT: 1e-09

RETURNS	DESCRIPTION
RachfordRiceResult	Result.

See Also

• residual_Rachford_Rice: related method to determine the residual and its derivative.

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

def solve_Rachford_Rice(
    K: FloatVector,
    z: FloatVector,
    beta0: float,
    *,
    maxiter: int = 50,
    atol_res: float = 1e-9,
) -> RachfordRiceResult:
    r"""Solve the Rachford-Rice flash residual equation.

    The numerical solution of the Rachford-Rice equation is carried out using
    a combination of the Newton and bisection methods, ensuring efficient and
    robust convergence.

    **References**

    *  M.L. Michelsen and J. Mollerup, "Thermodynamic models: Fundamentals and
       computational aspects", Tie-Line publications, 2nd edition, 2007.

    Parameters
    ----------
    K : FloatVector(N)
        K-values.
    z : FloatVector(N)
        Feed mole fractions [mol/mol].
    beta0 : float
        Initial guess for vapor phase fraction [mol/mol]. If `NaN`, an initial
        guess is automatically computed.
    maxiter : int
        Maximum number of iterations.
    atol_res : float
        Absolute tolerance for residual.

    Returns
    -------
    RachfordRiceResult
        Result.

    See Also
    --------
    * [`residual_Rachford_Rice`](residual_Rachford_Rice.md): related method to
      determine the residual and its derivative.
    """
    # Trivial subcooled and superheated cases
    F0 = residual_Rachford_Rice(0.0, K, z)[0]
    if F0 < 0.0:
        return RachfordRiceResult(True, 0, 0.0, F0)
    F1 = residual_Rachford_Rice(1.0, K, z)[0]
    if F1 > 0.0:
        return RachfordRiceResult(True, 0, 1.0, F1)

    # Bounds on beta
    beta_min = np.where(K > 1.0, (K * z - 1) / (K - 1), 0.0)
    beta_min = beta_min[(beta_min > 0.0) & (beta_min < 1.0)]
    beta_min = max(beta_min) if len(beta_min) > 0 else 0.0
    beta_max = np.where(K < 1.0, (1 - z) / (1 - K), 1.0)
    beta_max = beta_max[(beta_max > 0.0) & (beta_max < 1.0)]
    beta_max = min(beta_max) if len(beta_max) > 0 else 1.0

    # Initial guess
    if not np.isnan(beta0):
        beta = np.clip(beta0, beta_min, beta_max)
    else:
        beta = (beta_min + beta_max) / 2

    # Iteration loop
    success = False
    for iter in range(maxiter):

        F, dF = residual_Rachford_Rice(beta, K, z, derivative=True)

        # Check convergence
        if abs(F) <= atol_res:
            success = True
            break

        # Update bounds
        if F > 0.0:
            beta_min = beta
        else:
            beta_max = beta

        # Update beta (Newton or bisection)
        beta_new = beta - F / dF
        if (beta_new > beta_min) and (beta_new < beta_max):
            beta = beta_new
        else:
            beta = (beta_min + beta_max) / 2

    return RachfordRiceResult(success, iter + 1, beta, F)