stages-thermo

A staged-separation (distillation) learning library and fast steady-state column solver, built on vle-thermo.

import stages — walk the full pedagogical ladder of column methods (McCabe–Thiele → Ponchon–Savarit → Fenske–Underwood–Gilliland shortcut → rigorous MESH: bubble-point, sum-rates, inside-out, full Newton), each implemented from scratch and anchored to its textbook equations, with a granular, batch-capable API ("numpy for distillation columns"). All thermodynamics come from vle-thermo; this package adds none of its own.

The ladder's terminal target is an atmospheric crude distillation unit — hundreds of pseudocomponents, pumparounds, steam-stripped side strippers, and products specified on D86 95 % points and gaps rather than mole-fraction purity.

pip install stages-thermo           # import name is `stages`
pip install "stages-thermo[plot]"   # + matplotlib for the staircase diagrams

Status: 1.0 ships the first six rungs of the ladder and the batch layer over them — inside-out, the flagship rigorous solver (Milestone 11): inside_out solves any simple chain under any set of Specs (reflux ratio, product rates, boilup, duties, purities, recoveries, stage temperatures, the general Spec.stream_property, fixed and free side draws) by the Boston–Britt / Russell method — rigorous thermodynamics only in the outer loop, which fits a per-stage Kb/α/linear-enthalpy surrogate (fit_surrogate, Surrogate, StageSurrogate — itself a thermo provider), and Newton on the stripping factors inside; it lands on the same fixed points as Wang–Henke and sum-rates (≤ 7 × 10⁻⁶ K) and on Example 10.8 39× faster than Wang–Henke on the same SRK column, and its nonsmooth=True option solves columns whose stages run dry — and the two teaching solvers below it: the Wang–Henke bubble-point method (Milestone 6): wang_henke solves multi-feed columns with side draws under the classic reflux-ratio + distillate-rate specs on any thermo provider (Peng–Robinson, the SRK route, NRTL/van-Laar γ-φ, or the ideal teaching mocks), with trace mode for the convergence machinery, validated against Seader–Henley–Roper 4th-ed Examples 10.1 (exact) and 10.8 (0.04 K agreement with the book's CHEMCAD/SRK run); and the Burningham–Otto sum-rates method (Milestone 7): sum_rates solves absorbers and strippers (zero-degree-of-freedom columns) with the temperatures from Newton on all the energy balances at once (analytic Jacobian), seeded by seed_absorber, validated against Example 10.3 (the book's Aspen Plus SRK run — lean-gas temperature to 0.01 °F, component flows to ≤ 0.94 kmol/h). Below them sit the binary McCabe–Thiele (Milestone 1) and Ponchon–Savarit (Milestone 2) layers and the multicomponent FUG shortcut (Milestone 3): equilibrium and enthalpy–composition (H–x–y) curves from real thermodynamics, minimum reflux by geometric pinch detection (tangent pinches included), stage stepping with Murphree efficiency, total reflux, N(R), the energy-exact difference-point construction (with condenser/reboiler duties), the NRTL γ-φ model, per-phase enthalpies, Fenske–Underwood–Gilliland + Kirkbride + Winn with relative volatilities computed from vle-thermo K-values (fug), the multicomponent K-value / dew-point / flash surface on ThermoSystem, and the diagram plots (staircase, H–x–y, Gilliland). Milestone 5 added the MESH infrastructure every rigorous solver shares: Column (the multicomponent column as a stage graph — multi-feed incl. vapor-only, side draws, duties, pressure profiles), Profiles (the state/result object), mesh_residuals (through a ThermoSystem or the textbook IdealProvider), Spec incl. the general Spec.stream_property, seed_profiles / seed_from_fug, the thomas / thomas_component_balance kernels, SolveReport and ColumnSolution. Milestone 10 added the batch layer, solve_batch: one base column under many cases of overrides (spec targets, feed flows, pressures), solved by inside-out on every core with the GIL released, returned as numpy arrays with a failed case as a NaN row. 1.0 is inside-out as the flagship alongside Wang–Henke and sum-rates, plus the batch layer. The crude-column capability follows at 2.0, and Naphtali–Sandholm (the specialist for strongly nonideal columns) at 2.1 — see the repo's ROADMAP.md. From 1.0 the public API follows semantic versioning.

import stages

# Real thermodynamics from vle-thermo: Peng–Robinson benzene–toluene at 1 atm
# (light component first; units K / kPa absolute / mole fractions).
sys = stages.ThermoSystem.peng_robinson(["benzene", "toluene"])
curve = stages.EquilibriumCurve.from_thermo(sys, 101.325)

# Minimum reflux by pinch detection — tangent pinches included …
r = stages.rmin(curve, x_distillate=0.95, x_bottoms=0.05, z_feed=0.50, q=1.0)

# … then the full construction. Rich result objects, never bare numbers:
design = stages.mccabe_thiele(curve, 0.95, 0.05, 0.50, reflux=1.5 * r.r_min)
print(f"N = {design.n_stages:.2f} stages, feed stage {design.feed_stage}, "
      f"R_min = {r.r_min:.3f}")

design.stages          # every (x, y) stage corner
design.staircase       # the full polyline, ready to plot
design.rectifying      # operating lines as slope/intercept
r.pinch, r.tangent     # where the pinch sits, and whether it's a tangent pinch

# The staircase diagram (requires the [plot] extra):
# from stages import plotting
# plotting.plot_mccabe_thiele(design, curve, show_rmin=True)

Strongly non-ideal systems go through the γ-φ route — the same construction, different thermodynamics:

# van Laar methanol–water (the system validated in vle's Chapter IV notebooks).
mw = stages.ThermoSystem.van_laar(["methanol", "water"], 0.5853, 0.3458)
curve_mw = stages.EquilibriumCurve.from_thermo(mw, 101.325)

Rung 2 — Ponchon–Savarit — closes the energy balance on the enthalpy–composition (H–x–y) diagram, so it also returns the condenser and reboiler duties (which McCabe–Thiele cannot):

# H–x–y curve: saturated-liquid and -vapor enthalpies alongside y*(x).
ec = stages.EnthalpyCurve.from_thermo(sys, 101.325)
ps = stages.ponchon_savarit(ec, x_distillate=0.95, x_bottoms=0.05,
                            z_feed=0.50, reflux=1.5)
print(f"N = {ps.n_stages:.2f} stages, feed stage {ps.feed_stage}")
print(f"Q_C/F = {ps.q_condenser:,.0f}, Q_R/F = {ps.q_reboiler:,.0f} kJ/kmol feed")
ps.delta_d, ps.delta_b   # the two difference points (poles), (x, H) in kJ/kmol

# NRTL for strongly non-ideal aqueous-organic systems (a12/a21 in kJ/kmol):
aw = stages.ThermoSystem.nrtl(["ammonia", "water"], a12=-1800.0, a21=-1200.0, alpha=0.2)
# from stages import plotting; plotting.plot_ponchon_savarit(ps, ec)

Rung 3 — the FUG shortcut — is multicomponent. Every relative volatility comes from vle-thermo K-values at the two ends of the column (the "FUG(K)" form that a simulator uses to seed its rigorous solver), refined to self-consistency; each rung is also callable on its own (fenske_n_min, underwood_roots, gilliland_stages, kirkbride_split, winn_fit, …):

# A depropanizer: propane (LK) / isobutane (HK), 98 % key recoveries, R = 1.3 R_min.
c3 = stages.ThermoSystem.peng_robinson(["ethane", "propane", "isobutane", "n-butane", "n-pentane"])
fug = stages.fug(c3, 1500.0, feed=[5, 40, 20, 25, 10], light_key=1, heavy_key=2,
                 lk_recovery=0.98, hk_recovery=0.98, q=1.0, reflux_factor=1.3)
print(fug)   # FugResult(N_min=12.95, R_min=1.885, R=2.450, N=26.47, feed_stage=15, ...)
fug.underwood_roots, fug.n_min_winn, fug.d, fug.x_d, fug.t_top   # every intermediate
# from stages import plotting; plotting.plot_gilliland(fug)

Milestone 5's infrastructure is callable on its own — the multicomponent column as a stage graph, the shortcut-seeded starting state, and the MESH residuals of any state through the provider boundary:

import math

# The depropanizer above as a column: total condenser at stage 0, partial
# reboiler last, the feed on the FUG feed stage (0-based indices).
n = math.ceil(fug.n_stages) + 1
col = stages.Column.simple(n, 5, "total", "partial", 1500.0)
col = col.with_feed(fug.feed_stage, [5, 40, 20, 25, 10], "saturated_liquid")
col.degrees_of_freedom          # 2 — e.g. Spec.reflux_ratio + Spec.product_rate

# Seed: linear T, constant-molal-overflow flows, x_D → x_B blend, Kb, α, S.
prof = stages.seed_from_fug(col, c3, fug)
prof.t, prof.l, prof.v, prof.s   # flat per-stage arrays; prof.x_stage(j) per stage

# The MESH residuals of that state (K, H_L, H_V are written back into prof).
res = stages.mesh_residuals(col, c3, prof)
print(res)   # MeshResiduals(max|M|=..., max|E|=..., max|S|=..., max|H|=...)

# Specifications reduce to residual rows — including the general variant.
specs = [stages.Spec.reflux_ratio(fug.reflux),
         stages.Spec.stream_property("distillate", lambda s: s["flows"][1], 39.2, name="propane in D")]
stages.validate_specs(col, specs)
stages.spec_residuals(col, prof, specs)

And Milestone 6 drives those residuals to zero — the Wang–Henke bubble-point solver, here on the five-stage column of Seader–Henley–Roper 4th-ed Examples 10.1/10.8 (SRK, 100 psia = 689.476 kPa, R = 2, D = 50 kmol/h):

srk = stages.ThermoSystem.soave_redlich_kwong(["propane", "n-butane", "n-pentane"])
col = stages.Column.simple(5, 3, "total", "partial", 689.476)
col = col.with_feed(2, [30.0, 30.0, 40.0], "saturated_liquid")
seed = stages.seed_profiles(col, srk, t_top=291.5, t_bottom=347.0,
                            reflux_ratio=2.0, distillate_rate=50.0,
                            x_top=[0.55, 0.40, 0.05], x_bottom=[0.02, 0.25, 0.73])
sol = stages.wang_henke(col, srk, 2.0, 50.0, seed, tol_sum_dt2=1e-8, trace=True)
sol.report.converged                 # True (18 iterations)
sol.report.outer.residual_norms      # τ = Σ(ΔT)² per iteration [K²]
sol.report.final_residual            # rigorous scaled MESH ∞-norm at the answer
sol.profiles.t                       # converged stage temperatures [K]
sol.condenser_duty, sol.reboiler_duty
stages.product_stream(col, sol.profiles, "bottoms")["flows"]
# → [0.955, 12.357, 36.687] kmol/h — the book's CHEMCAD run prints
#   [0.955, 12.363, 36.683]; stage temperatures agree to 0.04 K.

wang_henke_composition_step exposes the tridiagonal step alone — the exact kernel the inside-out inner loop iterates on.

Milestone 7 flips the tearing for wide-boiling systems — the sum-rates solver on the six-stage absorber of S&H 4th-ed Example 10.3 (SRK; n-decane absorbent at 90 °F / 400 psia on top, a C₁–C₅ gas at 105 °F at the bottom; no condenser, no reboiler, nothing to specify):

srk = stages.ThermoSystem.soave_redlich_kwong(
    ["methane", "ethane", "propane", "n-butane", "n-pentane", "n-decane"])
psia, f_to_k = 6.894757, lambda f: (f - 32.0) / 1.8 + 273.15
col = stages.Column.simple(6, 6, None, None, [(400.0 + 0.6 * j) * psia for j in range(6)])
col = col.with_feed(0, [0, 0, 0, 0.05, 0.78, 164.17], "temperature", t=f_to_k(90.0))   # absorbent
col = col.with_feed(5, [160, 370, 240, 25, 5, 0], "temperature", t=f_to_k(105.0))       # feed gas
col.degrees_of_freedom               # 0 — an absorber takes no specs
seed = stages.seed_absorber(col, srk, f_to_k(90.0), f_to_k(105.0))   # oil down, gas up, linear T
sol = stages.sum_rates(col, srk, seed, tol_sum_dt2=1e-8, trace=True)
sol.report.converged, sol.report.outer.iterations, sol.report.inner.iterations  # True, 19, 44
stages.product_stream(col, sol.profiles, "vapor:0")["flows"]   # lean gas, kmol/h
# → [146.20, 269.56, 99.45, 1.39, 0.24, 0.84]; Aspen Plus (book): [146.53, 270.50, 99.96, 1.41, 0.24, 0.83]
(sol.profiles.t[0] - 273.15) * 1.8 + 32   # lean-gas temperature → 151.40 °F (book: 151.4)

IdealProvider.demo_absorber() is the wide-boiling teaching mock (gas / solute / oil) the notebook's exercises run on.

Milestone 11 is the flagship — inside-out on the same Figure-10.6 column, now under any specification set, with the two-level machinery reported:

srk = stages.ThermoSystem.soave_redlich_kwong(["propane", "n-butane", "n-pentane"])
col = stages.Column.simple(5, 3, "total", "partial", 689.476)
col = col.with_feed(2, [30.0, 30.0, 40.0], "saturated_liquid")
seed = stages.seed_profiles(col, srk, 291.5, 347.0, 2.0, 50.0, [0.55, 0.40, 0.05], [0.02, 0.25, 0.73])
specs = [stages.Spec.reflux_ratio(2.0), stages.Spec.purity("distillate", 0, 0.60)]  # any Spec kind
sol = stages.inside_out(col, srk, specs, seed, trace=True)
sol.report.converged, sol.report.outer.iterations, sol.report.inner.iterations   # True, 8, 43
sol.report.outer.residual_norms[-1]  # the RIGOROUS scaled MESH ∞-norm at the answer (≈ 3e-8)
sol.report.surrogate_drift[-1]       # how far the last refit moved the local models (≈ 2e-8)
stages.product_stream(col, sol.profiles, "distillate")["rate"]   # → 48.12 kmol/h for x_D(C3) = 0.60
sur = stages.fit_surrogate(col, srk, sol.profiles)   # the per-stage models the outer loop fits
sur.stage(2).b, sur.stage(2).alpha                    # the Kb slope [K] and the frozen α's
sur.stage(2).bubble_point([0.2, 0.5, 0.3])[0]         # a StageSurrogate is a thermo provider

Milestone 10 is the batch layer — "numpy for distillation columns": one call solves a whole sweep by inside-out on every core (GIL released), one numpy row per case, and a failed case is a NaN row rather than an exception:

import numpy as np

seed = {"t_top": 291.5, "t_bottom": 347.0, "reflux_ratio": 2.0, "distillate_rate": 50.0,
        "x_top": [0.55, 0.40, 0.05], "x_bottom": [0.02, 0.25, 0.73]}
specs = [stages.Spec.reflux_ratio(2.0), stages.Spec.product_rate("distillate", 50.0)]
r = np.linspace(1.5, 4.0, 12)
out = stages.solve_batch(col, srk, specs, seed, spec_values=np.column_stack([r, np.full_like(r, 50.0)]))
out["converged"].all()        # True
out["x_distillate"][:, 0]     # propane in the distillate per R: 0.578 … 0.586
out["reboiler_duty"]          # kJ/h per R: 2.66e6 … 5.08e6
# Also: feed_flows=(n_cases, n_comp), pressure=(n_cases,) or (n_cases, n_stages),
# n_threads=…, return_solutions=True for every ColumnSolution.

Measured in a 4-vCPU cloud container: a 40-stage, 10-component inside-out solve takes 2.25 ms on the ideal mock and 24.8 ms on rigorous SRK; a 64-case sweep of the ideal column runs 3.7× faster through solve_batch than in a loop (see the repo README's Performance section).

The executable learning path lives in the repo's notebooks/ — 01-mccabe-thiele.ipynb, 02-ponchon-savarit.ipynb and 03-shortcut-design.ipynb design benzene–toluene, methanol–water, acetone–water and ammonia–water binaries (the last on a reference chart digitized from the Pátek–Klomfar correlation) and a C2–C5 depropanizer end-to-end, 04-mesh-and-bubble-point.ipynb opens the rigorous machinery — the tridiagonal structure, the convergence history, the MESH residuals — 05-absorbers-sum-rates.ipynb flips it for absorbers (Example 10.3 live, the Friday–Smith pairing rule measured, Kremser vs rigorous), and 06-inside-out.ipynb opens the flagship (the Kb fit, drift vs rigorous residual, the same answer as Wang–Henke and the book's CHEMCAD run, a stripper that runs dry), and 07-numpy-for-columns.ipynb puts solve_batch to work (a 336-column rigorous N-vs-R map against Gilliland, feed-sensitivity heatmaps, failure as data), with exercises throughout.

The native core is a Rust crate (stages-thermo on crates.io) with PyO3 bindings; wheels are abi3 (cp310-abi3-*), so one wheel per (OS, arch) covers CPython 3.10+.

License

MIT © Miguel Roberto Jackson Ugueto