McCabe-Thiele stepper: fractional stages, pinches, optimum reflux

difflow.staged steps a McCabe-Thiele diagram with a fixed-length lax.scan and returns a fractional, differentiable stage count. This notebook (1) reproduces a hand calculation at constant alpha, (2) steps the real benzene-toluene curve and compares with the rigorous DistillationColumn, and (3) finds the optimum R/R_min by gradient.

The constant-alpha numbers are an illustration computed for the issue, not a published example. No textbook example is reproduced here.

import jax, jax.numpy as jnp, numpy as np
jax.config.update("jax_enable_x64", True)
from difflow import staged as S

alpha, z, xD, xB = 2.5, 0.4, 0.95, 0.05
curve = S.constant_alpha(alpha)
R_min = S.minimum_reflux(curve, z, 1.0, xD)
print("R_min =", round(float(R_min), 4), "  N_min (Fenske) =", round(float(S.total_reflux_stages(curve, xD, xB)), 4))

def column(curve, R, q=1.0, **kw):
    return S.mccabe_thiele(curve, S.rectifying_line(R, xD), S.stripping_line(R, xD, xB, z, q),
                           xD, xB, q_line=S.q_line(q, z), **kw)

res = column(curve, 1.5 * R_min)
print("N =", round(float(res["N"]), 3), " -> integer stages", int(np.ceil(float(res["N"]))),
      " feed stage", round(float(res["feed_stage"]), 2))

R_min = 1.4444   N_min (Fenske) = 6.4269

N = 11.722  -> integer stages 12  feed stage 6.0

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from difflow.visualization import plot_mccabe_thiele
R = 1.5 * R_min
ax = plot_mccabe_thiele(res, curve, top_line=S.rectifying_line(R, xD),
                        bottom_line=S.stripping_line(R, xD, xB, z, 1.0), q_line=S.q_line(1.0, z))
ax.figure.set_size_inches(5, 5)
plt.show()

The real benzene-toluene curve vs the rigorous column

xy_curve gives the Raoult’s-law curve (alpha varies 2.4-2.6 with x); tabulated builds a monotone C1 interpolant and its inverse.

from difflow.database import get_species_data
from difflow.thermo import IdealThermo
from difflow.phase_diagrams import xy_curve, bubble_T
from difflow import DistillationColumn, make_stream
from difflow.units.distillation import DistillationColumnParams
from difflow.streams import get_flows

sp = ["benzene", "toluene"]
th = IdealThermo({s: get_species_data(s) for s in sp})
d = xy_curve(th, sp, 101325.0, n=81)
bt = S.tabulated(d["x"], d["y"])
Rmin_bt = S.minimum_reflux(bt, z, 1.0, xD)
R = 1.5 * Rmin_bt
res_bt = column(bt, R)
print("R_min", round(float(Rmin_bt), 3), " N_min", round(float(S.total_reflux_stages(bt, xD, xB)), 2),
      " N", round(float(res_bt["N"]), 2), " feed stage", round(float(res_bt["feed_stage"]), 1))

Tb = bubble_T(th, sp, jnp.array([z, 1 - z]), 101325.0)
D = 100.0 * (z - xB) / (xD - xB)
for N in (11, 12, 13):
    params = DistillationColumnParams(species_order=sp, n_stages=N, feed_stage=6,
                                      condenser_type="total", P=101325.0)
    dist, bot, info = DistillationColumn(params, th)(
        make_stream({"benzene": 40.0, "toluene": 60.0}, T=float(Tb), P=101325.0), R=float(R), D_spec=D)
    fd, fb = get_flows(dist), get_flows(bot)
    print(N, "stages: x_D =", round(float(fd["benzene"] / (fd["benzene"] + fd["toluene"])), 4),
          " x_B =", round(float(fb["benzene"] / (fb["benzene"] + fb["toluene"])), 4))

R_min 1.48  N_min 6.63  N 11.89  feed stage 6.0

11 stages: x_D = 0.9311  x_B = 0.0621

12 stages: x_D = 0.9451  x_B = 0.0532

13 stages: x_D = 0.955  x_B = 0.0469

The stepper asks for about 12 stages with the feed on stage 6; the rigorous column (constant molar overflow, ideal K-values) reaches about 0.945 / 0.053 with 12 and meets the 0.95 / 0.05 specification with 13. They agree to about one stage; the stepper’s integer count is the usual graphical estimate (the specification is met exactly at the fractional N).

Optimum R/R_min by gradient

An integer stage count has no gradient; the fractional one does. A crude cost: stages x diameter (vapor load grows as R + 1) plus reboiler steam (also growing as R + 1), in arbitrary units. The optimum sits where the two trade off; compare the 1.2-1.5 rule. The weights are arbitrary, so the number is only a demonstration of the gradient, not a design result.

def cost(ratio, w_steam=1.0):
    r = ratio * Rmin_bt
    N = column(bt, r)["N"]
    vapor = (r + 1.0) * D
    return N * (vapor ** 0.5) + w_steam * vapor     # capital ~ stages x diameter, + steam

grad = jax.jit(jax.grad(cost))
lo, hi = 1.05, 3.0                       # bisection on the sign of the exact gradient
for _ in range(40):
    mid = 0.5 * (lo + hi)
    lo, hi = (mid, hi) if grad(mid) < 0 else (lo, mid)
ratio = 0.5 * (lo + hi)
print("optimum R/R_min ~", round(ratio, 3), " cost", round(float(cost(ratio)), 2))
print("grad there", float(grad(ratio)), "(a kink of the piecewise-smooth N: the gradient changes sign across it)")
for r_ in (1.05, 1.1, 1.2, 1.3, 1.5, 2.0):
    print(r_, round(float(cost(r_)), 2))

optimum R/R_min ~ 1.35  cost 256.94
grad there -19.82776432913315 (a kink of the piecewise-smooth N: the gradient changes sign across it)
1.05 301.28

1.1 279.81
1.2 263.57

1.3 258.02
1.5 258.22

2.0 278.7