Distillation & diagrams¶

Shortcut (Fenske-Underwood-Gilliland) and rigorous equilibrium-stage distillation columns, plus binary diagrams, azeotropes, and residue-curve maps.

Distillation columns¶

column ¶

Distillation columns: shortcut (Fenske-Underwood-Gilliland) and rigorous stages.

Two complementary, fully differentiable column models:

Shortcut: the Fenske-Underwood-Gilliland (FUG) method, giving minimum stages at total reflux (fenske_min_stages), minimum reflux (underwood_min_reflux), the actual stage count at a working reflux (gilliland_stages), and the feed-stage location (kirkbride_feed_stage), tied together by shortcut_column. Cheap, robust, and ideal for screening or as an initial guess for the rigorous model.
Rigorous: a multistage equilibrium-stage column solved by the Wang-Henke bubble-point method under constant molar overflow, with EOS K-values on every stage (solve_column). The converged profile is differentiable through the fixed-point iteration by implicit differentiation.

Both expose exact gradients of their outputs (stage count, reflux, product purities) with respect to the design variables, so a column can be embedded in a gradient-based optimisation alongside the rest of a flowsheet.

Component ordering convention: relative volatilities alpha are given relative to a common reference (any component); the light key lk is more volatile than the heavy key hk (alpha[lk] > alpha[hk]).

Classes:

Name	Description
`ShortcutResult`	Summary of a Fenske-Underwood-Gilliland shortcut design.
`ColumnResult`	Converged profile and products of a rigorous equilibrium-stage column.

Functions:

Name	Description
`relative_volatility`	Relative volatilities `alpha_i = K_i / K_ref` at `(T, P, z)` from the EOS.
`fenske_min_stages`	Fenske minimum number of equilibrium stages at total reflux.
`underwood_min_reflux`	Underwood minimum reflux ratio `R_min` (constant relative volatility).
`gilliland_stages`	Actual equilibrium-stage count from the Gilliland (Eduljee) correlation.
`kirkbride_feed_stage`	Kirkbride correlation for the number of stages above the feed.
`shortcut_column`	Full FUG shortcut design from a feed and a specified product split.
`solve_column`	Rigorous multistage column by the Wang-Henke bubble-point method (CMO).

ShortcutResult ¶

Bases: NamedTuple

Summary of a Fenske-Underwood-Gilliland shortcut design.

Attributes:

Name	Type	Description
`n_min`	`Array`	Minimum equilibrium stages at total reflux (Fenske).
`r_min`	`Array`	Minimum reflux ratio (Underwood).
`theta`	`Array`	Underwood common root.
`r`	`Array`	Working reflux ratio used for the Gilliland step.
`n_stages`	`Array`	Actual equilibrium stages (Gilliland).
`feed_stage`	`Array`	Number of stages above the feed (Kirkbride).

ColumnResult ¶

Bases: NamedTuple

Converged profile and products of a rigorous equilibrium-stage column.

Attributes:

Name	Type	Description
`t`	`Array`	Stage temperatures (K), top stage first.
`x`	`Array`	Liquid mole fractions, shape `(n_stages, n_components)`.
`y`	`Array`	Vapour mole fractions, shape `(n_stages, n_components)`.
`distillate`	`Stream`	Distillate product `Stream`.
`bottoms`	`Stream`	Bottoms product `Stream`.
`reflux`	`Array`	Reflux ratio used.
`condenser_duty`	`Array`	Condenser heat removed (W, positive).
`reboiler_duty`	`Array`	Reboiler heat added (W, positive).

relative_volatility ¶

relative_volatility(
    eos: CubicEOS,
    t: ArrayLike,
    p: ArrayLike,
    z: Array,
    tc: Array,
    pc: Array,
    omega: Array,
    *,
    ref: int,
    kij: Array | None = None,
) -> Array

Relative volatilities alpha_i = K_i / K_ref at (T, P, z) from the EOS.

K-values are evaluated as phi_i^L(z) / phi_i^V(z) at the given composition, a standard shortcut estimate that is well defined whether or not the feed is two-phase. Differentiable in (T, P, z).

fenske_min_stages ¶

fenske_min_stages(
    d: Array, b: Array, lk: int, hk: int, alpha: Array
) -> Array

Fenske minimum number of equilibrium stages at total reflux.

N_min = ln[(d_LK/d_HK)(b_HK/b_LK)] / ln(alpha_LK/alpha_HK) where d and b are the distillate and bottoms component molar flows. Includes the reboiler as a stage (the classic Fenske count).

underwood_min_reflux ¶

underwood_min_reflux(
    z: Array,
    x_d: Array,
    alpha: Array,
    q: ArrayLike,
    lk: int,
    hk: int,
    *,
    tol: float = 1e-12,
    max_iter: int = 200,
) -> tuple[Array, Array]

Underwood minimum reflux ratio R_min (constant relative volatility).

Solves the first Underwood equation sum_i alpha_i z_i / (alpha_i - theta) = 1 - q for the common root theta between alpha_HK and alpha_LK, then evaluates R_min + 1 = sum_i alpha_i x_{i,D} / (alpha_i - theta).

Returns (R_min, theta); both are differentiable in (z, x_d, alpha, q).

gilliland_stages ¶

gilliland_stages(
    n_min: ArrayLike, r_min: ArrayLike, r: ArrayLike
) -> Array

Actual equilibrium-stage count from the Gilliland (Eduljee) correlation.

With X = (R - R_min)/(R + 1) and Y = 0.75 (1 - X^0.5668), the stage count follows from (N - N_min)/(N + 1) = Y, i.e. N = (N_min + Y)/(1 - Y).

kirkbride_feed_stage ¶

kirkbride_feed_stage(
    n: ArrayLike,
    z: Array,
    x_d: Array,
    x_b: Array,
    d_total: ArrayLike,
    b_total: ArrayLike,
    lk: int,
    hk: int,
) -> Array

Kirkbride correlation for the number of stages above the feed.

log10(N_R/N_S) = 0.206 log10[(z_HK/z_LK)(x_{LK,B}/x_{HK,D})^2 (B/D)]; returns N_R given the total stage count N = N_R + N_S.

shortcut_column ¶

shortcut_column(
    z: Array,
    d: Array,
    b: Array,
    alpha: Array,
    q: ArrayLike,
    lk: int,
    hk: int,
    *,
    reflux: ArrayLike | None = None,
    reflux_factor: ArrayLike = 1.3,
) -> ShortcutResult

Full FUG shortcut design from a feed and a specified product split.

Parameters:

Name	Type	Description	Default
`z`	`Array`	Feed mole fractions.	required
`d`	`Array`	Distillate component molar flows (the chosen split of the feed).	required
`b`	`Array`	Bottoms component molar flows (`z * F - d` on a consistent basis).	required
`alpha`	`Array`	Relative volatilities (to any common reference).	required
`q`	`ArrayLike`	Feed thermal quality (1 = saturated liquid, 0 = saturated vapour).	required
`lk`	`int`	Light-key component index.	required
`hk`	`int`	Heavy-key component index.	required
`reflux`	`ArrayLike \| None`	Working reflux ratio. If `None`, `reflux_factor * R_min` is used.	`None`
`reflux_factor`	`ArrayLike`	Multiplier on `R_min` when `reflux` is not given.	`1.3`

Returns:

Type	Description
`ShortcutResult`	A `ShortcutResult`; every field is differentiable in the inputs.

solve_column ¶

solve_column(
    feed: Stream,
    n_stages: int,
    feed_stage: int,
    reflux: ArrayLike,
    distillate_rate: ArrayLike,
    *,
    eos: CubicEOS = PR,
    q: ArrayLike = 1.0,
    kij: Array | None = None,
    t_top: ArrayLike | None = None,
    t_bottom: ArrayLike | None = None,
    t_min: float = 100.0,
    t_max: float = 800.0,
    tol: float = 1e-09,
    max_iter: int = 400,
) -> ColumnResult

Rigorous multistage column by the Wang-Henke bubble-point method (CMO).

A total condenser sits above stage 1 (distillate and reflux share the stage-1 vapour composition) and a partial reboiler is stage n_stages; a single feed of quality q enters at feed_stage (1-indexed from the top). Under constant molar overflow the section flows are fixed by reflux and distillate_rate, and each outer sweep (i) solves the tridiagonal component balances for the liquid profile with the current EOS K-values, (ii) takes a bubble-point Newton step on every stage temperature, and (iii) refreshes the vapour compositions. The sweep is iterated to a fixed point by tear_solve, so the converged profile, products, and duties are all differentiable (implicit differentiation) with respect to reflux, distillate_rate, the feed, and model parameters.

Parameters:

Name	Type	Description	Default
`feed`	`Stream`	Feed `Stream` (its temperature sets the feed enthalpy used for the duty balance).	required
`n_stages`	`int`	Number of equilibrium stages including the reboiler.	required
`feed_stage`	`int`	1-indexed feed stage (`2 <= feed_stage <= n_stages - 1`).	required
`reflux`	`ArrayLike`	Reflux ratio `L/D`.	required
`distillate_rate`	`ArrayLike`	Distillate molar flow (mol/s); bottoms is the remainder.	required
`eos`	`CubicEOS`	Cubic equation of state (default Peng-Robinson).	`PR`
`q`	`ArrayLike`	Feed thermal quality (1 = saturated liquid).	`1.0`
`kij`	`Array \| None`	Optional binary interaction matrix.	`None`
`t_top`	`ArrayLike \| None`	Optional initial top-stage temperature for the linear starting profile (default `feed.t` minus a small spread).	`None`
`t_bottom`	`ArrayLike \| None`	Optional initial bottom-stage temperature for the linear starting profile (default `feed.t` plus a small spread).	`None`
`t_min`	`float`	Lower bracket clamp for the per-stage temperature updates.	`100.0`
`t_max`	`float`	Upper bracket clamp for the per-stage temperature updates.	`800.0`
`tol`	`float`	Convergence tolerance for the outer fixed point.	`1e-09`
`max_iter`	`int`	Maximum number of outer sweeps.	`400`

Returns:

Type	Description
`ColumnResult`	A `ColumnResult`.

Binary diagrams & residue curves¶

diagrams ¶

Binary phase-diagram data and azeotrope finding for an equilibrium model.

Conceptual design and column screening lean on the binary picture: the P-x-y and T-x-y envelopes and, above all, where the azeotrope is (it bounds what ordinary distillation can reach). These helpers turn any binary EquilibriumModel into:

pxy_diagram / txy_diagram: the bubble (liquid) and the equilibrium-vapour curves on a composition grid, from one bubble sweep each; and
azeotrope_pressure / azeotrope_temperature: the azeotropic composition (where y_1 = x_1) at fixed T or P, returned with an exists flag so a non-azeotropic system is reported, not faked.

All outputs are differentiable: the diagram arrays through the bubble solves, and the azeotrope locus through the bracketed root's implicit derivative, so an azeotrope's pressure sensitivity to a model parameter is a single jax.grad away.

Classes:

Name	Description
`PxyDiagram`	Isothermal P-x-y data for a binary.
`TxyDiagram`	Isobaric T-x-y data for a binary.
`AzeotropeResult`	A binary azeotrope locus (or the best bracketed point if none exists).
`ResidueCurve`	A simple-distillation residue curve: the still-liquid composition path.

Functions:

Name	Description
`pxy_diagram`	Isothermal P-x-y curves for a binary on an `n`-point composition grid.
`txy_diagram`	Isobaric T-x-y curves for a binary on an `n`-point composition grid.
`azeotrope_pressure`	Find the binary azeotrope at fixed temperature `t` (where `y_1 = x_1`).
`azeotrope_temperature`	Find the binary azeotrope at fixed pressure `p` (where `y_1 = x_1`).
`residue_curve`	Integrate one residue curve of open (Rayleigh) distillation at fixed `P`.
`residue_curve_map`	Trace a residue-curve map: one full curve through each starting composition.

PxyDiagram ¶

Bases: NamedTuple

Isothermal P-x-y data for a binary.

Attributes:

Name	Type	Description
`x1`	`Array`	Liquid mole fraction of component 1 (the grid).
`y1`	`Array`	Equilibrium vapour mole fraction of component 1 (the bubble vapour).
`p`	`Array`	Bubble pressure at each `x1` (Pa).
`t`	`Array`	The (fixed) temperature (K).

TxyDiagram ¶

Bases: NamedTuple

Isobaric T-x-y data for a binary.

Attributes:

Name	Type	Description
`x1`	`Array`	Liquid mole fraction of component 1 (the grid).
`y1`	`Array`	Equilibrium vapour mole fraction of component 1.
`t`	`Array`	Bubble temperature at each `x1` (K).
`p`	`Array`	The (fixed) pressure (Pa).

AzeotropeResult ¶

Bases: NamedTuple

A binary azeotrope locus (or the best bracketed point if none exists).

Attributes:

Name	Type	Description
`exists`	`Array`	`True` if `y_1 - x_1` changes sign on the search bracket (a genuine azeotrope); `False` means the returned point is just a bracket end and should be ignored.
`x1`	`Array`	Azeotropic composition (`y_1 = x_1` there).
`t`	`Array`	Temperature (K).
`p`	`Array`	Pressure (Pa).

ResidueCurve ¶

Bases: NamedTuple

A simple-distillation residue curve: the still-liquid composition path.

Attributes:

Name	Type	Description
`x`	`Array`	Liquid composition along the curve, shape `(steps + 1, n)` (each row a mole-fraction vector on the composition simplex).
`t`	`Array`	Bubble (boiling) temperature at each point (K), shape `(steps + 1,)`.
`p`	`Array`	The (fixed) pressure of the map (Pa).

pxy_diagram ¶

pxy_diagram(
    model: _BinaryModel,
    t: ArrayLike,
    *,
    n: int = 51,
    eps: float = 0.001,
) -> PxyDiagram

Isothermal P-x-y curves for a binary on an n-point composition grid.

The grid is clipped to [eps, 1 - eps] to avoid the pure-component limits. Both curves come from a single bubble-pressure sweep: plot (x1, p) for the liquid line and (y1, p) for the vapour line.

txy_diagram ¶

txy_diagram(
    model: _BinaryModel,
    p: ArrayLike,
    *,
    n: int = 51,
    eps: float = 0.001,
    t_min: float = 200.0,
    t_max: float = 600.0,
) -> TxyDiagram

Isobaric T-x-y curves for a binary on an n-point composition grid.

Each grid point solves a bubble temperature (bracketed in [t_min, t_max]) and returns the equilibrium vapour, giving both the liquid line (x1, t) and the vapour line (y1, t) from one sweep. Widen the bracket for very light or very heavy pairs.

azeotrope_pressure ¶

azeotrope_pressure(
    model: _BinaryModel,
    t: ArrayLike,
    *,
    x_lo: float = 0.001,
    x_hi: float = 1.0 - 0.001,
    tol: float = 1e-10,
    max_iter: int = 200,
) -> AzeotropeResult

Find the binary azeotrope at fixed temperature t (where y_1 = x_1).

Brackets the root of y_1(x_1) - x_1 from the bubble-pressure relation. The returned x1 and p are differentiable with respect to the model parameters; check exists before trusting them.

azeotrope_temperature ¶

azeotrope_temperature(
    model: _BinaryModel,
    p: ArrayLike,
    *,
    x_lo: float = 0.001,
    x_hi: float = 1.0 - 0.001,
    t_min: float = 200.0,
    t_max: float = 600.0,
    tol: float = 1e-10,
    max_iter: int = 200,
) -> AzeotropeResult

Find the binary azeotrope at fixed pressure p (where y_1 = x_1).

As azeotrope_pressure but at constant pressure: each residual evaluation solves a bubble temperature (bracketed in [t_min, t_max]).

residue_curve ¶

residue_curve(
    model: _BubbleTemperatureModel,
    x0: Array,
    p: ArrayLike,
    *,
    steps: int = 200,
    ds: float = 0.05,
    direction: float = 1.0,
    t_min: float = 150.0,
    t_max: float = 700.0,
) -> ResidueCurve

Integrate one residue curve of open (Rayleigh) distillation at fixed P.

Solves the residue-curve ODE dx_i/dtau = x_i - y_i(x), where y(x) is the bubble-point vapour in equilibrium with the still liquid x. With direction = +1 the curve advances toward the high-boiling stable node (the still-pot residue enriches in heavy components); direction = -1 traces it back toward the low-boiling node.

Integration is an explicit Euler march with the composition re-projected onto the simplex each step, so the path stays a valid composition. The expensive part, the bubble-point (T, y) at each point, is one compiled, differentiable solve (_bubble_ty) reused across every step and every curve.

Parameters:

Name	Type	Description	Default
`model`	`_BubbleTemperatureModel`	Any object with a `bubble_temperature(p, x)` method (e.g. an `EquilibriumModel` from `fugacio.sim.models`).	required
`x0`	`Array`	Starting liquid composition, shape `(n,)`.	required
`p`	`ArrayLike`	Pressure (Pa).	required
`steps`	`int`	Number of integration steps.	`200`
`ds`	`float`	Pseudo-time step size.	`0.05`
`direction`	`float`	`+1` toward heavies, `-1` toward lights.	`1.0`
`t_min`	`float`	Lower bound of the bubble-temperature search bracket (K).	`150.0`
`t_max`	`float`	Upper bound of the bubble-temperature search bracket (K).	`700.0`

Returns:

Type	Description
`ResidueCurve`	A `ResidueCurve` with the composition and temperature path
`ResidueCurve`	(`steps + 1` points).

residue_curve_map ¶

residue_curve_map(
    model: _BubbleTemperatureModel,
    starts: Array,
    p: ArrayLike,
    *,
    steps: int = 150,
    ds: float = 0.05,
    t_min: float = 150.0,
    t_max: float = 700.0,
) -> list[ResidueCurve]

Trace a residue-curve map: one full curve through each starting composition.

Each start is integrated both directions (toward the light and the heavy node) and the two halves are stitched into a single curve passing through the start. This is the ternary distillation designer's master diagram: distillation boundaries and reachable products fall out of the family of curves.

Parameters:

Name	Type	Description	Default
`model`	`_BubbleTemperatureModel`	Object exposing `bubble_temperature(p, x)`.	required
`starts`	`Array`	Starting compositions, shape `(k, n)`.	required
`p`	`ArrayLike`	Pressure (Pa).	required
`steps`	`int`	Steps taken in each direction (the curve has `2*steps + 1` points).	`150`
`ds`	`float`	Pseudo-time step size.	`0.05`
`t_min`	`float`	Lower bound of the bubble-temperature search bracket (K).	`150.0`
`t_max`	`float`	Upper bound of the bubble-temperature search bracket (K).	`700.0`

Returns:

Type	Description
`list[ResidueCurve]`	A list of `k` `ResidueCurve` objects.