""" silica_fluidization.py ====================== Hierarchical agglomerate model for gas fluidization of sub-micron silica spheres (nominal primary diameter 500 nm). WHY AGGLOMERATES ---------------- A 500 nm silica sphere has a van der Waals adhesion force ~10^6-10^7 times its own weight (see `granular_bond_number`). It is deep in Geldart Group C. Such a powder cannot be fluidized as discrete primary particles: it channels, plugs and lifts as a plug. What actually fluidizes is a *hierarchy of agglomerates* (Chaouki 1985; Wang 2002; Zhu 2005; Valverde & Castellanos 2007): primary sphere (d_p ~ 0.5 um) -> simple / "hard" agglomerate (d_s ~ 1-50 um) [level 1] -> complex / "soft" agglomerate (d_c ~ 100-800 um) [level 2] Level 2 is the entity that the gas actually fluidizes. Every hydrodynamic correlation in this module is therefore evaluated on (d_c, rho_c), never on (d_p, rho_p). MODEL CHAIN ----------- 1. Interparticle force : Hamaker vdW, optional Rabinovich asperity correction, optional capillary bridge. 2. Fractal density closure : rho = rho_solid * (d_big/d_small)^(D_f - 3) 3. Agglomerate strength : Rumpf (1962) tensile strength. 4. Break-up balance : sigma_tensile(d_c) = sigma_separating(d_c) -> self-consistent complex-agglomerate size. Separating stress = bed weight gradient, optionally amplified by vibration (Gamma) and by an imposed shear rate (stirrer / jet / sound). 5. Hydrodynamics on (d_c,rho_c): Ergun/Wen-Yu U_mf, Haider-Levenspiel u_t, Richardson-Zaki (Khan-Richardson n) expansion. 6. Diagnostics : operating window, elutriation, DEM feasibility. CALIBRATION ----------- Everything is first-principles except one lumped constant, `k_breakup`, which absorbs the geometry of the break-up event (collision vs. static stress, contact-number statistics, asperity distribution). `calibrate_k_breakup` fits it to a measured U_mf or a measured agglomerate size. Published agglomerate-fluidization models all carry an equivalent constant. UNITS: strict SI throughout (m, kg, s, Pa, N, K). KNOWN LIMITATIONS (read before trusting a number) ------------------------------------------------- 1. `d_simple` is an INPUT, not a prediction. The model closes the level-2 (complex) agglomerate size self-consistently but takes the level-1 size as given. Measure it (SEM, laser diffraction) rather than accepting the default. A consequence: primary-surface roughness changes primary-primary adhesion strongly (`f_vdw_rough`) but barely moves the fluidization result, because inter-agglomerate contact is dominated by primary-sphere curvature (~d_p/2 = 250 nm), which swamps nm-scale asperities. Roughness should physically act by changing `d_simple`; that coupling is NOT modelled. 2. Monodisperse agglomerates. Real beds have a broad, dynamically evolving size distribution; d_c here is a mass-mean equivalent. 3. Homogeneous (APF) expansion is assumed in `bed_voidage`. Nanostructured silica agglomerates usually do expand this way, but a bed that bubbles (ABF regime) will expand far less than Richardson-Zaki predicts. 4. `k_breakup` is fitted, not derived. Every published agglomerate- fluidization model carries an equivalent constant. Fit it to YOUR powder. 5. Electrostatics are ignored. Dry-gas fluidization of fine silica charges readily and tribocharging can dominate cohesion; if experiments disagree with the model, check grounding and humidity before adjusting parameters. 6. Wall effects, distributor design, and gas-channel formation are outside the scope of the 0-D model - that is what the CFD-DEM decks are for. References ---------- Rumpf H. (1962) in Agglomeration, Knepper (ed.), Interscience, 379-418. Ergun S. (1952) Chem. Eng. Prog. 48, 89. Wen C.Y. & Yu Y.H. (1966) AIChE J. 12, 610. Haider A. & Levenspiel O. (1989) Powder Technol. 58, 63. Khan A.R. & Richardson J.F. (1989) Chem. Eng. Commun. 78, 111. Richardson J.F. & Zaki W.N. (1954) Trans. IChemE 32, 35. Rabinovich Y.I. et al. (2000) J. Colloid Interface Sci. 232, 17. Geldart D. (1973) Powder Technol. 7, 285. Chaouki J. et al. (1985) Powder Technol. 43, 117. Zhu C. et al. (2005) AIChE J. 51, 426. Valverde J.M. & Castellanos A. (2007) AIChE J. 53, 1622. """ from __future__ import annotations from dataclasses import dataclass, field, replace import math import numpy as np from scipy.optimize import brentq # ---------------------------------------------------------------------------- # physical constants # ---------------------------------------------------------------------------- G = 9.80665 # m/s^2 R_UNIV = 8.314462 # J/mol/K M_AIR = 0.0289647 # kg/mol # ============================================================================ # 1. GAS PHASE # ============================================================================ @dataclass(frozen=True) class Gas: """Ideal gas with Sutherland viscosity. Defaults = dry air, 20 C, 1 atm.""" T: float = 293.15 # K P: float = 101325.0 # Pa M: float = M_AIR # kg/mol mu_ref: float = 1.716e-5 # Pa.s (air at T_ref) T_ref: float = 273.15 # K S: float = 110.4 # K Sutherland constant for air @property def rho(self) -> float: """Density [kg/m^3].""" return self.P * self.M / (R_UNIV * self.T) @property def mu(self) -> float: """Dynamic viscosity [Pa.s] via Sutherland's law.""" return (self.mu_ref * (self.T / self.T_ref) ** 1.5 * (self.T_ref + self.S) / (self.T + self.S)) @property def nu(self) -> float: """Kinematic viscosity [m^2/s].""" return self.mu / self.rho @property def mean_free_path(self) -> float: """Molecular mean free path [m] (hard-sphere kinetic estimate).""" return (self.mu / self.P) * math.sqrt( math.pi * R_UNIV * self.T / (2.0 * self.M)) def knudsen(d: float, gas: Gas) -> float: """Knudsen number 2*lambda/d. Slip matters when Kn > ~0.01.""" return 2.0 * gas.mean_free_path / d def cunningham(d: float, gas: Gas) -> float: """Cunningham slip correction C_c (Allen-Raabe coefficients).""" kn = knudsen(d, gas) return 1.0 + kn * (1.257 + 0.400 * math.exp(-1.100 / kn)) # ============================================================================ # 2. SOLID / POWDER SPECIFICATION # ============================================================================ @dataclass(frozen=True) class Powder: """Primary-particle and cohesion properties. Defaults describe monodisperse ~500 nm amorphous (Stoeber-type) SiO2. """ d_p: float = 500e-9 # m primary sphere diameter rho_p: float = 2200.0 # kg/m3 skeletal density of amorphous silica A_hamaker: float = 6.5e-20 # J SiO2-air-SiO2 (Bergstrom 1997: 6.3-6.6e-20) z0: float = 4.0e-10 # m vdW cut-off separation rms_roughness: float = 0.0 # m surface asperity rms; 0 = ideal smooth surface_tension: float = 0.0722 # N/m water, for capillary bridges contact_angle: float = 0.0 # rad perfectly wetting hydrophilic silica rel_humidity: float = 0.0 # 0-1; >0 activates capillary bridging # fractal dimensions of the two agglomeration levels Df_simple: float = 2.80 # primaries -> simple agglomerate (dense, compact) Df_complex: float = 2.55 # simple -> complex agglomerate (open, ramified) d_simple: float = 15e-6 # m simple-agglomerate size (measurable by SEM/laser) # mechanics (for DEM deck generation) E_young: float = 70e9 # Pa poisson: float = 0.17 restitution: float = 0.90 friction: float = 0.50 rolling_friction: float = 0.10 # ============================================================================ # 3. INTERPARTICLE FORCES # ============================================================================ def f_vdw_smooth(d1: float, d2: float, pw: Powder) -> float: """van der Waals force between two smooth spheres [N]. F = A * R_eff / (6 z0^2), R_eff = (r1 r2)/(r1 + r2). For equal spheres this reduces to A*d/(24 z0^2). """ r1, r2 = 0.5 * d1, 0.5 * d2 r_eff = r1 * r2 / (r1 + r2) return pw.A_hamaker * r_eff / (6.0 * pw.z0 ** 2) def f_vdw_rough(d1: float, d2: float, pw: Powder, rms: float | None = None) -> float: """vdW force with Rabinovich-type asperity correction [N]. Rabinovich et al. (2000) Eq. 12, two-term asperity model: F = (A/(6 z0^2)) * (R_eff r_a)/(R_eff + r_a) [asperity in contact] + (A R_eff)/(6 (z0 + r_a)^2) [bulk, held off by r_a] with asperity radius r_a = 1.485 * rms. Term 1 is a small sphere (the asperity tip) touching at the vdW cut-off. Term 2 is the main body, now separated by the asperity height instead of z0 - this is what roughness destroys, and it falls off as (z0+r_a)^-2. Limits: r_a -> 0 recovers the smooth result exactly (term1 -> 0, term2 -> F_smooth). The function is NON-MONOTONIC and that is real: a few nm of roughness cuts adhesion by >10x by propping the surfaces apart, but once asperities grow comparable to the particle they act as contacting spheres in their own right and adhesion recovers. The minimum near rms ~ 1-5 nm is the documented adhesion minimum. An earlier formulation here used a term that saturated at the full smooth force for large rms, which wrongly predicted that very rough particles are as sticky as perfectly smooth ones. """ lam = pw.rms_roughness if rms is None else rms if lam <= 0.0: return f_vdw_smooth(d1, d2, pw) r1, r2 = 0.5 * d1, 0.5 * d2 r_eff = r1 * r2 / (r1 + r2) r_a = 1.485 * lam t1 = (r_eff * r_a) / (r_eff + r_a) / pw.z0 ** 2 t2 = r_eff / (pw.z0 + r_a) ** 2 return pw.A_hamaker / 6.0 * (t1 + t2) def f_capillary(d1: float, d2: float, pw: Powder) -> float: """Liquid-bridge (capillary) force between two spheres [N]. Toroidal bridge, small-volume limit: F = 2 pi R_eff gamma (1 + cos theta). Only returned when rel_humidity > 0; above ~60% RH on hydrophilic silica this term DOMINATES vdW and is the usual reason 'dry' runs are irreproducible. """ if pw.rel_humidity <= 0.0: return 0.0 r1, r2 = 0.5 * d1, 0.5 * d2 r_eff = r1 * r2 / (r1 + r2) return 2.0 * math.pi * r_eff * pw.surface_tension * ( 1.0 + math.cos(pw.contact_angle)) def f_cohesive(d1: float, d2: float, pw: Powder, rms: float | None = None) -> float: """Total quasi-static cohesive force [N] = vdW(rough) + capillary.""" return f_vdw_rough(d1, d2, pw, rms) + f_capillary(d1, d2, pw) def particle_weight(d: float, rho: float) -> float: """Gravitational weight of one sphere [N].""" return math.pi / 6.0 * d ** 3 * rho * G def granular_bond_number(pw: Powder, d: float | None = None, rho: float | None = None) -> float: """Bo_g = cohesive force / particle weight. Bo_g >> 1 -> Geldart C, cohesive, will not fluidize as primaries. Bo_g ~ 1 -> fluidization becomes possible. """ d = pw.d_p if d is None else d rho = pw.rho_p if rho is None else rho return f_cohesive(d, d, pw) / particle_weight(d, rho) # ============================================================================ # 4. GELDART CLASSIFICATION # ============================================================================ def geldart_group(d: float, rho_solid: float, gas: Gas) -> str: """Geldart (1973) group from particle size and solid-gas density difference. Geldart drew his lines in mixed units (d in um, rho in g/cm^3): A/B : d * (rho_s - rho_g) = 225 [um * g/cm^3] B/D : d^2 * (rho_s - rho_g) = 1e6 [um^2 * g/cm^3] Converting to SI (d in m, rho in kg/m^3) multiplies the first group by 1e3 and the second by 1e9, giving the thresholds below. Getting this conversion wrong misclassifies every Group B powder as Group A. """ drho = rho_solid - gas.rho ab = drho * d * 1.0e3 # -> um * g/cm^3 bd = drho * d ** 2 * 1.0e9 # -> um^2 * g/cm^3 if d < 30e-6 and ab < 225.0: return "C" if bd > 1.0e6: return "D" if ab < 225.0: return "A" return "B" # ============================================================================ # 5. FRACTAL AGGLOMERATE STRUCTURE # ============================================================================ def fractal_solid_fraction(d_big: float, d_small: float, Df: float) -> float: """Solid volume fraction of a fractal cluster of `d_small` building blocks. phi = (d_big/d_small)^(Df - 3), clipped to (0, 1]. Df -> 3 gives a compact (non-fractal) cluster, phi -> 1. """ if d_big <= d_small: return 1.0 return float(min(1.0, (d_big / d_small) ** (Df - 3.0))) def simple_agglomerate(pw: Powder) -> tuple[float, float, float]: """Level-1 cluster: (d_s [m], rho_s [kg/m3], phi_s [-]).""" phi_s = fractal_solid_fraction(pw.d_simple, pw.d_p, pw.Df_simple) return pw.d_simple, pw.rho_p * phi_s, phi_s def complex_agglomerate(d_c: float, pw: Powder) -> tuple[float, float, float]: """Level-2 cluster: (rho_c [kg/m3], phi_cs [-], eps_c [-]). phi_cs = solid fraction of simple agglomerates inside the complex one. rho_c = rho_p * phi_s * phi_cs (bulk density of the fluidizing entity). """ d_s, rho_s, phi_s = simple_agglomerate(pw) phi_cs = fractal_solid_fraction(d_c, d_s, pw.Df_complex) return pw.rho_p * phi_s * phi_cs, phi_cs, 1.0 - phi_cs # ============================================================================ # 6. AGGLOMERATE STRENGTH AND BREAK-UP # ============================================================================ def rumpf_tensile_strength(phi: float, d_block: float, f_bond: float) -> float: """Rumpf (1962) tensile strength of a random packing of bonded spheres [Pa]. sigma_T = (phi / (1 - phi)) * F / d_block^2 using Rumpf's coordination-number closure k * eps ~ pi. """ phi = float(np.clip(phi, 1e-6, 1.0 - 1e-9)) return (phi / (1.0 - phi)) * f_bond / d_block ** 2 def effective_hamaker(A_bulk: float, phi: float) -> float: """Hamaker constant of a porous body, A_eff = A_bulk * phi^2. Standard two-body dilution scaling: the vdW interaction is proportional to the product of the number densities of polarizable matter in each body. """ return A_bulk * phi ** 2 def bond_between_simple_agglomerates(pw: Powder) -> float: """Cohesive force holding two simple agglomerates together [N]. Two effects weaken it far below the primary-primary force: (a) the bodies are porous -> A_eff = A * phi_s^2 (b) contact is made through surface primaries, i.e. the effective asperity scale is d_p/2, not the molecular cut-off. """ d_s, _, phi_s = simple_agglomerate(pw) pw_eff = replace(pw, A_hamaker=effective_hamaker(pw.A_hamaker, phi_s), rms_roughness=max(pw.rms_roughness, 0.5 * pw.d_p)) return f_cohesive(d_s, d_s, pw_eff) def separating_stress(d_c: float, pw: Powder, gas: Gas, vib_gamma: float = 0.0, shear_rate: float = 0.0, bed_solid_fraction: float = 0.45) -> float: """Stress available to tear a complex agglomerate apart [Pa]. Three additive contributions: inertial/gravitational : rho_bed * (1 + Gamma) * g * d_c (Gamma = dimensionless vibration acceleration) viscous shear : mu_g * shear_rate (stirrer, jet, acoustic field) dynamic pressure : 0.5 * rho_g * (shear_rate * d_c)^2 """ rho_c, _, _ = complex_agglomerate(d_c, pw) rho_bed = bed_solid_fraction * rho_c s_grav = rho_bed * (1.0 + vib_gamma) * G * d_c s_shear = gas.mu * shear_rate s_dyn = 0.5 * gas.rho * (shear_rate * d_c) ** 2 return s_grav + s_shear + s_dyn def equilibrium_agglomerate_size(pw: Powder, gas: Gas, k_breakup: float = 1.0, vib_gamma: float = 0.0, shear_rate: float = 0.0, bed_solid_fraction: float = 0.45, bracket: tuple[float, float] = (2e-6, 0.05) ) -> float: """Self-consistent complex-agglomerate diameter d_c [m]. Solves k_breakup * sigma_separating(d_c) = sigma_tensile(d_c). Below the root the cluster is stronger than the flow can tear -> it grows. Above it, the cluster is torn apart -> it shrinks. The root is the stable dynamic-equilibrium size. `k_breakup` lumps the break-up geometry (collisional amplification, contact statistics). Fit it with `calibrate_k_breakup`. """ d_s, _, _ = simple_agglomerate(pw) f_bond = bond_between_simple_agglomerates(pw) def residual(d_c: float) -> float: _, phi_cs, _ = complex_agglomerate(d_c, pw) sig_t = rumpf_tensile_strength(phi_cs, d_s, f_bond) sig_s = separating_stress(d_c, pw, gas, vib_gamma, shear_rate, bed_solid_fraction) return k_breakup * sig_s - sig_t lo, hi = max(bracket[0], 1.05 * d_s), bracket[1] if residual(lo) > 0.0: return lo # even the smallest cluster is torn apart if residual(hi) < 0.0: return math.nan # unbounded growth -> defluidization / plugging return float(brentq(residual, lo, hi, xtol=1e-10, rtol=1e-10)) def calibrate_k_breakup(pw: Powder, gas: Gas, *, target_d_c: float | None = None, target_u_mf: float | None = None, eps_mf: float = 0.55, vib_gamma: float = 0.0, shear_rate: float = 0.0, bed_solid_fraction: float = 0.45) -> float: """Fit `k_breakup` to a measured agglomerate size or a measured U_mf.""" if (target_d_c is None) == (target_u_mf is None): raise ValueError("give exactly one of target_d_c or target_u_mf") if target_u_mf is not None: def size_err(d_c: float) -> float: rho_c, _, _ = complex_agglomerate(d_c, pw) return u_mf(d_c, rho_c, gas, eps_mf=eps_mf) - target_u_mf target_d_c = float(brentq(size_err, 2e-6, 0.05, xtol=1e-12, rtol=1e-12)) def k_err(k: float) -> float: d = equilibrium_agglomerate_size(pw, gas, k, vib_gamma, shear_rate, bed_solid_fraction) return (1e6 if math.isnan(d) else d) - target_d_c return float(brentq(k_err, 1e-3, 1e5, xtol=1e-10, rtol=1e-10)) # ============================================================================ # 7. HYDRODYNAMICS (evaluated on the AGGLOMERATE, not the primary particle) # ============================================================================ def archimedes(d: float, rho_s: float, gas: Gas) -> float: """Ar = rho_g (rho_s - rho_g) g d^3 / mu^2.""" return gas.rho * (rho_s - gas.rho) * G * d ** 3 / gas.mu ** 2 def u_mf(d: float, rho_s: float, gas: Gas, eps_mf: float = 0.55, sphericity: float = 1.0) -> float: """Minimum fluidization velocity [m/s] from the Ergun force balance. Setting the Ergun pressure gradient equal to the buoyant bed weight gives a quadratic in Re_mf: (1.75/(eps^3 phi)) Re^2 + (150(1-eps)/(eps^3 phi^2)) Re - Ar = 0 Solved exactly (no Wen-Yu constant-eps approximation). """ Ar = archimedes(d, rho_s, gas) a = 1.75 / (eps_mf ** 3 * sphericity) b = 150.0 * (1.0 - eps_mf) / (eps_mf ** 3 * sphericity ** 2) Re = (-b + math.sqrt(b * b + 4.0 * a * Ar)) / (2.0 * a) return Re * gas.mu / (gas.rho * d) def u_mf_wen_yu(d: float, rho_s: float, gas: Gas) -> float: """Wen & Yu (1966) correlation. CROSS-CHECK ONLY - see caveat. Wen-Yu's constants (33.7, 0.0408) come from fitting the Ergun balance with a FIXED voidage. Solving eps^3/(150(1-eps)) = 1/1650 shows the embedded value is eps_mf = 0.383, i.e. a dense packing of hard grains. Verified numerically: `u_mf` reproduces this function to 0.1% when called with eps_mf = 0.383, which confirms both implementations. But agglomerate beds run at eps_mf = 0.55-0.75, where Wen-Yu under-predicts U_mf by 4-12x. Use `u_mf` with the true voidage for anything quantitative here. """ Ar = archimedes(d, rho_s, gas) Re = math.sqrt(33.7 ** 2 + 0.0408 * Ar) - 33.7 return Re * gas.mu / (gas.rho * d) def u_terminal(d: float, rho_s: float, gas: Gas, sphericity: float = 1.0) -> float: """Terminal settling velocity [m/s], Haider & Levenspiel (1989).""" Ar = archimedes(d, rho_s, gas) d_star = Ar ** (1.0 / 3.0) u_star = (18.0 / d_star ** 2 + (2.335 - 1.744 * sphericity) / math.sqrt(d_star)) ** -1 return u_star * (gas.rho ** 2 / (gas.mu * (rho_s - gas.rho) * G)) ** (-1 / 3) def richardson_zaki_n(d: float, rho_s: float, gas: Gas, D_bed: float = 0.05) -> float: """Richardson-Zaki index n from the Khan & Richardson (1989) correlation.""" Ar = archimedes(d, rho_s, gas) rhs = 0.043 * Ar ** 0.57 * (1.0 - 1.24 * (d / D_bed) ** 0.27) # (4.8 - n)/(n - 2.4) = rhs -> n = (4.8 + 2.4*rhs)/(1 + rhs) return (4.8 + 2.4 * rhs) / (1.0 + rhs) def bed_voidage(u: float, d: float, rho_s: float, gas: Gas, eps_mf: float = 0.55, D_bed: float = 0.05) -> float: """Expanded-bed voidage from Richardson-Zaki, eps = (u/u_t)^(1/n). Valid for homogeneous (particulate) expansion, which is exactly the APF regime that nanostructured silica agglomerates exhibit. Clipped at eps_mf below minimum fluidization. """ ut = u_terminal(d, rho_s, gas) n = richardson_zaki_n(d, rho_s, gas, D_bed) return float(np.clip((u / ut) ** (1.0 / n), eps_mf, 0.9995)) def bed_height(u: float, d: float, rho_s: float, gas: Gas, H_mf: float, eps_mf: float = 0.55, D_bed: float = 0.05) -> float: """Expanded bed height [m] from solids-volume conservation.""" eps = bed_voidage(u, d, rho_s, gas, eps_mf, D_bed) return H_mf * (1.0 - eps_mf) / (1.0 - eps) def ergun_pressure_drop(u: float, d: float, gas: Gas, eps: float, H: float, sphericity: float = 1.0) -> float: """Ergun (1952) packed-bed pressure drop [Pa].""" visc = (150.0 * gas.mu * (1.0 - eps) ** 2 * u / (eps ** 3 * (sphericity * d) ** 2)) inert = (1.75 * gas.rho * (1.0 - eps) * u ** 2 / (eps ** 3 * sphericity * d)) return (visc + inert) * H def pressure_drop_curve(u_array, d_c: float, rho_c: float, gas: Gas, H_mf: float, eps_mf: float = 0.55, D_bed: float = 0.05, cohesive_overshoot: float = 0.0): """Delta-p vs superficial velocity for the whole bed [Pa]. Fixed-bed branch = Ergun. Fluidized branch = buoyant bed weight per area (velocity independent). `cohesive_overshoot` (0-1) adds the classic Group-C pressure peak on the increasing-velocity branch: cohesive powders must have their inter-agglomerate bonds ruptured before the bed unlocks, so dp overshoots the bed weight by (1 + overshoot) before collapsing onto it. Note that on the fluidized branch dp is exactly velocity-independent: the bed expands (eps rises, H rises) but the product (1-eps)*H is fixed by solids conservation, so dp always equals the buoyant weight per unit area. With cohesive_overshoot > 0 the increasing-u branch stays PACKED past the nominal U_mf - the gas must first rupture the inter-agglomerate bonds - so the Ergun line continues to climb to (1+overshoot)*dp_fluid before the bed unlocks and collapses onto the plateau. That peak, and the hysteresis against the decreasing-u branch, is the experimental fingerprint of a cohesive (Group C / agglomerate) powder. """ u_array = np.atleast_1d(np.asarray(u_array, dtype=float)) dp_fluid = (1.0 - eps_mf) * (rho_c - gas.rho) * G * H_mf umf = u_mf(d_c, rho_c, gas, eps_mf=eps_mf) dp_peak = dp_fluid * (1.0 + cohesive_overshoot) dp = np.empty_like(u_array) for i, u in enumerate(u_array): dp_packed = ergun_pressure_drop(u, d_c, gas, eps_mf, H_mf) if cohesive_overshoot > 0.0: # stay on the packed branch until the peak stress is reached dp[i] = dp_packed if dp_packed < dp_peak else dp_fluid else: dp[i] = dp_packed if u < umf else dp_fluid return dp, dp_fluid, umf # ============================================================================ # 8. FULL CASE SOLVE # ============================================================================ @dataclass class Bed: """Column and inventory specification.""" D: float = 0.05 # m column inner diameter H_mf: float = 0.10 # m bed height at minimum fluidization eps_mf: float = 0.55 # - voidage between complex agglomerates at U_mf @dataclass class Case: powder: Powder = field(default_factory=Powder) gas: Gas = field(default_factory=Gas) bed: Bed = field(default_factory=Bed) k_breakup: float = 1.0 vib_gamma: float = 0.0 # dimensionless vibration acceleration A*omega^2/g shear_rate: float = 0.0 # 1/s, imposed by stirrer / jet / acoustic field def solve(case: Case) -> dict: """Run the full model chain and return every diagnostic.""" pw, gas, bed = case.powder, case.gas, case.bed # --- primary-particle diagnostics ------------------------------------- bo_smooth = f_vdw_smooth(pw.d_p, pw.d_p, pw) / particle_weight(pw.d_p, pw.rho_p) bo_total = granular_bond_number(pw) out = { "d_p": pw.d_p, "rho_p": pw.rho_p, "gas_rho": gas.rho, "gas_mu": gas.mu, "F_vdw_primary": f_vdw_smooth(pw.d_p, pw.d_p, pw), "F_cohesive_primary": f_cohesive(pw.d_p, pw.d_p, pw), "W_primary": particle_weight(pw.d_p, pw.rho_p), "Bo_smooth": bo_smooth, "Bo_total": bo_total, "Knudsen_primary": knudsen(pw.d_p, gas), "geldart_primary": geldart_group(pw.d_p, pw.rho_p, gas), "u_mf_primary_hypothetical": u_mf(pw.d_p, pw.rho_p, gas, bed.eps_mf), } # --- agglomerate structure -------------------------------------------- d_s, rho_s, phi_s = simple_agglomerate(pw) out.update({"d_simple": d_s, "rho_simple": rho_s, "phi_simple": phi_s, "F_bond_simple_simple": bond_between_simple_agglomerates(pw)}) d_c = equilibrium_agglomerate_size(pw, gas, case.k_breakup, case.vib_gamma, case.shear_rate, 1.0 - bed.eps_mf) out["d_complex"] = d_c if math.isnan(d_c): out["status"] = ("UNBOUNDED AGGLOMERATION - cohesion exceeds every " "available break-up stress; bed will channel or lift " "as a plug. Add assistance (vibration/stirring/sound).") return out rho_c, phi_cs, eps_c = complex_agglomerate(d_c, pw) n_prim = (rho_c / pw.rho_p) * (d_c / pw.d_p) ** 3 out.update({ "rho_complex": rho_c, "phi_complex_in_simple": phi_cs, "porosity_complex": 1.0 - rho_c / pw.rho_p, "n_primaries_per_complex": n_prim, "geldart_complex": geldart_group(d_c, rho_c, gas), "Bo_complex": (bond_between_simple_agglomerates(pw) / particle_weight(d_c, rho_c)), }) # --- hydrodynamics on the agglomerate --------------------------------- umf = u_mf(d_c, rho_c, gas, bed.eps_mf) ut = u_terminal(d_c, rho_c, gas) out.update({ "Ar_complex": archimedes(d_c, rho_c, gas), "u_mf": umf, "u_mf_wen_yu": u_mf_wen_yu(d_c, rho_c, gas), "u_t": ut, "Re_mf": gas.rho * umf * d_c / gas.mu, "turndown_ut_over_umf": ut / umf, "n_RZ": richardson_zaki_n(d_c, rho_c, gas, bed.D), "dp_fluidized": (1.0 - bed.eps_mf) * (rho_c - gas.rho) * G * bed.H_mf, "bed_mass": (math.pi / 4 * bed.D ** 2 * bed.H_mf * (1.0 - bed.eps_mf) * rho_c), "H_at_2umf": bed_height(2 * umf, d_c, rho_c, gas, bed.H_mf, bed.eps_mf, bed.D), "H_at_5umf": bed_height(5 * umf, d_c, rho_c, gas, bed.H_mf, bed.eps_mf, bed.D), }) # --- DEM feasibility --------------------------------------------------- V_bed = math.pi / 4 * bed.D ** 2 * bed.H_mf out["N_primaries_in_bed"] = (V_bed * (1 - bed.eps_mf) * rho_c / pw.rho_p) / (math.pi / 6 * pw.d_p ** 3) out["N_agglomerates_in_bed"] = (V_bed * (1 - bed.eps_mf) / (math.pi / 6 * d_c ** 3)) out["dt_rayleigh_complex"] = rayleigh_timestep(d_c, rho_c, pw) out["status"] = "OK" return out def rayleigh_timestep(d: float, rho: float, pw: Powder, safety: float = 0.2) -> float: """Rayleigh critical DEM time step [s] (surface-wave criterion).""" Gmod = pw.E_young / (2.0 * (1.0 + pw.poisson)) return safety * (math.pi * 0.5 * d * math.sqrt(rho / Gmod) / (0.1631 * pw.poisson + 0.8766)) def sjkr_cohesion_energy_density(pw: Powder, d_c: float, E_sim: float | None = None) -> float: """LIGGGHTS `cohesionEnergyDensity` k_c [J/m^3] for the SJKR model. SJKR adds F_coh = k_c * A_contact, with A_contact the Hertzian contact area. Matching that to the physical inter-agglomerate bond force gives k_c = F_bond / A_contact. MUST be evaluated at the Young's modulus the SIMULATION will actually use. DEM runs almost always soften the contact (see `stiffness_scaled_powder`) to afford a usable time step. Because A_contact ~ E^(-2/3), a k_c computed at the physical modulus but applied at a softened one inflates the cohesive force by exactly (E_true/E_sim)^(2/3) - silently turning a fluidizable bed into a solid plug. For 70 GPa -> 10 MPa that factor is 366x (not ~80x, which is the 70 GPa -> 100 MPa case). Use `cohesion_scaling_error` rather than quoting a remembered number. Pass E_sim to keep the *force* invariant. """ E = pw.E_young if E_sim is None else E_sim f_bond = bond_between_simple_agglomerates(pw) rho_c, _, _ = complex_agglomerate(d_c, pw) E_star = E / (2.0 * (1.0 - pw.poisson ** 2)) R_star = d_c / 4.0 F_load = max(particle_weight(d_c, rho_c), 1e-18) a_contact = (3.0 * F_load * R_star / (4.0 * E_star)) ** (1.0 / 3.0) A_contact = math.pi * a_contact ** 2 return f_bond / max(A_contact, 1e-24) def cohesion_scaling_error(pw: Powder, E_sim: float, E_true: float | None = None) -> float: """Factor by which cohesion is overstated if k_c is NOT rescaled. Applying a k_c computed at `E_true` in a run at `E_sim` multiplies the SJKR cohesive force by (E_true/E_sim)^(2/3), since F = k_c * A_contact and the Hertzian area scales as E^(-2/3). Returns that factor so documentation and generated decks quote a computed number instead of a recalled one. """ E_true = pw.E_young if E_true is None else E_true return (E_true / E_sim) ** (2.0 / 3.0) def stiffness_scaled_powder(pw: Powder, E_sim: float) -> Powder: """Return a copy of `pw` with a reduced Young's modulus for DEM. Standard DEM practice: the Rayleigh time step scales as E^(-1/2), so running silica at its true 70 GPa costs ~1.1e8 steps per simulated second. Softening to 1e7-1e8 Pa buys 1-2 orders of magnitude and leaves bulk behaviour essentially unchanged PROVIDED the contacts stay stiff relative to the driving stresses - check with `stiffness_adequacy`. """ return replace(pw, E_young=E_sim) def stiffness_adequacy(pw: Powder, d_c: float, E_sim: float) -> dict: """Check that a softened contact is still stiff enough to be meaningful. Returns the maximum Hertzian overlap under the bed's own weight as a fraction of particle diameter. Keep it below ~1%: beyond that the particles interpenetrate visibly and the packing structure is corrupted. """ rho_c, _, _ = complex_agglomerate(d_c, pw) E_star = E_sim / (2.0 * (1.0 - pw.poisson ** 2)) R_star = d_c / 4.0 F = particle_weight(d_c, rho_c) delta = (9.0 * F ** 2 / (16.0 * E_star ** 2 * R_star)) ** (1.0 / 3.0) return {"overlap_m": delta, "overlap_fraction_of_d": delta / d_c, "acceptable": bool(delta / d_c < 0.01), "dt_rayleigh": rayleigh_timestep(d_c, rho_c, replace(pw, E_young=E_sim))}