Two-Phase Flash¶

Problem statement¶

An isothermal-isobaric two-phase flash specifies \((T,P,\boldsymbol{z})\) and solves for liquid composition \(\boldsymbol{x}\), vapor composition \(\boldsymbol{y}\), and vapor fraction \(\beta\). The component balances are

\[ z_i=(1-\beta)x_i+\beta y_i. \]

Defining equilibrium ratios \(K_i=y_i/x_i\) gives

\[ x_i=\frac{z_i}{1+\beta(K_i-1)},\qquad y_i=K_i x_i. \]

Normalization of both phases produces the Rachford-Rice equation

\[ F(\beta) = \sum_i \frac{z_i(K_i-1)} {1+\beta(K_i-1)} =0. \]

The material-balance equation originates with Rachford and Rice (1952). Its admissible root must keep all denominators and phase compositions positive.

Coupling material balance and thermodynamics¶

For a \(\phi\)-\(\phi\) equation-of-state calculation,

\[ f_i^L=x_i\phi_i^L P,\qquad f_i^V=y_i\phi_i^V P. \]

The equilibrium condition \(f_i^L=f_i^V\) implies

\[ K_i=\frac{\phi_i^L}{\phi_i^V}. \]

Because each fugacity coefficient depends on its phase composition, the Rachford-Rice balance and the fugacity equations must be iterated together. torch-flash uses:

optional tangent-plane stability analysis;
explicit initial \(K\) values or a Wilson estimate when critical constants are available;
a bounded Rachford-Rice solve;
damped successive substitution in \(\ln K\); and
an autodifferentiated Newton correction with line search.

This follows the phase-split strategy of Michelsen, Part II. The reported residual is

\[ \max_i\left| \ln K_i-\left(\ln\phi_i^L-\ln\phi_i^V\right) \right|, \]

which is the maximum dimensionless log-fugacity mismatch at the current split.

Runnable example¶

"""Solve a known two-phase methane/n-butane TP state."""

import torch

from torch_flash import (
    ChemicalState,
    component_set,
    configure,
    peng_robinson_1978,
    two_phase_flash,
)

runtime = configure(device="cpu", dtype=torch.float64)
model = peng_robinson_1978(component_set(("methane", "n_butane")))
state = ChemicalState(
    temperature=runtime.tensor(270.0),  # K
    pressure=runtime.tensor(3.0e6),  # Pa
    composition=runtime.tensor([0.50, 0.50]),  # overall z
)

# This example deliberately solves a state already known to be on a
# two-phase branch. Leave check_stability=True when the phase count is unknown.
result = two_phase_flash(
    model,
    state,
    check_stability=False,
    tolerance=1.0e-10,
    raise_on_failure=True,
)

liquid, vapor = result.phases
print("phase fractions [liquid, vapor] =", result.phase_fractions)
print("liquid x =", liquid.composition)
print("vapor y =", vapor.composition)
print("maximum log-fugacity residual =", float(result.residual_norm))

# Always check material balance and the numerical result before using it.
reconstructed_feed = (
    result.phase_fractions[0] * liquid.composition + result.phase_fractions[1] * vapor.composition
)
torch.testing.assert_close(reconstructed_feed, state.composition, atol=1.0e-9, rtol=0.0)
assert result.converged

The explicit check_stability=False makes this a solve on a known two-phase branch. For a general feed state, retain the default stability test:

result = two_phase_flash(
    model,
    state,
    check_stability=True,
    tolerance=1.0e-8,
    raise_on_failure=True,
)

If the stability result is homogeneous, result.nphases == 1. If a split is returned, phase fractions are ordered with the liquid-like solver root first and the vapor-like root second; result.phase_kinds contains the separate physical-identification diagnostic.

Failure modes and interpretation¶

No Rachford-Rice bracket: the current \(K\) values do not support two positive phases for the specified feed.
Trivial solution: \(K_i\rightarrow1\) and \(\boldsymbol{x}\rightarrow\boldsymbol{y}\). Do not interpret this algebraic state as coexistence without a dedicated critical calculation.
Near-critical ill-conditioning: phase compositions and roots coalesce, so small residuals need not imply well-resolved phase properties.
Metastable split: equality of fugacities can hold even though another phase or phase count has lower Gibbs energy.
Non-convergence: use raise_on_failure=True for workflows that must not continue with an invalid result; otherwise inspect the emitted convergence warning and result.converged.

For batched states already known to be two phase, use batched_two_phase_flash. It intentionally does not perform per-state phase discovery.

Trust-region alternative for difficult states¶

two_phase_trust_region_flash minimizes the two-phase Gibbs energy in direct component mole amounts. The initially larger phase amount for each component is dependent, so every accepted iterate satisfies material balance exactly. The dense trust-region subproblem uses the exact PyTorch-autodifferentiated gradient and Hessian. tangent_plane_stability(..., minimizer="trust-region") selects the corresponding formal-mole-number TPD minimizer.

from torch_flash import two_phase_trust_region_flash

result = two_phase_trust_region_flash(
    model,
    state,
    check_stability=True,
    tolerance=1.0e-8,
    raise_on_failure=True,
)

The formulation follows Petitfrere and Nichita's restricted-step phase equilibrium method (doi:10.1016/j.fluid.2013.08.039). It is an opt-in path, not the default replacement for the \(\ln K\) flash: exact Hessians cost more than the small linear systems in ordinary well-conditioned two-phase states.

For a fixed three-or-more-phase hypothesis, multiphase_trust_region_flash applies the paper's per-component reference phase and direct-mole Gibbs formulation. It returns explicit residual, material-balance, accepted/rejected-step, and Hessian-curvature diagnostics. The requested phase count and roots remain scientific inputs; convergence of this local minimum is not proof of global stability.