Absorption and stripping: acetone scrubber, TCE air stripper, and the L/G optimum by gradient

Absorption and stripping: acetone scrubber, TCE air stripper, and the L/G optimum by gradient#

difflow.units.absorption provides a generic physical-solvent Absorber and Stripper (Kremser, stage-by-stage or packed HTU/NTU). This notebook (1) scrubs acetone from air with water, staged and packed; (2) designs an air stripper for TCE in groundwater; (3) finds the solvent rate that minimises a cost by jax.grad through an absorber-stripper loop with recycle.

Data caveat. The Henry’s constants in HENRY_WATER are transcribed from memory (order of magnitude, +/-20-30 %), and no published textbook example is reproduced: this notebook checks the model against its own closed forms, not against a book.

import jax, jax.numpy as jnp, numpy as np
jax.config.update("jax_enable_x64", True)
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from difflow import Absorber, Stripper, AbsorptionParams, Flowsheet, make_stream
from difflow.flowsheet import Unit
from difflow.mass_transfer import diffusion_volume, fuller_diffusivity, wilke_chang_diffusivity
from difflow.units.absorption import HENRY_WATER, henry_params, minimum_solvent_rate
from difflow.staged import linear_equilibrium, kremser_fraction

T0, P0 = 298.15, 101325.0

1. Acetone scrubbed from air with water#

1 mol/s of air carrying 2 mol% acetone meets 2 mol/s of water. Kremser, the rigorous stage-by-stage mode and a closed-form check of the dilute limit.

henry, henry_C = henry_params(["acetone"])
gas = make_stream({"air": 1.0, "acetone": 0.02}, T0, P0)
water = make_stream({"H2O": 2.0}, T0, P0)
K = HENRY_WATER["acetone"].H298 / P0
print(f"K(acetone, 25 C) = {K:.2f}  (from-memory table)")
for model in ("kremser", "stages"):
    ab = Absorber(AbsorptionParams(henry=henry, henry_C=henry_C, model=model, n_stages=6))
    treated, rich, info = ab(gas, water)
    print(f"{model:8s} recovery {float(info['recovery']['acetone']):.4f}   "
          f"L/L_min {float(info['L_over_Lmin']):.3f}   dT_adiabatic {float(info['delta_T_adiabatic']):.2f} K")
# dilute limit against the closed form
dil = make_stream({"air": 1.0, "acetone": 1e-9}, T0, P0)
_, _, i = Absorber(AbsorptionParams(henry=henry, henry_C=henry_C, n_stages=6, variable_flow=False))(dil, water)
A = 2.0 / K
print("dilute recovery", float(i["recovery"]["acetone"]), " closed form", float((A**7 - A) / (A**7 - 1)))
K(acetone, 25 C) = 1.97  (from-memory table)
kremser  recovery 0.8598   L/L_min 1.167   dT_adiabatic 4.71 K
stages   recovery 0.8602   L/L_min 1.166   dT_adiabatic 4.71 K
dilute recovery 0.8627217799271575  closed form 0.862721780346219

The solute-free operating line and the equilibrium curve (the unit returns both for plotting), and the stage count versus solvent rate, which blows up at the minimum rate.

_, _, info = Absorber(AbsorptionParams(henry=henry, henry_C=henry_C, n_stages=6))(gas, water)
fig, ax = plt.subplots(1, 2, figsize=(9, 3.6))
ax[0].plot(info["X"], info["Y"], label="operating line")
ax[0].plot(info["X"], info["Y_eq"], "--", label="equilibrium")
ax[0].set_xlabel("X (mol acetone / mol water)"); ax[0].set_ylabel("Y (mol acetone / mol air)"); ax[0].legend()
Lmin = 0.95 * K * 1.0   # dilute: A_min = recovery
Ls = Lmin * np.array([1.02, 1.05, 1.1, 1.2, 1.5, 2.0, 3.0])
Ns = []
for L in Ls:
    _, _, d = Absorber(AbsorptionParams(henry=henry, henry_C=henry_C, recovery=0.95, variable_flow=False))(
        make_stream({"air": 1.0, "acetone": 1e-9}, T0, P0), make_stream({"H2O": float(L)}, T0, P0))
    Ns.append(float(d["n_stages"]))
ax[1].plot(Ls / Lmin, Ns, "o-"); ax[1].set_xlabel("L / L_min"); ax[1].set_ylabel("stages for 95 %")
plt.tight_layout(); plt.savefig("/tmp/_abs.png", dpi=60)
print([f"{r:.2f}: {n:.1f}" for r, n in zip(Ls / Lmin, Ns)])
['1.02: 29.7', '1.05: 19.5', '1.10: 13.6', '1.20: 9.2', '1.50: 5.4', '2.00: 3.6', '3.00: 2.5']

Packed column with Onda films, rating and design (Pall rings 50 mm, 0.8 m diameter).

D_G = float(fuller_diffusivity(T0, P0, 58.08, 28.96, diffusion_volume({"C": 3, "H": 6, "O": 1}), 19.7))
D_L = float(wilke_chang_diffusivity(T0, 8.9e-4, 18.015, 74.0, 2.6))
gas_p = make_stream({"air": 20.0, "acetone": 0.4}, T0, P0)
water_p = make_stream({"H2O": 50.0}, T0, P0)
col = dict(henry=henry, henry_C=henry_C, model="packed", packing="pall_metal_50", diameter=0.8,
           D_G={"acetone": D_G}, D_L={"acetone": D_L})
_, _, d = Absorber(AbsorptionParams(recovery=0.95, **col))(gas_p, water_p)
Z = float(d["height"])
print(f"design: Z = {Z:.2f} m = H_OG {float(d['H_OG']['acetone']):.3f} m x N_OG {float(d['N_OG']['acetone']):.2f}")
_, _, r = Absorber(AbsorptionParams(height=Z, **col))(gas_p, water_p)
print("rating at that height: recovery", float(r["recovery"]["acetone"]))
_, _, c = Absorber(AbsorptionParams(recovery=0.95, ntu_method="integral", **col))(gas_p, water_p)
print(f"Colburn height {Z:.4f} m, numerical N_OG integral height {float(c['height']):.4f} m")
design: Z = 9.49 m = H_OG 1.226 m x N_OG 7.74
rating at that height: recovery 0.95
Colburn height 9.4937 m, numerical N_OG integral height 9.4165 m

2. Air stripping of TCE from groundwater#

A stripper with stripping factor S = K V / L. With clean air, the fraction left is (S-1)/(S^(N+1)-1). TCE has a large Henry’s constant (from memory), so a small air/water ratio works.

h_t, C_t = henry_params(["TCE"])
K_t = HENRY_WATER["TCE"].H298 / P0
L = 1000.0
for S in (1.5, 3.0, 5.0):
    air = make_stream({"air": S * L / K_t}, T0, P0)
    gw = make_stream({"H2O": L, "TCE": 1e-3}, T0, P0)
    _, _, i = Stripper(AbsorptionParams(henry=h_t, henry_C=C_t, recovery=0.99, variable_flow=False))(air, gw)
    print(f"S = {S}: air/water = {S/K_t:.4f} mol/mol, stages for 99 % = {float(i['n_stages']):.2f}, "
          f"G/G_min = {float(i['G_over_Gmin']):.2f}")
S = 1.5: air/water = 0.0028 mol/mol, stages for 99 % = 8.70, G/G_min = 1.52
S = 3.0: air/water = 0.0055 mol/mol, stages for 99 % = 3.83, G/G_min = 3.03
S = 5.0: air/water = 0.0092 mol/mol, stages for 99 % = 2.72, G/G_min = 5.05

3. The L/G optimum by gradient#

An absorber (acetone from the air) and a stripper (hot, fixed air flow) share a recycled water loop; a flow controller holds the water rate at L. The treated-gas acetone loss is differentiated by jax.grad through Flowsheet.solve (the recycle converges, then is differentiated implicitly).

class SolventMakeup:
    """Flow-controlled makeup/purge: rescales the recycled liquid to a water rate L."""
    def __init__(self, L): self.L = L
    def __call__(self, inlet):
        k = self.L / inlet["F_H2O"]
        return make_stream({"H2O": self.L, "acetone": inlet["F_acetone"] * k}, inlet["T"], inlet["P"])

def loop_loss(L):
    fs = Flowsheet(species_order=["acetone", "H2O"])
    fs.add_feed("gas", make_stream({"air": 1.0, "acetone": 0.02}, T0, P0))
    fs.add_feed("strip_air", make_stream({"air2": 0.6}, T0, P0))
    fs.add_unit(Unit("absorber", Absorber(AbsorptionParams(henry=henry, henry_C=henry_C, n_stages=6.0, T=T0)),
                     ["gas", "lean"], ["treated", "rich"]))
    fs.add_unit(Unit("stripper", Stripper(AbsorptionParams(henry=henry, henry_C=henry_C, n_stages=5.0, T=340.0)),
                     ["strip_air", "rich"], ["offgas", "lean_hot"]))
    fs.add_unit(Unit("makeup", SolventMakeup(L), ["lean_hot"], ["lean_set"]))
    fs.add_recycle("lean_set", "lean")
    s = fs.solve(tear_initial={"lean": make_stream({"acetone": 0.0, "H2O": L}, T0, P0)}, tol=1e-12, max_iter=200)
    return s["treated"]["F_acetone"]

Ls = np.linspace(0.8, 6.0, 7)
loss = np.array([float(loop_loss(L)) for L in Ls])
price_loss, price_solvent = 40.0, 0.15          # arbitrary weights, illustration only
cost = price_loss * loss + price_solvent * Ls
k = int(np.argmin(cost)); L0 = float(Ls[k])
g = float(jax.grad(loop_loss)(L0))
fd = (float(loop_loss(L0 + 1e-4)) - float(loop_loss(L0 - 1e-4))) / 2e-4
print(f"best grid L = {L0:.2f} mol/s;  d(loss)/dL = {g:.4e} (finite difference {fd:.4e})")
print(f"stationarity of the cost: price_loss*g + price_solvent = {price_loss*g + price_solvent:.4f}")
best grid L = 2.53 mol/s;  d(loss)/dL = -1.8820e-03 (finite difference -1.8820e-03)
stationarity of the cost: price_loss*g + price_solvent = 0.0747

Because the cost is differentiable in L, the optimum is the root of the stationarity condition price_loss * d(loss)/dL + price_solvent = 0, found by bisection on jax.grad (no scan needed).

dloss = jax.grad(loop_loss)
lo, hi = 0.8, 6.0
for _ in range(8):
    mid = 0.5 * (lo + hi)
    if price_loss * float(dloss(mid)) + price_solvent < 0.0:   # still improving as L grows
        lo = mid
    else:
        hi = mid
L = 0.5 * (lo + hi)
print(f"gradient-based optimum L = {L:.3f} mol/s (grid best {L0:.2f}); loss {float(loop_loss(L)):.5f} mol/s acetone")
gradient-based optimum L = 2.090 mol/s (grid best 2.53); loss 0.00245 mol/s acetone