Mass transfer: diffusivities, films, Onda packing and column height

Contents

Mass transfer: diffusivities, films, Onda packing and column height#

difflow.mass_transfer collects the building blocks of Ch 14-15: diffusivity correlations, film theory, Sherwood correlations, two-resistance coefficients and the packed-column Onda / HTU / NTU machinery. Everything is differentiable.

Honesty note: no textbook worked example is reproduced here. The checks are analytic oracles (stagnant film vs its closed form, Colburn NTU vs quadrature) and the one data point I am confident of (water vapour in air, about 2.5e-5 m^2/s).

import jax, jax.numpy as jnp, numpy as np
from scipy.integrate import quad
from difflow import mass_transfer as mt
from difflow.units.evaporator import water_saturation_pressure

1. Diffusivities and a stagnant-film evaporation#

P, T = 101325.0, 298.15
D = mt.fuller_diffusivity(T, P, 18.015, 28.96, mt.FULLER_VOLUMES["H2O"], mt.FULLER_VOLUMES["air"])
print(f"Fuller D(H2O-air) = {float(D):.3e} m^2/s  (data ~2.5e-5)")

c = P / (mt.R_GAS * T)
y_sat = water_saturation_pressure(T) / P
for delta in (0.01, 0.1):
    N = mt.stagnant_film_flux(D, c, y_sat, 0.0, delta)
    ref = D * c / delta * np.log(1.0 / (1.0 - float(y_sat)))      # closed form
    print(f"film {delta*100:4.0f} cm: N = {float(N):.4e} mol/m2/s  closed form {ref:.4e}  "
          f"equimolar {float(mt.equimolar_flux(D, c, y_sat, 0.0, delta)):.4e}")

Fuller D(H2O-air) = 2.509e-05 m^2/s  (data ~2.5e-5)

film    1 cm: N = 3.2594e-03 mol/m2/s  closed form 3.2594e-03  equimolar 3.2081e-03
film   10 cm: N = 3.2594e-04 mol/m2/s  closed form 3.2594e-04  equimolar 3.2081e-04

2. Sphere: Sh -> 2 and the Ranz-Marshall growth#

for Re in (0.0, 1.0, 100.0, 1000.0):
    print(f"Re = {Re:7.1f}: Sh = {float(mt.ranz_marshall(Re, 0.6)):.3f}")

Re =     0.0: Sh = 2.000
Re =     1.0: Sh = 2.506
Re =   100.0: Sh = 7.061
Re =  1000.0: Sh = 18.003

3. Onda film coefficients, before and after the correction#

kw = dict(u_L=0.01, rho_L=1010.0, mu_L=2e-3, D_L=1.5e-9, sigma=0.072)
k_L = mt.onda_liquid_coefficient(**kw, a_p=250.0, d_p=0.05, sigma_c=0.075)
a_w = mt.onda_wetted_area(0.01, 1010.0, 2e-3, 0.072, a_p=250.0, sigma_c=0.075)
print(f"corrected k_L = {float(k_L):.2e} m/s (packed columns: ~1e-4), a_w = {float(a_w):.0f} m2/m3")
print("old difflow_cc value for dry-area Re and (rho g/mu)^(1/3): 0.487 m/s (roughly 4600x too large)")
k_c = mt.onda_gas_coefficient(1.0, 1.2, 1.8e-5, 1.6e-5, a_p=250.0, d_p=0.05)
print(f"k_c = {float(k_c):.3e} m/s")

corrected k_L = 1.06e-04 m/s (packed columns: ~1e-4), a_w = 155 m2/m3
old difflow_cc value for dry-area Re and (rho g/mu)^(1/3): 0.487 m/s (roughly 4600x too large)
k_c = 6.540e-03 m/s

4. Column height by Colburn and by quadrature, and its gradient#

For a straight equilibrium line the analytic Colburn NTU equals the numerical integral; the height H_OG N_OG is then differentiable in L/G, the film coefficients and a_w.

y_in, y_out, m = 0.02, 0.002, 0.8
for lam in (0.5, 1.0, 1.0 + 1e-9):
    Nc = float(mt.ntu_og_colburn(y_in, y_out, 0.0, lam, lam))   # m=lam*L/G ; x_in=0
    Ni = float(mt.ntu_og_integral(y_in, y_out, lambda y: lam * (y - y_out), n=64))
    Nq, _ = quad(lambda y: 1.0 / (y - lam * (y - y_out)), y_out, y_in)
    print(f"lam = {lam:.9f}: Colburn {Nc:.8f}  Gauss-Legendre {Ni:.8f}  quad {Nq:.8f}")

lam = 0.500000000: Colburn 3.40949618  Gauss-Legendre 3.40949618  quad 3.40949618
lam = 1.000000000: Colburn 9.00000000  Gauss-Legendre 9.00000000  quad 9.00000000
lam = 1.000000001: Colburn 9.00000004  Gauss-Legendre 9.00000004  quad 9.00000004

pk = mt.PACKINGS["pall_metal_50"]      # table from memory, unverified
G_m = 40.0
a_w = mt.onda_wetted_area(0.01, 1000.0, 1e-3, 0.072, pk)

def height(LG, scale_kG=1.0, scale_kL=1.0, scale_a=1.0):
    k_y = scale_kG * mt.mass_transfer_coefficient_conversions(
        mt.onda_gas_coefficient(1.0, 1.2, 1.8e-5, 1.6e-5, pk), "c", "y", 300.0, P)
    k_x = scale_kL * mt.k_x_from_k_L(mt.onda_liquid_coefficient(0.01, 1000.0, 1e-3, 1.5e-9, 0.072, pk), 55500.0)
    a = scale_a * a_w
    lam = m / LG
    H = mt.htu_og_from_films(mt.htu_og(G_m, k_y, a), mt.htu_ol(G_m * LG, k_x, a), lam)
    return H * mt.ntu_og_colburn(y_in, y_out, 0.0, m, lam)

print("height at L/G = 1.5:", float(height(1.5)), "m")
# elasticities d ln(Z)/d ln(x): which of k_G, k_L, a_w, L/G dominates?
Z0 = height(1.5)
for name, arg in [("L/G", 0), ("k_G", 1), ("k_L", 2), ("a_w", 3)]:
    g = jax.grad(height, argnums=arg)(1.5, 1.0, 1.0, 1.0)
    base = 1.5 if arg == 0 else 1.0
    print(f"{name:4s}: dlnZ/dlnx = {float(g * base / Z0):+.3f}")

height at L/G = 1.5: 1.7544181631502638 m

L/G : dlnZ/dlnx = -0.583
k_G : dlnZ/dlnx = -0.932
k_L : dlnZ/dlnx = -0.068
a_w : dlnZ/dlnx = -1.000

5. HETP is continuous through lambda = 1#

for lam in (0.99, 0.999999, 1.0, 1.000001, 1.01):
    print(lam, float(mt.hetp_from_htu(0.5, lam)))
print("d HETP / d lam at 1:", float(jax.grad(mt.hetp_from_htu, argnums=1)(0.5, 1.0)), "(expected -H/2 = -0.25)")

99 0.5025167926750721
999999 0.5000002500001667
0 0.5
000001 0.49999975000016667
01 0.49751654265840417
d HETP / d lam at 1: -0.25 (expected -H/2 = -0.25)