Batch distillation: Rayleigh from scratch, and a still schedule by gradient

difflow.units.batch_distillation solves simple (Rayleigh) distillation and a batch rectifier at constant reflux or constant distillate composition. This notebook (1) builds the Rayleigh equation from scratch and compares it with rayleigh, (2) checks a batch rectifier against an independent time integration, and (3) schedules a heads / hearts / tails run of an ethanol-water still and tunes the hearts reflux with jax.grad.

Constant-alpha numbers are illustrations, not published examples; no textbook worked example is reproduced. The ethanol-water NRTL and Antoine constants are those of tests/test_phase_diagrams.py (see its docstring for provenance).

import jax, jax.numpy as jnp, numpy as np
from scipy import integrate
jax.config.update("jax_enable_x64", True)
from difflow.staged import constant_alpha
from difflow.units.batch_distillation import (
    BatchDistillation, BatchDistillationParams, rayleigh, batch_cuts, equilibrium_curve)

1. Rayleigh from scratch

Differential balance on a still of W mol: d(W x) = y dW  =>  ln(F/W) = int_{x_W}^{x_F} dx/(y*-x). Constant alpha has the closed form; the library also integrates any tabulated curve.

alpha, xF, WF = 2.5, 0.5, 0.4
y = lambda x: alpha*x/(1+(alpha-1)*x)
# from scratch: find x_W with ln(F/W) = integral, using scipy
from scipy.optimize import brentq
lnFW = lambda xW: integrate.quad(lambda x: 1/(y(x)-x), xW, xF)[0]
xW_scratch = brentq(lambda xW: lnFW(xW) - np.log(1/WF), 1e-6, xF)
r = rayleigh(constant_alpha(alpha), xF, W_over_F=WF)
print("from scratch x_W =", round(xW_scratch, 8), "  library x_W =", round(float(r["x_W"]), 8))
print("average distillate x_D,avg =", round(float(r["x_D_avg"]), 5))

from scratch x_W = 0.29675619   library x_W = 0.29675619
average distillate x_D,avg = 0.6355

2. A rectifier against an independent time integration

Three stages above the still, reflux 3, boil-up 2 mol/s, 100 mol charge. The library integrates in x_W; here the batch ODEs are integrated in time with solve_ivp and a NumPy tray-to-tray march.

from scipy.optimize import brentq
def march(xD, L, N):
    yv = xD
    for _ in range(N):
        x = yv/(alpha-(alpha-1)*yv); yv = L*x+(1-L)*xD
    return yv/(alpha-(alpha-1)*yv)
def rhs(t, s, R=3.0, V=2.0, N=3):
    W, M = s; xW = M/W
    xd = brentq(lambda xd: march(xd, R/(R+1), N)-xW, xW, 1-1e-12)
    D = V/(R+1); return [-D, -D*xd]
ref = integrate.solve_ivp(rhs, [0, 30.0], [100.0, 50.0], rtol=1e-10, atol=1e-12)
W_ref, M_ref = ref.y[:, -1]

unit = BatchDistillation(BatchDistillationParams(species_order=["a", "b"], n_stages=3, R=3.0, boilup=2.0),
                         curve=constant_alpha(alpha))
charge = {"F_a": 50.0, "F_b": 50.0, "T": 350.0, "P": 1e5}      # moles, not mol/s
still, dist, info = unit(charge, t_final=30.0, n_traj=11)
print("W      :", float(info["W_final"]), W_ref)
print("x_W    :", float(info["x_W_final"]), M_ref/W_ref)
print("balance:", float(still["F_a"] + dist["F_a"]), "(charge 50.0)")

W      : 85.00000000000001 85.00000000000001
x_W    : 0.4229442087588032 0.42294420875324995
balance: 50.0 (charge 50.0)

3. Whisky-still schedule

Ethanol-water (NRTL) with 6 stages above the still, a 100 mol charge of 20 mol% ethanol. Heads: 4 mol at high reflux; hearts: the next 12 mol; tails: to the end. The azeotrope (about 0.89) caps what any reflux can do, and the unit says so in info.

from difflow.thermo import IdealThermo, SpeciesData
from difflow.units.lle import NRTLParams
sd = lambda n, a, Tc: SpeciesData(name=n, MW=1.0, Cp_coeffs=(0,0,0,0), Hvap_coeffs=(1.0, 0.38, Tc), antoine_coeffs=a)
thermo = IdealThermo({"ethanol": sd("ethanol", (10.24677, 1598.673, -46.424), 514.0),
                      "water": sd("water", (9.6543, 1435.264, -64.848), 647.0)})
nrtl = NRTLParams(species=("ethanol", "water"), a=jnp.array([[0., -0.8009], [3.4578, 0.]]),
                  b=jnp.array([[0., 246.18], [-586.0809, 0.]]), alpha=jnp.array([[0., .3], [.3, 0.]]))
p = BatchDistillationParams(species_order=["ethanol", "water"], n_stages=6, R=8.0, boilup=1.0)
still_unit = BatchDistillation(p, thermo, nrtl)
ch = {"F_ethanol": 20.0, "F_water": 80.0, "T": 350.0, "P": 101325.0}
res = batch_cuts(still_unit, ch, [{"distillate": 4.0},
                                  {"stop": {"distillate": 12.0}, "R": 4.0},
                                  {"stop": {"x_W": 0.02}, "R": 2.0}])
for name, c in zip(["heads", "hearts", "tails"], res["cuts"]):
    print(f"{name:7s} {float(c['moles']):6.2f} mol  x_EtOH = {float(c['x_key']):.3f}  t = {float(c['time']):7.1f} s  "
          f"azeotrope-limited: {bool(c['info']['azeotrope_limited'])}")
print("balance (should be ~0):", {k: float(v) for k, v in res["balance"].items()})
print("x_azeotrope:", float(res['cuts'][0]['info']['x_azeotrope']))

heads     4.00 mol  x_EtOH = 0.820  t =    36.0 s  azeotrope-limited: False
hearts   12.00 mol  x_EtOH = 0.805  t =    60.0 s  azeotrope-limited: False
tails     7.53 mol  x_EtOH = 0.735  t =    22.6 s  azeotrope-limited: False
balance (should be ~0): {'ethanol': -1.9984014443252818e-15, 'water': 0.0}
x_azeotrope: 0.893451014802547

Tuning the hearts reflux by gradient

Question: what is the smallest reflux (the shortest cut) that makes the hearts at least 0.82 mol fraction ethanol? Purity rises with reflux and so does time (t = D_total (R+1)/V). jax.grad runs through the whole cut sequence, stopping compositions included, so a Newton iteration on x_hearts(R) = 0.82 needs no finite differences.

def hearts(R_hearts):
    res = batch_cuts(still_unit, ch, [{"distillate": 4.0}, {"stop": {"distillate": 12.0}, "R": R_hearts}])
    c = res["cuts"][1]
    return c["x_key"], c["time"]
vg = jax.jit(jax.value_and_grad(lambda R: hearts(R)[0]))
R = 4.0
for k in range(6):
    x, dx = vg(R)
    R = R - (float(x) - 0.82)/float(dx)
    print(f"iter {k}: R = {R:8.4f}   (x_hearts before the step {float(x):.5f}, dx/dR = {float(dx):.5f})")
x, t = hearts(R)
print(f"R = {R:.3f} gives x_hearts = {float(x):.5f} in {float(t):.1f} s;  R = 4 gives {float(hearts(4.0)[0]):.5f} in {float(hearts(4.0)[1]):.1f} s")

iter 0: R =   6.4951   (x_hearts before the step 0.80463, dx/dR = 0.00616)
iter 1: R =   8.8038   (x_hearts before the step 0.81434, dx/dR = 0.00245)
iter 2: R =   9.8628   (x_hearts before the step 0.81856, dx/dR = 0.00136)
iter 3: R =  10.0011   (x_hearts before the step 0.81985, dx/dR = 0.00109)

iter 4: R =  10.0030   (x_hearts before the step 0.82000, dx/dR = 0.00106)
iter 5: R =  10.0030   (x_hearts before the step 0.82000, dx/dR = 0.00106)

R = 10.003 gives x_hearts = 0.82000 in 132.0 s;  R = 4 gives 0.80463 in 60.0 s