#!/usr/bin/env python3 """ Self-contained driver for SR-2026-001. Regenerates every number in the continuum-model tables of the report from the input parameters alone. No imports beyond the standard library, NumPy and SciPy. python3 driver.py Verified against the report with NumPy 2.5.1 / SciPy 1.18.0 on Python 3.12.13. """ import math from scipy.optimize import brentq G = 9.80665 # ---------------------------------------------------------------- inputs d_p = 5.00e-7 # primary diameter [m] rho_p = 2200.0 # skeletal density [kg/m^3] A_ham = 6.50e-20 # Hamaker constant, silica across air [J] z0 = 4.00e-10 # contact separation [m] E_true = 7.00e10 # Young modulus [Pa] nu = 0.17 # Poisson ratio [-] d_s = 1.50e-5 # simple-agglomerate diameter (INPUT) [m] Df1 = 2.80 # level-1 fractal dimension [-] Df2 = 2.55 # level-2 fractal dimension [-] rho_g = 1.2041 # gas density (ideal gas, 293.15 K) [kg/m^3] mu_g = 1.8135e-5 # gas viscosity (Sutherland, 293.15 K) [Pa s] eps_mf = 0.65 # voidage at minimum fluidization [-] D_bed = 0.05 # column diameter [m] H_mf = 0.10 # settled bed height [m] E_sim = 1.0e7 # softened DEM contact stiffness [Pa] d_c_target = 250e-6 # calibration target for d_c [m] # ------------------------------------------------------- structure/cohesion def f_vdw(d1, d2, A=A_ham): """van der Waals force between two smooth spheres in contact [N].""" r_eff = (0.5 * d1) * (0.5 * d2) / (0.5 * d1 + 0.5 * d2) return A * r_eff / (6.0 * z0 ** 2) def f_vdw_rough(d1, d2, A, rms): """Rabinovich et al. (2000) two-term asperity model [N]. Term 1 is the asperity tip in contact at the van der Waals cut-off; term 2 is the bulk body, propped away by the asperity height. rms -> 0 recovers the smooth result exactly. """ if rms <= 0.0: return f_vdw(d1, d2, A) r1, r2 = 0.5 * d1, 0.5 * d2 r_eff = r1 * r2 / (r1 + r2) r_a = 1.485 * rms t1 = (r_eff * r_a) / (r_eff + r_a) / z0 ** 2 t2 = r_eff / (z0 + r_a) ** 2 return A / 6.0 * (t1 + t2) def weight(d, rho): return (math.pi / 6.0) * d ** 3 * rho * G phi_s = (d_s / d_p) ** (Df1 - 3.0) # solid fraction of simple agglomerate rho_s = phi_s * rho_p def complex_agglomerate(d_c): """(solid fraction of simple agglomerates inside, effective density).""" phi_cs = (d_c / d_s) ** (Df2 - 3.0) return phi_cs, rho_p * phi_s * phi_cs # Bond between two simple agglomerates. Two effects weaken it far below the # primary-primary force: the bodies are porous, so A_eff = A*phi_s^2, and the # contact is made through a surface primary, so the asperity scale is d_p/2 # rather than the molecular cut-off. F_bond = f_vdw_rough(d_s, d_s, A_ham * phi_s ** 2, 0.5 * d_p) def sigma_T(d_c): """Rumpf tensile strength of the complex agglomerate [Pa].""" phi_c, _ = complex_agglomerate(d_c) return phi_c / (1.0 - phi_c) * F_bond / d_s ** 2 def sigma_S(d_c, solid_frac, vib_gamma=0.0): """Separating stress available in the bed: bed-weight gradient [Pa]. Buoyancy is not subtracted here: the bed bulk density is solid_frac * rho_c, and the gas contributes negligibly at these densities. """ _, rho_c = complex_agglomerate(d_c) return solid_frac * rho_c * (1.0 + vib_gamma) * G * d_c def d_complex(k, solid_frac, lo=None, hi=0.20): """Root of k*sigma_S = sigma_T. Bracket is widened until it straddles.""" lo = 1.05 * d_s if lo is None else lo f = lambda x: k * sigma_S(x, solid_frac) - sigma_T(x) return brentq(f, lo, hi, xtol=1e-16, rtol=8.9e-16) def calibrate_k(target, solid_frac): """Secant search on log k so that d_complex(k) == target.""" f = lambda lk: d_complex(math.exp(lk), solid_frac) - target return math.exp(brentq(f, math.log(0.5), math.log(1e3), xtol=1e-14)) # ---------------------------------------------------------- hydrodynamics def archimedes(d, rho): return rho_g * (rho - rho_g) * G * d ** 3 / mu_g ** 2 def u_mf(d, rho, eps): """Exact root of the Ergun balance written in Re.""" a = 1.75 / eps ** 3 b = 150.0 * (1.0 - eps) / eps ** 3 Ar = archimedes(d, rho) Re = (-b + math.sqrt(b * b + 4.0 * a * Ar)) / (2.0 * a) return Re * mu_g / (rho_g * d) def u_terminal(d, rho): """Haider-Levenspiel (1989) for a sphere.""" dstar = d * (rho_g * (rho - rho_g) * G / mu_g ** 2) ** (1.0 / 3.0) ustar = (18.0 / dstar ** 2 + 0.591 / math.sqrt(dstar)) ** -1.0 return ustar / (rho_g ** 2 / (mu_g * (rho - rho_g) * G)) ** (1.0 / 3.0) def rz_index(d, rho, D=D_bed): """Khan & Richardson (1989) Richardson-Zaki index, with wall correction.""" Ar = archimedes(d, rho) rhs = 0.043 * Ar ** 0.57 * (1.0 - 1.24 * (d / D) ** 0.27) return (4.8 + 2.4 * rhs) / (1.0 + rhs) # ------------------------------------------------------------ DEM mapping def rayleigh_dt(d, rho, E, safety=0.2): """Rayleigh surface-wave step [s]. safety=1.0 returns the raw stability limit; the production runs use safety=0.2, i.e. one fifth of that limit. """ Gmod = E / (2.0 * (1.0 + nu)) return safety * math.pi * 0.5 * d * math.sqrt(rho / Gmod) / (0.1631 * nu + 0.8766) def hertz_area(d_c, rho_c, E): E_star = E / (2.0 * (1.0 - nu ** 2)) R_star = d_c / 4.0 a = (3.0 * weight(d_c, rho_c) * R_star / (4.0 * E_star)) ** (1.0 / 3.0) return math.pi * a ** 2 def jkr_gamma(d_c): """Surface energy reproducing F_bond as the JKR pull-off force.""" return F_bond / (1.5 * math.pi * (d_c / 4.0)) def tabor(d_c, gamma, E): E_star = E / (2.0 * (1.0 - nu ** 2)) return ((d_c / 4.0) * gamma ** 2 / (E_star ** 2 * z0 ** 3)) ** (1.0 / 3.0) # ------------------------------------------------------------------- main if __name__ == "__main__": solid_frac = 1.0 - eps_mf k = calibrate_k(d_c_target, solid_frac) d_c = d_complex(k, solid_frac) phi_c, rho_c = complex_agglomerate(d_c) Umf = u_mf(d_c, rho_c, eps_mf) ut = u_terminal(d_c, rho_c) gam = jkr_gamma(d_c) print(f"Bo_primary {f_vdw(d_p,d_p)/weight(d_p,rho_p):.4g}") print(f"F_vdW {f_vdw(d_p,d_p):.4g} N") print(f"W_primary {weight(d_p,rho_p):.4g} N") print(f"phi_s {phi_s:.4g}") print(f"k (fitted) {k:.4f}") print(f"d_c {d_c*1e6:.4g} um") print(f"rho_c {rho_c:.4g} kg/m3") print(f"porosity {(1-phi_c)*100:.3f} %") print(f"Bo_complex {F_bond/weight(d_c,rho_c):.4g}") print(f"Ar {archimedes(d_c,rho_c):.4g}") print(f"U_mf {Umf*1e3:.4g} mm/s") print(f"Re_mf {rho_g*Umf*d_c/mu_g:.4g}") print(f"u_t {ut*1e3:.4g} mm/s") print(f"turndown {ut/Umf:.4g}") print(f"n_RZ {rz_index(d_c,rho_c):.4g}") print(f"dp_bed {solid_frac*(rho_c-rho_g)*G*H_mf:.4g} Pa") print(f"F_bond {F_bond:.4g} N") print(f"gamma_JKR {gam:.4g} J/m2") print(f"Tabor mu {tabor(d_c,gam,E_sim):.4g}") print(f"dt_R raw (70 GPa) {rayleigh_dt(d_c,rho_c,E_true,safety=1.0):.4g} s") print(f"dt_R raw (10 MPa) {rayleigh_dt(d_c,rho_c,E_sim,safety=1.0):.4g} s") print(f"dt production (x0.2) {rayleigh_dt(d_c,rho_c,E_sim):.4g} s") print(f"dt_R 0.2x (70 GPa) {rayleigh_dt(d_c,rho_c,E_true):.4g} s") print(f"k_c at E_sim {F_bond/hertz_area(d_c,rho_c,E_sim):.4g} J/m3") print(f"chi {(E_true/E_sim)**(2.0/3.0):.4g}")