# Chemical Unit Operations This document provides comprehensive documentation for all chemical unit operations available in Difflow. --- (reactors)= ## Reactors (cstr-continuous-stirred-tank-reactor)= ### CSTR (Continuous Stirred-Tank Reactor) **Location**: `difflow/units/cstr.py` **Class**: `CSTR` **Description**: Models an ideal continuous stirred-tank reactor with perfect mixing. The reactor contents are assumed to be at uniform temperature and composition, equal to the outlet conditions. #### Process Role CSTRs are widely used in chemical processes for: - Liquid-phase reactions - Polymerization reactions - Fermentation (as idealized model) - Processes requiring uniform conditions #### Parameters ```python @dataclass class CSTRParams: V: float # Reactor volume (m³) rate_fn: Callable # Rate function: rate_fn(C, T, rate_params) -> r [n_reactions] stoich: Array # Stoichiometric matrix [n_species × n_reactions] rate_params: dict # Parameters passed to rate_fn (k_ref, E_a, ...) species_order: list # Species names, in the row order of `stoich` dH_rxn: Array = None # Heat of reaction for each reaction (J/mol) molar_density: float = None # Constant molar density (mol/m³), see below eos: Any = None # Cubic EOS for the reaction-phase density reaction_phase: str = None # 'liquid' or 'vapor' (required with eos) T_damping: float = ... # Damping on the adiabatic/duty temperature solve ``` #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `inlet` | Stream | - | Inlet stream with species flows, T, P | | `T_spec` | float | K | Target outlet temperature (isothermal mode) | | `Q_spec` | float | W | Specified heat duty (specified_duty mode) | | `volumetric_flow` | float | m³/s | Volumetric flow rate (optional) | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `outlet` | Stream | - | Outlet stream | | `info['Q']` | float | W | Heat duty (positive = heating) | | `info['rates']` | Array | mol/m³/s | Reaction rates | | `info['conversion']` | dict | - | Conversion of each species | #### Operating Modes The mode is chosen when the unit is built (`CSTR(params, thermo, mode=...)`): 1. **Isothermal** (`mode='isothermal'`, `T_spec` given at call time): Outlet temperature is fixed, heat duty calculated 2. **Adiabatic** (`mode='adiabatic'`): Q = 0, outlet temperature calculated 3. **Specified Duty** (`mode='specified_duty'`, `Q_spec` given at call time): Heat duty fixed, outlet temperature calculated #### Governing Equations **Material Balance** (steady-state): $$F_{i,out} = F_{i,in} + V \sum_j \nu_{ij} r_j$$ Where: - $F_{i,out}$: Outlet molar flow of species i (mol/s) - $F_{i,in}$: Inlet molar flow of species i (mol/s) - $V$: Reactor volume (m³) - $\nu_{ij}$: Stoichiometric coefficient of species i in reaction j - $r_j$: Rate of reaction j (mol/m³/s) **Reaction Rate** (Arrhenius kinetics): $$r_j = k_{ref} \exp\left[\frac{E_a}{R}\left(\frac{1}{T_{ref}} - \frac{1}{T}\right)\right] \prod_i C_i^{n_i}$$ Where: - $k_{ref}$: Rate constant at reference temperature - $E_a$: Activation energy (J/mol) - $R$: Gas constant (8.314 J/mol/K) - $C_i$: Concentration of species i (mol/m³) - $n_i$: Reaction order with respect to species i **Energy Balance**: $$Q = \dot{H}_{out} - \dot{H}_{in} + V \sum_j r_j \Delta H_{rxn,j}$$ Where: - $Q$: Heat duty (W) - $\dot{H}$: Enthalpy flow rate (W) - $\Delta H_{rxn,j}$: Heat of reaction j (J/mol) **Conversion**: $$X = \frac{F_{A,in} - F_{A,out}}{F_{A,in}}$$ #### Example Usage ```python from difflow import CSTR, CSTRParams, IdealThermo, SpeciesData, make_stream import jax.numpy as jnp # Two-species thermodynamics (liquid Cp, Antoine and Hvap constants) species_data = { n: SpeciesData(n, MW=100.0, Cp_coeffs=(150.0, 0.0, 0.0, 0.0), Hvap_coeffs=(35000.0, 0.38, 600.0), antoine_coeffs=(10.0, 3000.0, -50.0)) for n in ('A', 'B') } thermo = IdealThermo(species_data) # A -> B (first-order, exothermic) def rate_fn(C, T, p): k = p['k_ref'] * jnp.exp(p['E_a'] / 8.314 * (1 / p['T_ref'] - 1 / T)) return jnp.array([k * C['A']]) params = CSTRParams( V=2.0, # m³ rate_fn=rate_fn, stoich=jnp.array([[-1.0], [1.0]]), # A -> B rate_params={'k_ref': 1e-3, 'E_a': 50000.0, 'T_ref': 350.0}, # 1/s at 350 K species_order=['A', 'B'], dH_rxn=jnp.array([-80000.0]), # J/mol (exothermic) molar_density=1000.0, # mol/m³ ) inlet = make_stream({'A': 1.0, 'B': 0.0}, T=350.0, P=101325.0) # Isothermal operation: the mode is set on the unit cstr = CSTR(params, thermo, mode='isothermal') outlet, info = cstr(inlet, T_spec=350.0) print(f"Conversion of A: {info['conversion']['A']:.2%}") print(f"Heat duty: {info['Q']:.2f} W") # Adiabatic operation cstr_adiab = CSTR(params, thermo, mode='adiabatic') outlet_adiab, info_adiab = cstr_adiab(inlet) print(f"Outlet temperature: {outlet_adiab['T']:.1f} K") ``` #### Concentration Basis (Molar Density) The rate law is evaluated at concentrations, so the reactor needs a molar density: $C_i = F_i / \dot{V}$ with $\dot{V} = F_{total}/\rho$, which makes the residence time $\tau = V\rho/F_{total}$. An error in $\rho$ is a proportional error in $\tau$, and so in the conversion. Three ways to set it, in the order the CSTR resolves them: 1. **An equation of state** -- `eos=` with `reaction_phase='liquid'` or `'vapor'`. Concentration is then the real molarity at reactor $(T, P, y)$, $C_i = y_i\,\rho_{EOS}(T, P, y)$, recomputed inside the solve, so the reactor shares the flash's thermodynamics. `reaction_phase` is required with `eos`: the liquid and vapor molar densities differ by two orders of magnitude, so there is no defensible default. Pair it with a `CubicThermo` to make the *enthalpy* real-gas too; a `CubicThermo` passed as `thermo` also supplies the EOS itself when `reaction_phase` is set and no `eos=` is given. 2. **A constant** -- `molar_density=`. 3. **Neither**, in which case the reactor falls back to 55500 mol/m^3 (liquid water) and raises a `CSTRDensityWarning`. That fallback is right only for aqueous systems: it is ~8x high for a C3-C8 hydrocarbon liquid and ~200x high for a gas, and it silently inflates residence time. Treat the warning as a request to say which basis you meant. Passing `volumetric_flow=` to the call sets $\dot{V}$ outright and bypasses all three; `info['molar_density']` then reports the density that flow implies. ```python from difflow.eos import PengRobinson from difflow.database import get_critical_props species = ['propane', 'n-butane'] eos = PengRobinson({c: get_critical_props(c) for c in species}) params_eos = CSTRParams(V=2.0, rate_fn=lambda C, T, p: jnp.array([p['k'] * C['propane']]), stoich=jnp.array([[-1.0], [1.0]]), rate_params={'k': 1e-3}, species_order=species, eos=eos, reaction_phase='liquid') ``` #### Design Considerations - **Residence Time**: $\tau = V/Q_{vol}$ determines conversion - **Heat Transfer**: Large exothermic reactions may require cooling coils or jackets - **Mixing**: Perfect mixing assumption requires adequate agitation - **Multiple CSTRs**: Series arrangement approaches PFR behavior --- (pfr-plug-flow-reactor)= ### PFR (Plug Flow Reactor) **Location**: `difflow/units/pfr.py` **Class**: `PFR` **Description**: Models an ideal plug flow reactor where all fluid elements have the same residence time. No axial mixing, but perfect radial mixing is assumed. #### Process Role PFRs are preferred for: - Gas-phase reactions - High conversion requirements - Reactions where selectivity depends on conversion - Fast reactions #### Parameters ```python @dataclass class PFRParams: V: float # Total reactor volume (m³) rate_fn: Callable # Rate function: rate_fn(C, T, rate_params) -> r stoich: Array # Stoichiometry matrix [n_species × n_reactions] rate_params: dict # Parameters passed to rate_fn species_order: list # List of species names dH_rxn: Array = None # Heat of reaction (J/mol), required for non-isothermal dP_dV: float = None # Pressure gradient along the reactor (Pa/m³), optional rtol: float = 1e-6 # Relative tolerance for ODE solver atol: float = 1e-8 # Absolute tolerance for ODE solver n_save_points: int = 101 # Points to save in profile output ``` #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `inlet` | Stream | - | Inlet stream | | `T_spec` | float | K | Outlet temperature (isothermal mode) | | `volumetric_flow` | float | m³/s | Volumetric flow rate | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `outlet` | Stream | - | Outlet stream | | `info['conversion']` | dict | - | Conversion of each species | | `info['profiles']['V']` | Array | m³ | Volume along reactor | | `info['profiles']['F']` | Array | mol/s | Molar flows along reactor, `[n_save_points, n_species]` | | `info['profiles']['T']` | Array | K | Temperature along reactor | #### Governing Equations **Material Balance** (differential): $$\frac{dF_i}{dV} = \sum_j \nu_{ij} r_j$$ **Energy Balance** (adiabatic): $$\frac{dT}{dV} = \frac{-\sum_j r_j \Delta H_{rxn,j}}{F_{total} C_{p,mix}}$$ **Integration Method**: Adaptive ODE integration using diffrax (Tsit5 or Dopri5 solvers) #### Example Usage ```python from difflow.units.pfr import PFR, PFRParams from difflow import make_stream import jax.numpy as jnp # Define rate function: A -> B (first-order) def pfr_rate_fn(C, T, params): """Rate function: r = k * C_A with Arrhenius temperature dependence.""" k_ref, E_a, T_ref = params['k_ref'], params['E_a'], params['T_ref'] R = 8.314 k = k_ref * jnp.exp(E_a / R * (1/T_ref - 1/T)) return jnp.array([k * C['A']]) pfr_params = PFRParams( V=0.005, rate_fn=pfr_rate_fn, stoich=jnp.array([[-1.0], [1.0]]), # A -> B rate_params={'k_ref': 0.5, 'E_a': 60000.0, 'T_ref': 400.0}, species_order=['A', 'B'], dH_rxn=jnp.array([-50000.0]), n_save_points=201 ) pfr = PFR(pfr_params, thermo) pfr_inlet = make_stream({'A': 2.0, 'B': 0.0}, T=400.0, P=200000.0) outlet, info = pfr(pfr_inlet, volumetric_flow=0.001) # m³/s # Plot conversion profile import matplotlib.pyplot as plt prof = info['profiles'] plt.plot(prof['V'], 1 - prof['F'][:, 0] / pfr_inlet['F_A']) plt.xlabel('Volume (m³)') plt.ylabel('Conversion') ``` --- (gaspfr-gas-phase-pfr-with-pressure-drop)= ### GasPFR (Gas-Phase PFR with Pressure Drop) **Location**: `difflow/units/pfr.py` **Class**: `GasPFR` **Description**: Extended PFR model for gas-phase reactions accounting for: - Pressure drop (Ergun equation) - Variable volumetric flow due to mole change and pressure/temperature effects #### Additional Parameters ```python @dataclass class GasPFRParams(PFRParams): alpha: float # Pressure drop parameter (1/m³); fold the Ergun # terms (diameter, void fraction, particle size) into it ``` #### Additional Equations **Pressure Drop** (Ergun equation): $$\frac{dP}{dV} = -\alpha \frac{P_0}{P} \frac{T}{T_0} \frac{F_{total}}{F_{total,0}}$$ Where $\alpha$ combines Ergun parameters: $$\alpha = \frac{G}{\rho_0 g_c D_p A_c} \left[\frac{150(1-\phi)\mu}{D_p} + 1.75 G\right] \frac{(1-\phi)}{\phi^3}$$ **Variable Volumetric Flow**: $$Q = Q_0 \frac{F_{total}}{F_{total,0}} \frac{P_0}{P} \frac{T}{T_0}$$ #### Outputs Additional outputs compared to PFR: | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `info['profiles']['P']` | Array | Pa | Pressure along reactor | | `info['pressure_drop']` | float | Pa | Total pressure drop | --- (fedbatchreactor)= ### FedBatchReactor **Location**: `difflow/units/fed_batch.py` **Class**: `FedBatchReactor`, `SemiBatchReactor` **Description**: Models fed-batch (semi-batch) reactors with time-varying feed profiles. Commonly used when reactant addition rate affects selectivity or safety. #### Process Role Fed-batch reactors are used for: - Controlling exothermic reactions - Improving selectivity by maintaining low reactant concentration - Fermentation with substrate feeding - Polymerization with monomer addition #### Parameters ```python @dataclass class FedBatchParams: V0: float # Initial reactor volume (m³) rate_fn: Callable # Rate function: rate_fn(C, T, rate_params) -> r stoich: Array # Stoichiometry matrix [n_species × n_reactions] rate_params: dict # Parameters passed to rate_fn species_order: list # List of species names dH_rxn: Array = None # Heat of reaction (J/mol), None for isothermal ``` #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `inlet` | Stream | - | Feed stream composition | | `feed_flow` | Callable | mol/s | Feed rate as function of time: `F(t)` | | `T_profile` | Callable | K | Temperature as function of time: `T(t)` | #### Governing Equations **Volume Change**: $$\frac{dV}{dt} = Q_{feed}$$ **Material Balance**: $$\frac{d(V C_i)}{dt} = F_{in} C_{in,i} + V \sum_j \nu_{ij} r_j$$ Or equivalently: $$\frac{dN_i}{dt} = F_{in,i} + V \sum_j \nu_{ij} r_j$$ **Energy Balance**: $$\frac{d(V \rho C_p T)}{dt} = F_{in} \rho_{in} C_{p,in} T_{in} + V \sum_j r_j (-\Delta H_{rxn,j}) + Q$$ #### Utility Functions ```python from difflow.units.fed_batch import ( FedBatchParams, batch_time_for_conversion, optimal_feed_profile, ) fb_params = FedBatchParams( V0=1.0, rate_fn=lambda C, T, p: jnp.array([p['k'] * C['A']]), stoich=jnp.array([[-1.0], [1.0]]), # A -> B rate_params={'k': 1e-3}, species_order=['A', 'B'], ) C0 = {'A': 1000.0, 'B': 0.0} # mol/m³ # Calculate batch time for target conversion (here 95 % of A) t_batch = batch_time_for_conversion(fb_params, C0, 350.0, 'A', 0.95) # Generate an optimal piecewise-constant feed profile (maximize product B) feed_fn, t_opt = optimal_feed_profile( 'max_yield', fb_params, C0, 350.0, 'B', feed_composition={'A': 2000.0, 'B': 0.0}, V_max=2.0, t_max=3600.0, n_intervals=4, n_sim_steps=50, ) ``` --- ### SemiBatchReactor **Location**: `difflow/units/fed_batch.py` **Class**: `SemiBatchReactor` **Description**: `FedBatchReactor` under the name the process usually goes by. Same model, same parameters, same equations. A semi-batch reactor is a batch vessel into which one reactant is fed gradually --- to cap the heat release, to hold a reactant concentration low for selectivity, or to keep a hazardous intermediate from accumulating. That is exactly the fed-batch model above, so `SemiBatchReactor` is a subclass of `FedBatchReactor` that changes nothing but its display symbol: ```python from difflow.units.fed_batch import SemiBatchReactor, FedBatchParams reactor = SemiBatchReactor(fb_params, thermo, mode="isothermal") # fb_params: FedBatchParams ``` Everything in [FedBatchReactor](#fedbatchreactor) --- parameters, feed and temperature profiles, the material and energy balances, the utility functions --- applies unchanged. Pick the name that makes the flowsheet read correctly; there is no modelling difference to weigh. --- (separators)= ## Declarative Kinetics The reactors above take `rate_fn` as a Python callable. That is expressive, but it is *code*, not data — it cannot be written to a file, built from a form, or round-tripped through a GUI. `mass_action_kinetics` builds the callable from plain dictionaries instead, so a reaction network can be stored, edited and shared as data. ```python from difflow import CSTR, CSTRParams, mass_action_kinetics reactions = [{ "equation": "A -> B", "reactants": {"A": 1.0}, "products": {"B": 1.0}, "rate_params": {"A": 1.0e6, "Ea": 50_000.0, "n": 0.0}, }] kin = mass_action_kinetics(reactions, species_order=["A", "B"]) cstr = CSTR(CSTRParams(V=1.5, **kin.params_kwargs())) ``` `params_kwargs()` supplies `rate_fn`, `stoich`, `rate_params` and `species_order` — every rate-law field the reactors need. The result is numerically identical to the equivalent hand-written callable. The dictionary format is exactly what [`import_reactions`](thermodynamics.md) returns from a Cantera YAML file, so a published mechanism goes straight into a reactor: ```python from difflow import import_reactions reactions = import_reactions("mech.yaml") kin = mass_action_kinetics(reactions, reverse="forward_only") ``` ### The rate law $$k_j(T) = A_j \, T^{n_j} \exp\!\left(\frac{-E_{a,j}}{R T}\right), \qquad r_j = k_j \prod_i C_i^{\alpha_{ji}}$$ Orders $\alpha$ come from the reactant stoichiometry unless given explicitly via `orders=`, which covers empirical rate laws where the order is not the coefficient. A reversible reaction subtracts the reverse term scaled by its equilibrium constant: $$r_j = k_j \left( \prod_i C_i^{\alpha_{ji}} - \frac{1}{K_{eq,j}} \prod_i C_i^{\beta_{ji}} \right)$$ **Parameters:** - `reactions` — one dict per reaction with `reactants`, `products` and `rate_params` (`A`, `Ea`, `n`); optionally `equation`, `reversible`, `type` and `K_eq`. - `species_order` — fixes the rows of `stoich`; defaults to the sorted union of every species mentioned. - `reverse` — `"error"` (default), `"forward_only"`, or `"equilibrium"`. - `orders` — per-reaction `{species: order}` overrides. `None` keeps stoichiometric orders; `{}` means zeroth order in everything. ### What it refuses, and why Two specifications raise `KineticsSpecError` rather than being approximated, because in both cases a guess produces a plausible number that is wrong: - **Reversible reactions, by default.** Forward Arrhenius parameters alone do not determine the reverse rate. Pass `reverse="equilibrium"` with a `K_eq` on each reaction, or `reverse="forward_only"` to drop the reverse term as an explicit, recorded approximation. - **Three-body, falloff and other pressure-dependent types.** These need their own rate law; mass action covers elementary reactions only. ### Units Concentrations mol/m³, temperature K, activation energy J/mol, rates mol/m³/s. Cantera files declare their own units in a `units:` block — a mechanism written in cm³ or kcal/mol imports numerically unchanged and will be **silently wrong**, so check that block before trusting a rate constant. --- ## Separators (flash-drum)= ### Flash Drum **Location**: `difflow/units/flash.py` **Classes**: `Flash`, `EOSFlash`, `PHFlash` **Description**: Performs vapor-liquid equilibrium (VLE) separation. The feed is separated into vapor and liquid phases at equilibrium conditions. #### Process Role Flash drums are used for: - Separating light and heavy components - Pressure reduction with phase separation - Overhead condensers - Feed preparation for distillation #### Parameters ```python @dataclass class FlashParams: species_order: list[str] # List of species names for array ordering @dataclass class EOSFlashParams: species_order: list[str] # List of species names for array ordering eos_type: str = "PR" # "PR" (Peng-Robinson) or "SRK" ``` #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `inlet` | Stream | - | Feed stream | | `T` | float | K | Flash temperature (optional override) | | `P` | float | Pa | Flash pressure (optional override) | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `liquid` | Stream | - | Liquid product | | `vapor` | Stream | - | Vapor product | | `info['V_frac']` | float | - | Vapor fraction | | `info['K']` | dict | - | K-values for each species | | `info['x']` | dict | - | Liquid mole fractions | | `info['y']` | dict | - | Vapor mole fractions | #### Governing Equations **Rachford-Rice Equation**: $$f(V) = \sum_i \frac{z_i (K_i - 1)}{1 + V(K_i - 1)} = 0$$ Where: - $z_i$: Feed mole fraction of species i - $K_i$: Equilibrium ratio (K-value) = $y_i / x_i$ - $V$: Vapor fraction (moles vapor / moles feed) **Phase Compositions**: $$x_i = \frac{z_i}{1 + V(K_i - 1)}$$ $$y_i = \frac{K_i z_i}{1 + V(K_i - 1)}$$ **K-Value Calculation** (Raoult's Law for `Flash`): $$K_i = \frac{P_i^{sat}(T)}{P}$$ **K-Value Calculation** (Fugacity-based for `EOSFlash`): $$K_i = \frac{\phi_i^L}{\phi_i^V}$$ Where $\phi_i$ are fugacity coefficients from Peng-Robinson or SRK equation of state. **Material Balance**: $$F = L + V$$ $$F z_i = L x_i + V y_i$$ #### Flash Classes ##### Flash (Ideal) Uses Raoult's law K-values from IdealThermo. Suitable for ideal or near-ideal mixtures. ```python from difflow import IdealThermo, make_stream from difflow.database import get_species_data from difflow.units.flash import Flash, FlashParams thermo = IdealThermo({n: get_species_data(n) for n in ['methane', 'ethane', 'propane']}) flash = Flash(FlashParams(species_order=['methane', 'ethane', 'propane']), thermo) feed = make_stream({'methane': 0.5, 'ethane': 0.3, 'propane': 0.2}, T=200.0, P=500000.0) liquid, vapor, info = flash(feed) print(f"Vapor fraction: {info['V_frac']:.3f}") ``` (eosflash)= ##### EOSFlash (Non-Ideal) Uses fugacity coefficients from cubic equations of state (Peng-Robinson or SRK) for non-ideal VLE. ```python from difflow.units.flash import EOSFlash, EOSFlashParams from difflow.eos import PengRobinson, CriticalProperties # Define species with critical properties species_data = { "methane": CriticalProperties("methane", 190.6, 4.6e6, 0.011), "ethane": CriticalProperties("ethane", 305.4, 4.9e6, 0.099), "propane": CriticalProperties("propane", 369.8, 4.2e6, 0.152), } eos = PengRobinson(species_data) params = EOSFlashParams(species_order=["methane", "ethane", "propane"], eos_type="PR") flash = EOSFlash(params, eos) feed = make_stream({'methane': 40.0, 'ethane': 30.0, 'propane': 30.0}, T=250.0, P=2e6) liquid, vapor, info = flash(feed) ``` (phflash)= ##### PHFlash (Isenthalpic) Performs adiabatic flash at constant pressure and enthalpy. Solves for flash temperature. ```python from difflow.units.flash import PHFlash, FlashParams thermo_lh = IdealThermo({n: get_species_data(n) for n in ['n_pentane', 'n_heptane']}) ph_flash = PHFlash(FlashParams(species_order=['n_pentane', 'n_heptane']), thermo_lh) # Hot liquid feed, flash to lower pressure feed = make_stream({'n_pentane': 50.0, 'n_heptane': 50.0}, T=380.0, P=101325.0) liquid, vapor, info = ph_flash(feed, P=30000.0) print(f"Flash temperature: {info['T_flash']:.1f} K") print(f"Vapor fraction: {info['V_frac']:.3f}") ``` #### Bubble and Dew Point Methods The `Flash` class provides methods for calculating phase boundaries: ```python flash = Flash(FlashParams(species_order=['n_pentane', 'n_heptane']), thermo_lh) feed = make_stream({'n_pentane': 50.0, 'n_heptane': 50.0}, T=350.0, P=50000.0) # Pressure calculations (at specified T) P_bubble = flash.bubble_point_pressure(feed, T=350.0) # First bubble forms P_dew = flash.dew_point_pressure(feed, T=350.0) # Last drop condenses # Temperature calculations (at specified P) T_bubble = flash.bubble_point_temperature(feed, P=50000.0) T_dew = flash.dew_point_temperature(feed, P=50000.0) print(f"Bubble point: T={T_bubble:.1f} K, P={P_bubble:.0f} Pa") print(f"Dew point: T={T_dew:.1f} K, P={P_dew:.0f} Pa") ``` --- ### Mixer **Location**: `difflow/units/flash.py` **Class**: `Mixer` **Description**: Combines multiple inlet streams into a single outlet stream by summing molar flows. #### Governing Equations **Mass Balance**: $$F_{out,i} = \sum_k F_{k,i}$$ **Energy Balance** (adiabatic mixing): $$T_{out} = \frac{\sum_k F_k C_{p,k} T_k}{\sum_k F_k C_{p,k}}$$ (Simplified for ideal mixing with similar heat capacities) #### Example Usage ```python from difflow.units.flash import Mixer mixer = Mixer(['A', 'B', 'C']) # optional: thermo=..., phase='liquid' stream1 = make_stream({'A': 1.0, 'B': 0.5, 'C': 0.0}, T=350.0, P=101325.0) stream2 = make_stream({'A': 0.0, 'B': 0.3, 'C': 0.2}, T=360.0, P=101325.0) # Streams are passed as separate arguments; each carries every species outlet, info = mixer(stream1, stream2) ``` --- ### Splitter **Location**: `difflow/units/flash.py` **Class**: `Splitter` **Description**: Divides a single inlet stream into multiple outlet streams based on split fractions. #### Governing Equations $$F_{out,k} = \alpha_k F_{in}$$ Where $\sum_k \alpha_k = 1$ All outlet streams have the same composition and temperature as the inlet. #### Example Usage ```python from difflow.units.flash import Splitter splitter = Splitter(species_order=['A', 'B']) inlet = make_stream({'A': 1.0, 'B': 0.5}, T=350.0, P=101325.0) # Split into 3 streams: 50%, 30%, 20% out1, out2, out3, info = splitter(inlet, split_frac=[0.5, 0.3, 0.2]) ``` --- ## Distillation ### ShortcutColumn **Location**: `difflow/units/distillation.py` **Class**: `ShortcutColumn` **Description**: Uses shortcut methods (Fenske-Underwood-Gilliland) for rapid distillation column design and rating calculations. #### Process Role Shortcut methods are used for: - Initial column design estimates - Optimization studies - Screening alternatives - Quick sensitivity analysis #### Parameters ```python @dataclass class ShortcutColumnParams: species_order: list[str] # Species names; sets array ordering light_key: str # Name of the light key component heavy_key: str # Name of the heavy key component x_D_LK: float = 0.99 # Fractional recovery of LK in the distillate x_B_HK: float = 0.99 # Fractional recovery of HK in the bottoms ``` `x_D_LK` and `x_B_HK` are **recoveries**, not mole fractions, despite the names: `x_D_LK=0.99` sends 99 % of the feed's light key overhead. The distillate's actual light-key mole fraction comes back in `info["x_D"]`. The reflux ratio is not a parameter — it is an argument to the call, along with the column pressure and the feed quality: ```python distillate, bottoms, info = column(feed, R=3.0, P=101325.0, q=1.0) ``` #### Governing Equations **Relative Volatility**: $$\alpha_{ij} = \frac{K_i}{K_j}$$ The second equality usually written here, $\alpha_{ij} = P_i^{sat}/P_j^{sat}$, holds only under Raoult's law. On a `CubicThermo` the K-values are fugacity coefficient ratios and the vapor pressures cancel out of nothing — which is the whole reason the volatilities differ between the two packages. **Average Relative Volatility** (geometric mean): $$\bar{\alpha} = (\alpha_{top} \cdot \alpha_{bottom})^{0.5}$$ The two ends are the column's actual ends, not estimates around the feed: the top is a total condenser, so $T_{top}$ is the **bubble point of the distillate**, and $T_{bot}$ is the **bubble point of the bottoms** in the reboiler. Those are the temperatures the product streams come out at, reported as `info["T_condenser"]` and `info["T_reboiler"]`. They also make the design a small fixed point, since $\bar\alpha$ sets the product split through Hengstebeck-Geddes and the split sets the bubble points in turn. The column sweeps it a few times from the feed's own bubble point; it settles to under a hundredth of a degree by the third sweep, and the whole loop is unrolled, so `jax.grad` runs through it. **Non-key distribution (Hengstebeck-Geddes).** The keys go where their recoveries put them; every other species follows the straight line in $(\log\alpha, \log d/b)$ through both keys: $$\log\frac{d_i}{b_i} = A + C\log\alpha_i,\qquad A = \log\left(\frac{d}{b}\right)_{HK},\qquad C = \frac{\log(d/b)_{LK} - \log(d/b)_{HK}}{\log\alpha_{LK}}$$ with $\alpha$ relative to the heavy key, so $\alpha_{HK} = 1$ (Geddes, *AIChE J.* 4, 389 (1958); Hengstebeck, *Distillation*, Reinhold (1961); page and equation numbers unverified). The distillate share $d_i/(d_i + b_i)$ is the logistic function of that line, so each species' balance closes exactly. Before this was corrected the constants were $A = \log(d/b)_{LK} - \log(d/b)_{HK}$ and $C = \log(d/b)_{LK}/\log\alpha_{LK}$, a line that misses the heavy key: on a propane/isobutane depropanizer it sent 99.9 % of the n-butane overhead. `tests/test_distillation.py::TestShortcutColumnNonKeyDistribution` pins the line through both keys. Two consequences worth knowing, because the old estimate had neither: the end temperatures no longer move when you feed the same mixture in hotter, and they do move with column pressure. **Fenske Equation** (minimum stages): $$N_{min} = \frac{\ln\left[\frac{x_{D,LK}}{x_{B,LK}} \cdot \frac{x_{B,HK}}{x_{D,HK}}\right]}{\ln \bar{\alpha}_{LK/HK}}$$ **Underwood Equations** (minimum reflux): For each component i: $$\sum_i \frac{\alpha_i x_{F,i}}{\alpha_i - \theta} = 1 - q$$ $$R_{min} + 1 = \sum_i \frac{\alpha_i x_{D,i}}{\alpha_i - \theta}$$ Where: - $\theta$: Root between $\alpha_{HK}$ and $\alpha_{LK}$ - $q$: Feed quality (1 for saturated liquid, 0 for saturated vapor) **Gilliland Correlation** (actual stages): $$\frac{N - N_{min}}{N + 1} = 1 - \exp\left[\frac{(1 + 54.4X)(X - 1)}{(11 + 117.2X)(X^{0.5})}\right]$$ Where: $$X = \frac{R - R_{min}}{R + 1}$$ **Feed Stage Location** (Kirkbride correlation): $$\frac{N_R}{N_S} = \left[\frac{B}{D} \cdot \frac{x_{F,HK}}{x_{F,LK}} \cdot \left(\frac{x_{B,LK}}{x_{D,HK}}\right)^2\right]^{0.206}$$ #### Outputs | Key | Type | Units | Description | |-----|------|-------|-------------| | `distillate` | Stream | - | Overhead product, at the condenser temperature | | `bottoms` | Stream | - | Bottom product, at the reboiler temperature | | `info['N_min']` | float | - | Minimum stages (Fenske) | | `info['R_min']` | float | - | Minimum reflux ratio (Underwood) | | `info['N']` | float | - | Actual stages (Gilliland) | | `info['N_feed']` | float | - | Feed stage (Kirkbride) | | `info['D']`, `info['B']` | float | mol/s | Distillate and bottoms flow | | `info['x_D']`, `info['x_B']` | dict | - | Product compositions by species | | `info['T_condenser']` / `info['T_top']` | float | K | Condenser temperature = bubble point of $x_D$ | | `info['T_reboiler']` / `info['T_bot']` | float | K | Reboiler temperature = bubble point of $x_B$ | | `info['Q_condenser']` | float | W | Condenser duty (negative: heat removed) | | `info['Q_reboiler']` | float | W | Reboiler duty (positive: heat added) | | `info['alpha']` | dict | - | Relative volatilities vs the heavy key | | `info['alpha_LK']` | float | - | Light key's relative volatility | | `info['alpha_top']`, `info['alpha_bot']` | dict | - | Volatilities at each column end | | `info['alpha_variation']` | dict | - | Relative spread between the two ends | | `info['alpha_varies_significantly']` | bool | - | True if that spread exceeds 0.3 | | `info['theta']` | float | - | Underwood root | | `info['close_boiling']` | bool | - | True if $\bar\alpha \approx 1$ capped $N_{min}$ | | `info['near_min_reflux']` | bool | - | True if $R \approx R_{min}$ | | `info['negative_flows_detected']` | bool | - | True if the split produced a negative flow | | `info['feasible']` | bool | - | All three checks above passed | #### Example Usage ```python from difflow import IdealThermo, make_stream from difflow.database import get_species_data from difflow.units.distillation import ShortcutColumn, ShortcutColumnParams names = ['benzene', 'toluene', 'ethylbenzene'] thermo = IdealThermo({s: get_species_data(s) for s in names}) params = ShortcutColumnParams( species_order=names, light_key='benzene', heavy_key='toluene', x_D_LK=0.99, # 99 % of the benzene recovered overhead x_B_HK=0.99, # 99 % of the toluene recovered in the bottoms ) column = ShortcutColumn(params, thermo) feed = make_stream({'benzene': 40.0, 'toluene': 35.0, 'ethylbenzene': 25.0}, T=370.0, P=101325.0) distillate, bottoms, info = column(feed, R=3.0, P=101325.0, q=1.0) print(f"Minimum stages: {float(info['N_min']):.1f}") print(f"Minimum reflux: {float(info['R_min']):.2f}") print(f"Actual stages: {float(info['N']):.1f}") print(f"Feed stage: {float(info['N_feed']):.1f}") print(f"Condenser: {float(info['T_condenser']):.1f} K, " f"{float(info['Q_condenser'])/1e6:.2f} MW") print(f"Reboiler: {float(info['T_reboiler']):.1f} K, " f"{float(info['Q_reboiler'])/1e6:.2f} MW") ``` #### Utility Functions The design correlations are also available on their own, taking plain numbers rather than a column object. Note the argument order — each takes the quantities in the order the correlation is written, not the order the column computes them: ```python from difflow.units.distillation import ( fenske_stages, # (x_D_LK, x_B_LK, alpha) minimum_reflux_ratio, # (z_LK, z_HK, x_D_LK, alpha, q=1.0) gilliland_stages, # (R, R_min, N_min) column_diameter, # (V, rho_V, rho_L, sigma=0.02, tray_spacing=0.6) ) N_min = fenske_stages(x_D_LK=0.98, x_B_LK=0.02, alpha=2.4) R_min = minimum_reflux_ratio(z_LK=0.45, z_HK=0.55, x_D_LK=0.98, alpha=2.4, q=1.0) N = gilliland_stages(R=1.3 * R_min, R_min=R_min, N_min=N_min) diameter = column_diameter(V=50.0, rho_V=3.0, rho_L=800.0) # mol/s, kg/m^3 ``` Relative volatility is a method on the column, not a free function (`column.relative_volatility(T, P, x=None)`), because it needs the key components from the column's parameters — and, on a `CubicThermo`, the phase composition. --- ### DistillationColumn (Rigorous) **Location**: `difflow/units/distillation.py` **Class**: `DistillationColumn` **Description**: Stage-by-stage calculation using MESH equations (Material balance, Equilibrium, Summation, enthalpy balance). #### Parameters ```python @dataclass class DistillationColumnParams: species_order: list[str] # Species names; sets array ordering n_stages: int # Equilibrium stages: the trays plus the reboiler feed_stage: int # Zero-based stage index from the bottom condenser_type: str = 'total' # Only 'total' is implemented P: float = 101325.0 # One pressure for the whole column q: float = 1.0 # Feed quality (1 = saturated liquid) ``` **Stage numbering is zero-based from the bottom**: stage 0 is the reboiler, stage `n_stages - 1` is the top tray, and every profile in `info` is in that order. The condenser is *not* a stage — it sits outside the cascade — so `n_stages` counts the trays plus the reboiler, and `feed_stage=10` in a 20-stage column puts the feed halfway up. Stages below the feed stage are the stripping section and stages above it the rectifying section. There is one column pressure: no tray pressure drop. The feed stage itself is split by the feed, which contributes `q F` to the liquid running down and `(1 - q) F` to the vapour running up. Its liquid is therefore a **stripping**-section flow and its vapour a **rectifying**-section one: $$L_j = \begin{cases} L' = L + qF & j \le n_F \\ L = R D & j > n_F \end{cases} \qquad V_j = \begin{cases} V' = V - (1-q)F & j < n_F \\ V = (R+1) D & j \ge n_F \end{cases}$$ That is one convention, not two: both flows follow from asking whether the horizontal cut a stream crosses has the feed above it. `_is_rectifying_cut` and `_cmo_section_flows` in `difflow/units/distillation.py` are the single definition, and `_cmo_flows` — the only thing that builds an L/V profile — is a thin wrapper over them, so both solver paths read the same boundary. `q` sets those rates on both paths, and it is also the feed's thermal condition in the energy balance: $$h_F = q\, h^L(z, T_F) + (1 - q)\, H^V(z, T_F)$$ both phase enthalpies at the feed stream's own temperature. So a saturated-vapour feed arrives with its latent heat already in it and the reboiler is not charged for it: going from `q = 1` to `q = 0` drops `Q_reboiler` by `F (H^V - h^L)` at the feed temperature (2.4 MW on a 100 mol/s equimolar benzene/toluene feed at 380 K), and moves the ~`F` step in the converged profile from `L` to `V` across the feed stage. The shortcut column forms its feed enthalpy the same way, from the `q` passed to the call. `q` outside `[0, 1]` is allowed, and means what a textbook means by it: `q > 1` is a subcooled feed, `q < 0` a superheated one. It has to be allowed, because on the CMO path (`use_mesh=False`) `q` is the *only* place either can be said — `_cmo_section_rates` is the whole model there, and `T_feed` never reaches it. `L_strip = L_rect + q F` with `q > 1` is exactly how the extra internal reflux of a subcooled feed is written. Underwood's equation takes it as written too. The one place it is clamped is the **feed enthalpy**, and that clamp is the physics rather than a guard. Both phase enthalpies are evaluated at the feed's own temperature, so at `q = 1.3` the honest answer is `h_liquid(z, T_feed)`: an all-liquid feed below its bubble point, with the subcooling carried by `T_feed`. Forming `1.3 h^L - 0.3 H^V` would subtract three tenths of a latent heat that is not there, and count the departure from saturation twice. `condenser_type='partial'` raises `NotImplementedError` rather than being silently solved as a total condenser. The reflux ratio and the product split are arguments to the call, not parameters — give exactly one of `D_spec` or `B_spec`: ```python distillate, bottoms, info = column(feed, R=2.0, B_spec=40.0) ``` #### Governing Equations (MESH) For each stage j: **Material Balance**: $$L_{j-1} x_{i,j-1} + V_{j+1} y_{i,j+1} + F_j z_{i,j} = L_j x_{i,j} + V_j y_{i,j}$$ **Equilibrium**: $$y_{i,j} = K_{i,j} x_{i,j}$$ **Summation**: $$\sum_i x_{i,j} = 1$$ $$\sum_i y_{i,j} = 1$$ **Enthalpy Balance**: $$L_{j-1} H^L_{j-1} + V_{j+1} H^V_{j+1} + F_j H^F_j = L_j H^L_j + V_j H^V_j + Q_j$$ #### Product temperatures The bottoms leaves the reboiler, which is a stage, so it is at `info["T_profile"][0]`. The distillate does not leave a stage — it leaves the condenser. A total condenser condenses the whole of the top stage's vapor, so the distillate (and the reflux returned with it) is a **saturated liquid of composition $x_D$ at its own bubble point**, reported as `info["T_condenser"]`. That is not the top stage temperature. The top stage sits at the bubble point of its liquid $x_{top}$, equivalently the dew point of the vapor $y_{top} = x_D$ it sends up, and a mixture's dew point is above its bubble point. The gap is the boiling range of the distillate itself: | distillate | top stage $T$ | condenser $T$ | gap | |---|---|---|---| | 99.9999 % benzene / toluene, 1 atm | 368.665 K | 368.665 K | 0.0001 K | | C3-C8 cut, 10 bar (Peng-Robinson) | 400.6 K | 362.9 K | 38 K | A one-component distillate has no boiling range and so no gap, which is why the binary row reads as zero: at $\alpha \approx 8$ over 15 stages that column takes 27 µmol/s of toluene overhead and nothing more. The gap is a property of the cut, not of the column. The same distinction runs through the energy balance: $Q_{cond}$ takes the top stage vapor down to that condensed state, and the reflux re-enters the top stage subcooled, at the condenser temperature rather than the tray's. #### Thermodynamics: ideal K-values or a cubic EOS The column takes either an `IdealThermo` or a [`CubicThermo`](thermodynamics.md), and the choice is the whole of the difference between a near-ideal separation and a hydrocarbon one: | | `IdealThermo` | `CubicThermo` | |---|---|---| | $K_i$ | $P^{sat}_i(T)/P$ (Raoult) | $\hat\phi^L_i(T,P,x)\,/\,\hat\phi^V_i(T,P,y)$ (PR or SRK) | | stage enthalpy | ideal-gas $C_p$ + Watson $H_{vap}$ | ideal-gas $C_p$ + EOS departure | | $K_i$ depends on composition | no | yes | ```python from difflow import IdealThermo, CubicThermo, PengRobinson, make_stream from difflow.database import get_critical_props, get_species_data from difflow.units.distillation import DistillationColumn, DistillationColumnParams names = ["propane", "isobutane", "n_butane", "isopentane", "n_pentane", "n_hexane", "n_heptane", "n_octane"] ideal = IdealThermo({s: get_species_data(s) for s in names}) eos = PengRobinson({s: get_critical_props(s) for s in names}) column = DistillationColumn( DistillationColumnParams(species_order=names, n_stages=20, feed_stage=10, condenser_type="total", P=10e5), thermo=CubicThermo(ideal, eos), # Peng-Robinson K-values and enthalpies ) feed = make_stream( dict(zip(names, [10.0, 7.0, 7.0, 8.0, 8.0, 20.0, 10.0, 30.0])), T=380.0, P=10e5, ) distillate, bottoms, info = column(feed, R=2.0, B_spec=40.0) ``` Use the EOS for light hydrocarbons at pressure. Raoult's law fails there in a one-sided way: a C3-C8 cut at 10 bar puts its heavy end within a degree of the EOS answer and its light end tens of degrees off, so an ideal-K column looks plausible at the reboiler and is wrong at the condenser. Two things about the EOS path are worth knowing: - **The K-values depend on composition**, so they are a fixed point rather than a formula. `CubicThermo.K_values_array(T, P, x)` runs a short successive substitution on $y$ internally to close it at the stage's own composition; pass both `x` and `y` if you already have a consistent pair. - **They only exist where the cubic has two roots.** Away from the bubble point -- a subcooled liquid, a superheated vapor -- there is one root, both phases take it, and $K_i$ comes back identically 1. That is the EOS reporting a single phase, but it makes $\sum_i K_i x_i - 1$ flat, so the column solves each stage's bubble point in two passes: a Newton solve on ideal K-values to land inside the two-root window, then a step-capped Newton on the EOS K-values that bisects back if a step leaves it. Nothing about this is visible in the API, but it is why the bubble point is not one `optimistix` call. The ideal pass is itself in two parts, for a reason that has nothing to do with the EOS. Vapor pressure is exponential in $-1/T$, so a couple of hundred degrees below the bubble point both $\sum_i K_i x_i$ and its slope are round-off away from zero, and a Newton step there divides one tiny number by another and lands tens of thousands of degrees away. So the solve first takes a few damped steps on $\log \sum_i K_i x_i$ -- nearly linear in $1/T$, and well scaled over the whole range -- and only then runs Newton on the residual itself. This is why a column can be handed a feed far below its own boiling point and still find its ends. (solver-paths)= #### Solver Paths Both paths run the Wang-Henke bubble-point iteration, which solves the component material balances for the whole column as a tridiagonal system. They differ only in where the L/V profiles come from: | `use_mesh` | L/V profiles | Cost | |---|---|---| | `True` (default) | corrected each iteration by the stage enthalpy balances | ~30 % more | | `False` | frozen at their constant-molar-overflow values, $L' = L + qF$ and $V' = V - (1-q)F$ | cheaper | The CMO path is a shortcut in the *energy* balance, not in the material balance: because the tridiagonal solve is the component balance, summing it over the stages telescopes to $D x_{D,i} + B x_{B,i} = F z_i$. Both paths satisfy that only to within their iteration count, so both report the residual: ```python from difflow import IdealThermo, make_stream from difflow.database import get_species_data from difflow.units.distillation import DistillationColumn, DistillationColumnParams names = ['n_pentane', 'n_hexane', 'n_heptane'] column = DistillationColumn( DistillationColumnParams(species_order=names, n_stages=20, feed_stage=10, P=101325.0), thermo=IdealThermo({s: get_species_data(s) for s in names}), ) feed = make_stream({'n_pentane': 30.0, 'n_hexane': 40.0, 'n_heptane': 30.0}, T=360.0, P=101325.0) distillate, bottoms, info = column(feed, R=2.0, B_spec=40.0, use_mesh=False) print(info['balance_error']) # (n_species,) D_i + B_i - F_i, mol/s print(info['balance_error_rel']) # max |error| / F_total ``` If `balance_error_rel` is larger than your problem tolerates, raise `cmo_iter` (default 30). The residual falls geometrically with it, but at a rate the column sets, so treat the reported number as the answer rather than assuming a count is enough: | case | `cmo_iter` 30 | 60 | 100 | |---|---|---|---| | 20-stage ternary, $R = 2$ | 1.6e-4 | 5.0e-10 | — | | 12-stage binary, $R = 1.2$ (near $R_{min}$) | 1.9e-4 | 1.8e-4 | 1.7e-4 | Near minimum reflux the fixed-point iteration converges with a rate close to one, and more sweeps buy almost nothing. That is a property of the column, not a defect in the solver — but it is exactly why the residual is reported instead of asserted. #### Outputs | Key | Type | Units | Description | |-----|------|-------|-------------| | `distillate` | Stream | - | Overhead product, at the condenser temperature | | `bottoms` | Stream | - | Bottom product, at the reboiler temperature | | `info['T_profile']` | (n,) array | K | Stage temperatures, reboiler first | | `info['x_profile']` | (n, nc) array | - | Liquid compositions per stage | | `info['y_profile']` | (n, nc) array | - | Vapor compositions per stage | | `info['T_condenser']` | float | K | Condenser temperature = distillate T | | `info['T_reboiler']` | float | K | Reboiler temperature = bottoms T = `T_profile[0]` | | `info['D']`, `info['B']` | float | mol/s | Distillate and bottoms flow | | `info['Q_condenser']` | float | W | Condenser duty (negative: heat removed) | | `info['Q_reboiler']` | float | W | Reboiler duty (positive: heat added) | | `info['L_profile']`, `info['V_profile']` | (n,) array | mol/s | Internal flows — `use_mesh=True` only | | `info['L_rect']`, `info['V_rect']` | float | mol/s | Rectifying flows — `use_mesh=False` only | | `info['balance_error']` | (nc,) array | mol/s | Component balance residual $D_i + B_i - F_i$ | | `info['balance_error_rel']` | float | - | `max abs(balance_error) / F_total` | #### Example Usage ```python from difflow import IdealThermo, make_stream from difflow.database import get_species_data from difflow.units.distillation import DistillationColumn, DistillationColumnParams names = ['n_pentane', 'n_hexane', 'n_heptane'] thermo = IdealThermo({s: get_species_data(s) for s in names}) column = DistillationColumn( DistillationColumnParams( species_order=names, n_stages=20, # 19 trays plus the reboiler feed_stage=10, # halfway up, counting from the reboiler at 0 P=101325.0, ), thermo=thermo, ) feed = make_stream({'n_pentane': 30.0, 'n_hexane': 40.0, 'n_heptane': 30.0}, T=360.0, P=101325.0) distillate, bottoms, info = column(feed, R=2.0, B_spec=40.0) print(f"Distillate: {float(distillate['T']):.1f} K, " f"{float(distillate['F_n_pentane']):.1f} mol/s n-pentane") print(f"Bottoms: {float(bottoms['T']):.1f} K, " f"{float(bottoms['F_n_heptane']):.1f} mol/s n-heptane") print(f"Duties: {float(info['Q_condenser'])/1e6:.2f} / " f"{float(info['Q_reboiler'])/1e6:.2f} MW") ``` Pass `use_mesh=False` for the faster constant-molar-overflow solution, which skips the energy-balance correction to the internal flows but solves the same component material balances — see [Solver Paths](#solver-paths) above. --- ## Heat Exchangers ### Heater **Location**: `difflow/units/heat_exchanger.py` **Class**: `Heater` **Description**: Single-stream heater that increases stream temperature using an external heat source (steam, hot oil, electric). #### Parameters ```python @dataclass class HeaterParams: duty: float = None # Heat duty (W); the mode is set by which of T_out: float = None # Outlet temperature (K); these three you give UA: float = None # Overall HTC x Area (W/K), with T_utility T_utility: float = None # Utility temperature (K) Cp: float = None # Constant heat capacity (J/mol/K) phase: str = None # Force 'liquid'/'vapor' for the thermo enthalpy ``` The operating mode is implied by which parameter is set -- `duty`, `T_out`, or `UA` together with `T_utility` -- there is no `mode` field. #### Energy models The heater has two, and the choice matters more than any other parameter: **Constant Cp** -- `Heater(HeaterParams(T_out=400.0, Cp=75.0))`: $$Q = F_{total} C_p (T_{out} - T_{in})$$ Sensible heat only. A constant $C_p$ cannot carry latent heat, so this is wrong -- often by a factor of several -- for any stream that vaporizes or condenses across the unit. **Thermo** -- `Heater(HeaterParams(T_out=400.0), thermo=thermo)`: $$Q = H(T_{out}, P) - H(T_{in}, P)$$ with $H$ from the thermo's stream enthalpy. A [`CubicThermo`](thermodynamics.md) supplies a two-phase flash enthalpy, so the duty carries the real temperature dependence of the heat capacity *and* the latent heat of any phase change. Use this whenever the stream may change phase, and whenever you are comparing against a rigorous simulator. Pass `phase='liquid'` or `phase='vapor'` to force a single-phase enthalpy instead (which is what an `IdealThermo` provides). With neither `Cp` nor `thermo`, the unit falls back to `DEFAULT_CP` (75 J/mol/K, roughly liquid water) and raises a `DefaultCpWarning`. Turn that into an error to make the fallback fatal: ```python import warnings from difflow import DefaultCpWarning warnings.simplefilter("error", DefaultCpWarning) ``` **LMTD Rating** (for utility heating): $$Q = UA \cdot LMTD$$ $$LMTD = \frac{(T_U - T_{in}) - (T_U - T_{out})}{\ln\left(\frac{T_U - T_{in}}{T_U - T_{out}}\right)}$$ Where $T_U$ is the utility (steam) temperature. On the constant-Cp path this is solved in closed form by effectiveness-NTU with an infinite-capacity utility; with a thermo it is a damped fixed point on $Q$, with $T_{out}$ from inverting the enthalpy. #### Example Usage ```python from difflow import IdealThermo, make_stream from difflow.database import get_species_data from difflow.units.heat_exchanger import Heater, HeaterParams thermo = IdealThermo({'water': get_species_data('water')}) inlet = make_stream({'water': 10.0}, T=300.0, P=101325.0) # Specified duty heater = Heater(HeaterParams(duty=50000.0, Cp=75.0)) outlet, info = heater(inlet) # Specified outlet temperature, duty from the thermo (carries latent heat) heater = Heater(HeaterParams(T_out=400.0), thermo=thermo) outlet, info = heater(inlet) info["Q"] # W # Rating against a steam utility heater = Heater(HeaterParams(UA=5000.0, T_utility=450.0), thermo=thermo) outlet, info = heater(inlet) info["LMTD"], info["UA_required"] ``` --- ### Cooler **Location**: `difflow/units/heat_exchanger.py` **Class**: `Cooler` **Description**: Single-stream cooler that decreases stream temperature using cooling water or refrigeration. #### Parameters Same as Heater with appropriate utility temperatures. #### Governing Equations Same as Heater, including both energy models and the `DefaultCpWarning` fallback, with the duty sign reversed: `Q > 0` means heat removed. ```python from difflow.units.heat_exchanger import Cooler, CoolerParams cooler = Cooler(CoolerParams(T_out=280.0), thermo=thermo) outlet, info = cooler(inlet) ``` --- ### CounterCurrentHX **Location**: `difflow/units/heat_exchanger.py` **Class**: `CounterCurrentHX` **Description**: Two-stream heat exchanger with counter-current flow arrangement. Provides maximum temperature driving force. #### Process Role Counter-current heat exchangers are preferred for: - Maximum heat recovery - Heating/cooling to approach inlet temperature of other stream - Most efficient use of heat transfer area #### Parameters ```python @dataclass class HeatExchangerParams: UA: float = None # Overall HTC × Area (W/K) Cp_hot: float = None # Hot side heat capacity (J/mol·K), default 75 Cp_cold: float = None # Cold side heat capacity (J/mol·K), default 75 min_approach: float = 10.0 # Minimum temperature approach (K) ``` `CounterCurrentHX(params)` takes only the parameters: there is no thermo object and no species list, and `UA` can be overridden per call. #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_stream` | Stream | - | Hot fluid inlet | | `cold_stream` | Stream | - | Cold fluid inlet | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_outlet` | Stream | - | Hot fluid outlet | | `cold_outlet` | Stream | - | Cold fluid outlet | | `info['Q']` | float | W | Heat duty transferred | | `info['LMTD']` | float | K | Log mean temperature difference | | `info['UA']` | float | W/K | UA used for the rating | #### Governing Equations **Energy Balance**: $$Q = \dot{m}_h C_{p,h} (T_{h,in} - T_{h,out}) = \dot{m}_c C_{p,c} (T_{c,out} - T_{c,in})$$ **LMTD** (Counter-current): $$LMTD = \frac{\Delta T_1 - \Delta T_2}{\ln(\Delta T_1 / \Delta T_2)}$$ Where: - $\Delta T_1 = T_{h,in} - T_{c,out}$ - $\Delta T_2 = T_{h,out} - T_{c,in}$ **Heat Transfer Rate**: $$Q = UA \cdot LMTD$$ **Effectiveness-NTU Method**: $$\epsilon = \frac{Q}{Q_{max}} = \frac{Q}{C_{min}(T_{h,in} - T_{c,in})}$$ $$\epsilon = \frac{1 - \exp[-NTU(1 - C_r)]}{1 - C_r \exp[-NTU(1 - C_r)]}$$ Where: - $C_r = C_{min}/C_{max}$ - $NTU = UA/C_{min}$ - $C = \dot{m} C_p$ (heat capacity rate) #### Example Usage ```python from difflow.units.heat_exchanger import CounterCurrentHX, HeatExchangerParams hx = CounterCurrentHX(HeatExchangerParams(UA=50.0, Cp_hot=75.0, Cp_cold=75.0)) hot_in = make_stream({'A': 1.0}, T=450.0, P=101325.0) cold_in = make_stream({'B': 0.8}, T=300.0, P=101325.0) hot_out, cold_out, info = hx(hot_in, cold_in) print(f"Heat duty: {info['Q']/1000:.2f} kW") print(f"LMTD: {info['LMTD']:.2f} K") ``` --- ### CoCurrentHX **Location**: `difflow/units/heat_exchanger.py` **Class**: `CoCurrentHX` **Description**: Two-stream heat exchanger with co-current (parallel) flow arrangement. #### Governing Equations **LMTD** (Co-current): $$LMTD = \frac{\Delta T_1 - \Delta T_2}{\ln(\Delta T_1 / \Delta T_2)}$$ Where: - $\Delta T_1 = T_{h,in} - T_{c,in}$ - $\Delta T_2 = T_{h,out} - T_{c,out}$ **Effectiveness** (Co-current): $$\epsilon = \frac{1 - \exp[-NTU(1 + C_r)]}{1 + C_r}$$ Note: Co-current flow cannot achieve temperature cross ($T_{c,out} > T_{h,out}$). --- ### CrossFlowHX **Location**: `difflow/units/heat_exchanger.py` **Class**: `CrossFlowHX` **Description**: Two-stream heat exchanger with cross-flow arrangement where fluids flow perpendicular to each other. Effectiveness depends on mixing configuration. #### Process Role Cross-flow heat exchangers are used for: - Air-to-liquid heat transfer (HVAC systems) - Car radiators and automotive cooling - Finned-tube heat exchangers - Applications where cross-flow geometry is advantageous #### Parameters ```python @dataclass class HeatExchangerParams: UA: float = None # Overall HTC × Area (W/K) for rating Cp_hot: float = None # Hot side heat capacity (J/mol·K) Cp_cold: float = None # Cold side heat capacity (J/mol·K) min_approach: float = 10.0 # Minimum approach temperature (K) # CrossFlowHX constructor CrossFlowHX(params: HeatExchangerParams, mixing: str = "both_unmixed") ``` #### Mixing Configurations The `mixing` parameter specifies the flow arrangement: | Configuration | Description | Common Applications | |--------------|-------------|---------------------| | `both_unmixed` | Both fluids flow through separate channels (default) | Car radiators, finned-tube HX | | `cmax_mixed` | Larger heat capacity stream is mixed | Shell-and-tube with mixing in shell | | `cmin_mixed` | Smaller heat capacity stream is mixed | Special geometries | | `both_mixed` | Both fluids can mix in flow direction | Compact heat exchangers | #### Governing Equations **Energy Balance** (same as other HX types): $$Q = \dot{m}_h C_{p,h} (T_{h,in} - T_{h,out}) = \dot{m}_c C_{p,c} (T_{c,out} - T_{c,in})$$ **Effectiveness** (both unmixed): $$\epsilon = 1 - \exp\left[\frac{NTU^{0.22}}{C_r}\left(\exp(-C_r \cdot NTU^{0.78}) - 1\right)\right]$$ **Effectiveness** (Cmax mixed, Cmin unmixed): $$\epsilon = \frac{1}{C_r}\left[1 - \exp\left(-C_r(1 - \exp(-NTU))\right)\right]$$ **Effectiveness** (Cmin mixed, Cmax unmixed): $$\epsilon = 1 - \exp\left[-\frac{1}{C_r}(1 - \exp(-C_r \cdot NTU))\right]$$ **Effectiveness** (both mixed): $$\frac{1}{\epsilon} = \frac{1}{1-\exp(-NTU)} + \frac{C_r}{1-\exp(-C_r \cdot NTU)} - \frac{1}{NTU}$$ Where: - $NTU = UA/C_{min}$ (Number of Transfer Units) - $C_r = C_{min}/C_{max}$ (Heat capacity ratio) **Physical Constraint**: Cross-flow effectiveness is capped at the counter-current value to maintain physical consistency, as cross-flow should never exceed counter-current performance. #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_stream` | Stream | - | Hot fluid inlet | | `cold_stream` | Stream | - | Cold fluid inlet | | `UA` | float | W/K | Override UA value | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_outlet` | Stream | - | Hot fluid outlet | | `cold_outlet` | Stream | - | Cold fluid outlet | | `info['Q']` | float | W | Heat duty transferred | | `info['effectiveness']` | float | - | Heat exchanger effectiveness | | `info['NTU']` | float | - | Number of transfer units | | `info['LMTD']` | float | K | Log mean temperature difference | | `info['mixing']` | str | - | Mixing configuration used | | `info['approach']` | float | K | Minimum temperature approach | #### Example Usage ```python from difflow import CrossFlowHX, HeatExchangerParams, make_stream # Car radiator example (both unmixed - most common) hx = CrossFlowHX( HeatExchangerParams(UA=2000.0, Cp_hot=75.0, Cp_cold=30.0), mixing="both_unmixed" # Default ) hot_coolant = make_stream({"ethylene_glycol": 10.0}, T=368.0, P=101325.0) # 95°C cold_air = make_stream({"air": 50.0}, T=298.0, P=101325.0) # 25°C hot_out, cold_out, info = hx(hot_coolant, cold_air) print(f"Heat rejected: {info['Q']/1000:.2f} kW") print(f"Effectiveness: {info['effectiveness']:.3f}") print(f"Air outlet temp: {cold_out['T']:.1f} K") # Compare mixing configurations for config in ["both_unmixed", "cmax_mixed", "cmin_mixed", "both_mixed"]: hx = CrossFlowHX(HeatExchangerParams(UA=2000.0), mixing=config) _, _, info = hx(hot_coolant, cold_air) print(f"{config:15s}: ε = {info['effectiveness']:.4f}, Q = {info['Q']/1000:.2f} kW") ``` #### Performance Comparison For the same UA and inlet conditions, effectiveness ranking: 1. Counter-current (highest) 2. Cross-flow with mixed streams 3. Cross-flow (both unmixed) 4. Co-current (lowest) Cross-flow heat exchangers offer intermediate performance between counter-current (most efficient) and co-current (simplest), making them practical for applications where perpendicular flow geometry is advantageous. --- ### ShellAndTubeHX **Location**: `difflow/units/heat_exchanger.py` **Class**: `ShellAndTubeHX` **Description**: Multi-pass shell-and-tube exchanger: counter-current duty reduced by the LMTD correction factor $F$. #### Process Role A 1-2N TEMA exchanger is not counter-current. Each tube pass runs with the shell fluid and the next against it, so part of the area works against a smaller driving force than the terminal temperatures suggest. The standard allowance for that is the correction factor $F \le 1$ applied to the counter-current LMTD, and this unit is `CounterCurrentHX` with that factor in place. Use it when the geometry is a real shell-and-tube bundle rather than an idealised two-stream exchanger; use `CounterCurrentHX` when $F$ would be 1 by construction (a true counter-current double-pipe or a 1-1 arrangement). #### Parameters ```python @dataclass class ShellAndTubeHXParams: UA: float # Overall HTC x area (W/K) Cp_hot: float = None # Hot side heat capacity (J/mol/K), default 75 Cp_cold: float = None # Cold side heat capacity (J/mol/K), default 75 min_approach: float = 10.0 # Minimum temperature approach (K) n_shell_passes: int = 1 # Shell passes: 1-2, 2-4, ... arrangements ``` #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_inlet` | Stream | - | Hot fluid inlet | | `cold_inlet` | Stream | - | Cold fluid inlet | | `UA` | float | W/K | Optional override of `params.UA` | #### Outputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_outlet` | Stream | - | Hot fluid outlet | | `cold_outlet` | Stream | - | Cold fluid outlet | | `info['Q']` | float | W | Heat duty | | `info['F_correction']` | float | - | LMTD correction factor | | `info['R']`, `info['P_param']` | float | - | The two parameters $F$ depends on | | `info['F_too_low']` | bool | - | `F < 0.75`: add shell passes or re-split the duty | | `info['effectiveness']` | float | - | $\epsilon_{CC} F$ | | `info['LMTD']` | float | K | Counter-current LMTD | #### Governing Equations **Corrected duty**: $$Q = UA \cdot F(P, R) \cdot LMTD_{CC}$$ **Correction-factor arguments** (temperature effectiveness and capacity ratio): $$P = \frac{T_{c,out} - T_{c,in}}{T_{h,in} - T_{c,in}}, \qquad R = \frac{T_{h,in} - T_{h,out}}{T_{c,out} - T_{c,in}}$$ **Correction factor** (1-2N TEMA, Bowman-Mueller-Nagle): $$F = \frac{\sqrt{R^2+1}\,\ln\!\frac{1-P}{1-RP}} {(R-1)\,\ln\!\frac{2-P(R+1-\sqrt{R^2+1})}{2-P(R+1+\sqrt{R^2+1})}}$$ At $R = 1$ that expression is $0/0$; `lmtd_correction_factor` takes the L'Hopital limit there, so $F$ and its derivative are finite for balanced flows instead of producing a `nan` in the middle of the usual operating range. $F$ falls as $P$ rises: a multi-pass exchanger asked for a close approach loses area effectiveness quickly, which is why `F_too_low` is reported rather than silently accepted. The conventional design rule is to keep $F > 0.75$ and add shell passes otherwise (`n_shell_passes`, which re-maps $P$ to the equivalent single-shell value before the formula above). #### Example Usage ```python from difflow import make_stream from difflow.units.heat_exchanger import ShellAndTubeHX, ShellAndTubeHXParams hx = ShellAndTubeHX(ShellAndTubeHXParams(UA=500.0, Cp_hot=75.0, Cp_cold=75.0)) hot = make_stream({'water': 2.0}, T=400.0, P=101325.0) cold = make_stream({'water': 1.0}, T=300.0, P=101325.0) hot_out, cold_out, info = hx(hot, cold) print(f"Q = {float(info['Q'])/1000:.2f} kW") print(f"F = {float(info['F_correction']):.3f}") ``` --- ### EnthalpyCounterCurrentHX **Location**: `difflow/units/heat_exchanger.py` **Class**: `EnthalpyCounterCurrentHX` **Description**: Counter-current exchanger closed on real, flash-based stream enthalpies, so it stays correct through a phase change. #### Process Role `CounterCurrentHX` assumes one constant $C_p$ per side, which makes its effectiveness-NTU solution closed-form and its answer wrong wherever the heat capacity is not constant --- above all where a stream boils or condenses, since latent heat is an infinite apparent $C_p$ that a constant-$C_p$ model has no way to represent. This unit closes an enthalpy balance per side instead, using a thermo object's two-phase stream enthalpy, and therefore matches an equation-oriented exchanger (IDAES's, say) through the phase change. Use it for condensers, reboiler-side service, cryogenic/NGL duty, and anywhere a vapour fraction changes across the exchanger. It costs a coupled solve rather than a formula, so prefer `CounterCurrentHX` for single-phase service where the two agree. #### Parameters ```python @dataclass class EnthalpyHXParams: UA: float = None # Overall HTC x area (W/K) max_iter: int = 80 # Iterations of the outer fixed point on Q damping: float = 0.5 # Damping of the Q update (0 < d <= 1) ``` The constructor takes the thermo alongside the params --- `EnthalpyCounterCurrentHX(params, thermo)` --- and the thermo must provide `stream_enthalpy_flash(flows, T, P)`. A `CubicThermo` built from an `IdealThermo` and a `PengRobinson`/`SRK` EOS does ([Thermodynamics](thermodynamics.md)). #### Inputs | Parameter | Type | Units | Description | |-----------|------|-------|-------------| | `hot_inlet` | Stream | - | Hot fluid inlet | | `cold_inlet` | Stream | - | Cold fluid inlet | | `UA` | float | W/K | Optional override of `params.UA` | #### Outputs Same keys as `CounterCurrentHX` where they mean the same thing: `info['Q']`, `info['LMTD']`, the four terminal temperatures, `info['approach']` (the smaller of the two terminal differences) and `info['flow_arrangement'] == 'counter_current_enthalpy'`. #### Governing Equations **Enthalpy balance per side** (not a $\dot{m} C_p \Delta T$ balance): $$H_{h,out} = H_{h,in} - Q, \qquad H_{c,out} = H_{c,in} + Q$$ **Heat transfer**: $$Q = UA \cdot LMTD(T_{h,in}, T_{h,out}, T_{c,in}, T_{c,out})$$ The unknowns $(Q, T_{h,out}, T_{c,out})$ are coupled: $Q$ sets the outlet enthalpies, the enthalpies set the outlet temperatures, and the temperatures set the LMTD that sets $Q$. The unit solves it as a damped fixed point on $Q$ with a one-dimensional enthalpy inversion per side (enthalpy is monotone in $T$, so each inversion is a well-posed root find). Every solve is an `optimistix` root find, so gradients come from the implicit function theorem at the converged result rather than from differentiating the iteration. #### Example Usage ```python from difflow import make_stream from difflow.eos import PengRobinson, CriticalProperties from difflow.thermo import IdealThermo, CubicThermo, SpeciesData from difflow.units.heat_exchanger import ( EnthalpyCounterCurrentHX, EnthalpyHXParams, ) species = { "propane": SpeciesData(name="propane", MW=44.10, Cp_coeffs=(73.0, 0.0, 0.0, 0.0), Hvap_coeffs=(18000.0, 0.38, 369.8), antoine_coeffs=(13.72, 1872.5, -25.16)), "butane": SpeciesData(name="butane", MW=58.12, Cp_coeffs=(98.0, 0.0, 0.0, 0.0), Hvap_coeffs=(22000.0, 0.38, 425.1), antoine_coeffs=(13.98, 2292.4, -27.86)), } crit = { "propane": CriticalProperties(name="propane", Tc=369.8, Pc=4.25e6, omega=0.152, MW=44.10), "butane": CriticalProperties(name="butane", Tc=425.1, Pc=3.80e6, omega=0.200, MW=58.12), } thermo = CubicThermo(IdealThermo(species), PengRobinson(crit)) hot = make_stream({"propane": 1.0, "butane": 1.0}, T=400.0, P=3e5) cold = make_stream({"propane": 1.0, "butane": 1.0}, T=300.0, P=3e5) hx = EnthalpyCounterCurrentHX(EnthalpyHXParams(UA=200.0), thermo) hot_out, cold_out, info = hx(hot, cold) print(f"Q = {float(info['Q'])/1000:.2f} kW") ``` #### Design Considerations - **Share the thermo object.** The coupled solve is JIT-compiled and cached on the identity of the thermo, so passing the same `CubicThermo` to every exchanger compiles once instead of once per unit. The first call takes seconds; later ones do not. - **`damping` is a convergence knob, not a model parameter.** Gradients are exact at the fixed point regardless of its value, so lowering it costs iterations and nothing else. --- ### Heat Exchanger Utility Functions ```python from difflow.units.heat_exchanger import ( log_mean_temperature_difference, effectiveness_counter_current, effectiveness_co_current, effectiveness_crossflow_both_unmixed, effectiveness_crossflow_cmax_mixed, effectiveness_crossflow_cmin_mixed, effectiveness_crossflow_both_mixed, heat_capacity_rate, design_heat_exchanger, size_heat_exchanger ) # Calculate LMTD with numerical stability lmtd = log_mean_temperature_difference(dT1=50.0, dT2=30.0) # Calculate effectiveness for different flow configurations eps_counter = effectiveness_counter_current(NTU=2.0, Cr=0.5) eps_co = effectiveness_co_current(NTU=2.0, Cr=0.5) eps_cross = effectiveness_crossflow_both_unmixed(NTU=2.0, Cr=0.5) # Design for specified duty design = design_heat_exchanger(Q=100000.0, T_hot_in=400.0, T_hot_out=350.0, T_cold_in=300.0, T_cold_out=340.0, U=500.0) design['UA'], design['A'], design['LMTD'] # W/K, m², K ``` --- ## Liquid-Liquid Extraction ### LLEEquilibrium **Location**: `difflow/units/lle.py` **Class**: `LLEEquilibrium` **Description**: The two-phase equilibrium model every extraction unit takes as a parameter: which species transfer, which carry each phase, and how the distribution coefficients are computed. #### Process Role `LLEEquilibrium` is a model object, not a unit operation: it has no inlets and no outlets, and it never appears in a flowsheet on its own. It is the `equilibrium` field of `CascadeParams` and `DifferentialContactorParams`, and what it decides is the part of an extraction calculation that is thermodynamics rather than cascade bookkeeping: - **which species are solutes** (they partition between the phases) and which are the two **carriers** (they do not, so their flows set the phase ratio at every stage); - **where the distribution coefficients come from** --- either tabulated $K$ values, optionally with a van't Hoff temperature dependence, or an activity-coefficient model (NRTL or UNIQUAC) from which $K_i = \gamma_i^{aq} / \gamma_i^{org}$ follows. Because it is a model object, the editor cannot construct one from a form alone: `solutes`, `aqueous_carrier` and `organic_carrier` have no defaults, and the palette reports them as unmet requirements until the code context binds them. #### Parameters ```python @dataclass class LLEEquilibrium: solutes: list[str] # Species that transfer between phases aqueous_carrier: str # Species that stays aqueous organic_carrier: str # Species that stays organic K_coeffs: DistributionCoeffs = None # For activity_model='K' nrtl_params: NRTLParams = None # For activity_model='NRTL' uniquac_params: UNIQUACParams = None# For activity_model='UNIQUAC' activity_model: str = 'K' # 'K' | 'NRTL' | 'UNIQUAC' mutual_solubility: dict = None # Optional carrier cross-solubility ``` `DistributionCoeffs` holds `species`, `K0` at `Tref`, and optionally the heats of extraction `dH` that make $K$ temperature-dependent: $$K_i(T) = K_{0,i} \exp\left[-\frac{\Delta H_i}{R} \left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right]$$ #### Governing Equations **Distribution coefficient** (the sign convention: $K > 1$ favours the organic/extract phase): $$K_i = \frac{y_i^{org}}{x_i^{aq}}$$ **Component balance across a contact**: $$F z_i = E\, y_i^{org} + R\, x_i^{aq}$$ **Isoactivity**, when an activity model is used instead of tabulated $K$: $$\gamma_i^{aq} x_i^{aq} = \gamma_i^{org} y_i^{org} \qquad \Rightarrow \qquad K_i = \frac{\gamma_i^{aq}}{\gamma_i^{org}}$$ #### Example Usage ```python from difflow.units.lle import LLEEquilibrium, DistributionCoeffs equilibrium = LLEEquilibrium( solutes=["acetic_acid"], aqueous_carrier="water", organic_carrier="butanol", K_coeffs=DistributionCoeffs(species=("acetic_acid",), K0=(2.5,)), ) # The one thing it computes, at a temperature: equilibrium.get_distribution_coefficients({}, {}, 298.15) # {'acetic_acid': 2.5} ``` With `activity_model='NRTL'` the two composition arguments are used (the coefficients depend on them); with `'K'` they are ignored, which is why they can be empty above. --- ### MultistageCascade **Location**: `difflow/units/lle.py` **Class**: `MultistageCascade` **Description**: Counter-current multistage extraction cascade for liquid-liquid separation. #### Process Role LLE is used for: - Separation of heat-sensitive compounds - Aromatics extraction (BTX) - Pharmaceutical purification - Metal extraction (hydrometallurgy) #### Parameters ```python @dataclass class CascadeParams: n_stages: int | float # Equilibrium stages (continuous, for optimization) equilibrium: LLEEquilibrium # The equilibrium model (see above) flow_config: str = 'counter_current' # or 'co_current' stage_efficiency: float = 0.8 # Murphree efficiency, co-current only ``` `n_stages` is deliberately continuous: the Kremser solution below is smooth in it, so stage count is an ordinary design variable that `jax.grad` can differentiate rather than an integer to enumerate. #### Governing Equations **Distribution Coefficient**: $$K_i = \frac{C_{i,extract}}{C_{i,raffinate}} = \frac{y_i}{x_i}$$ **Material Balance** (stage j): $$R_{j-1} x_{i,j-1} + E_{j+1} y_{i,j+1} = R_j x_{i,j} + E_j y_{i,j}$$ **Equilibrium**: $$y_{i,j} = K_i(T) x_{i,j}$$ **Kremser Equation** (for dilute systems): $$\frac{x_{in} - x_{out}}{x_{in} - x_{out}^*} = \frac{A^{N+1} - A}{A^{N+1} - 1}$$ Where $A = KE/R$ is the extraction factor. #### Activity Coefficient Models **NRTL**: $$\ln \gamma_i = \frac{\sum_j x_j \tau_{ji} G_{ji}}{\sum_k x_k G_{ki}} + \sum_j \frac{x_j G_{ij}}{\sum_k x_k G_{kj}} \left(\tau_{ij} - \frac{\sum_m x_m \tau_{mj} G_{mj}}{\sum_k x_k G_{kj}}\right)$$ Where: - $G_{ij} = \exp(-\alpha_{ij} \tau_{ij})$ - $\tau_{ij} = (g_{ij} - g_{jj})/RT$ **UNIQUAC**: $$\ln \gamma_i = \ln \gamma_i^C + \ln \gamma_i^R$$ Combinatorial and residual contributions based on molecular size and interaction parameters. #### Example Usage ```python from difflow import make_stream from difflow.streams import get_flows from difflow.units.lle import ( LLEEquilibrium, DistributionCoeffs, MultistageCascade, CascadeParams, ) equilibrium = LLEEquilibrium( solutes=["acetic_acid"], aqueous_carrier="water", organic_carrier="butanol", K_coeffs=DistributionCoeffs(species=("acetic_acid",), K0=(2.5,)), ) cascade = MultistageCascade(CascadeParams(n_stages=5, equilibrium=equilibrium)) feed = make_stream({'water': 100.0, 'acetic_acid': 10.0, 'butanol': 0.0}, T=298.15, P=101325.0) solvent = make_stream({'water': 0.0, 'acetic_acid': 0.0, 'butanol': 50.0}, T=298.15, P=101325.0) raffinate, extract, info = cascade(feed, solvent) recovery = 1.0 - get_flows(raffinate)['acetic_acid'] / 10.0 print(f"Recovery: {float(recovery):.2%}") # Recovery: 91.12% print(info['profiles'].keys()) # stage profiles: x, y, carrier_transfer, ... ``` #### Utility Functions ```python from difflow.units.lle import ( DistributionCoeffs, get_K_values, nrtl_activity_coefficients, uniquac_activity_coefficients, separation_factor, minimum_solvent_ratio, stages_for_recovery ) # Distribution coefficients at a temperature (tabulated K0 at Tref) coeffs = DistributionCoeffs(species=("acetic_acid",), K0=(2.5,)) K = get_K_values(coeffs, 298.15)["acetic_acid"] # Minimum solvent-to-feed ratio for a given recovery S_min = minimum_solvent_ratio(K, recovery=0.95) # Stages needed (Kremser) at an actual solvent-to-feed ratio above the minimum N = stages_for_recovery(K, 1.5 * S_min, recovery=0.95) ``` --- ### DifferentialContactor **Location**: `difflow/units/lle.py` **Class**: `DifferentialContactor` **Description**: Continuous differential contact extraction column (spray, packed, or rotating disc). #### Governing Equations **Height of Transfer Unit (HTU)**: $$HTU = \frac{R}{K_{OC} a A}$$ **Number of Transfer Units (NTU)**: $$NTU = \int_{x_{out}}^{x_{in}} \frac{dx}{x - x^*}$$ **Column Height**: $$H = HTU \times NTU$$ --- ## Pressure-Change & EOS-Consistent Units These units close energy balances on the **cubic-EOS enthalpy and entropy** (ideal-gas Cp + Peng-Robinson/SRK departures) rather than on ideal-K or constant-Cp models, so they are correct for real gases and near-cryogenic / gas-processing service (expander plants, NGL recovery, refrigeration). Each takes a [`CubicThermo`](thermodynamics.md) built from an `IdealThermo` (for the ideal-gas Cp) and a `PengRobinson`/`SRK` EOS. All internal temperature solves use `optimistix` root finds on the two-phase enthalpy/entropy, so every outlet temperature, duty and shaft work is differentiable with respect to feed conditions, discharge pressures and efficiencies. ```python from difflow import ( IdealThermo, CubicThermo, PengRobinson, Turboexpander, TurboexpanderParams, Compressor, CompressorParams, JTValve, JTValveParams, ComponentSeparator, ComponentSeparatorParams, ) from difflow.database import get_critical_props, get_species_data from difflow.streams import make_stream names = ["nitrogen", "methane", "ethane", "propane", "n_butane"] ideal = IdealThermo({c: get_species_data(c) for c in names}) eos = PengRobinson({c: get_critical_props(c) for c in names}) thermo = CubicThermo(ideal, eos) feed = make_stream({"nitrogen": 0.5, "methane": 86.0, "ethane": 7.0, "propane": 3.0, "n_butane": 1.0}, T=305.0, P=60e5) ``` ### Turboexpander Adiabatic expansion to `P_out` with an isentropic efficiency. The reversible outlet is found by matching entropy, then the efficiency is applied to the enthalpy drop: $$S(T_\text{isen}, P_\text{out}) = S(T_\text{in}, P_\text{in}), \qquad H_\text{out} = H_\text{in} + \eta\,(H_\text{isen} - H_\text{in})$$ The extracted shaft work is $W = H_\text{in} - H_\text{out} > 0$. Both enthalpy and entropy are two-phase aware, so an expander whose outlet partly condenses (common in cryogenic service) is handled correctly. ```python exp = Turboexpander(TurboexpanderParams(P_out=20e5, eta_isentropic=0.80), thermo) outlet, info = exp(feed) # info: {"W", "T_isen", "T_out", "H_in", "H_out"} ``` ### EOSCompressor **Exported as** `Compressor`; registered in the catalog --- and so labelled in the editor's palette --- as **`EOSCompressor`**, because `Compressor` is taken by the gas plugin's fixed-ratio compressor station ([Gas networks](unit-operations-gas.md)) and the two are different models of different things. Adiabatic compression to `P_out` with an isentropic efficiency. Same entropy match, but the efficiency **inflates** the enthalpy rise (an inefficient machine needs more work than the reversible one): $$H_\text{out} = H_\text{in} + \frac{H_\text{isen} - H_\text{in}}{\eta}$$ The required shaft work is $W = H_\text{out} - H_\text{in} > 0$. ```python comp = Compressor(CompressorParams(P_out=90e5, eta_isentropic=0.75), thermo) outlet, info = comp(feed) ``` ### JTValve (Joule-Thomson valve) Adiabatic, **isenthalpic** pressure letdown. Holds the two-phase EOS enthalpy constant across the pressure drop and solves for the outlet temperature, $H(T_\text{out}, P_\text{out}) = H(T_\text{in}, P_\text{in})$. On a real gas this produces the Joule-Thomson temperature change that an ideal-gas or ideal-K valve misses. Because no work is extracted, the same pressure drop cools **less** than a turboexpander. ```python valve = JTValve(JTValveParams(P_out=20e5), thermo) outlet, info = valve(feed) # info: {"T_out", "H"} ``` ### ComponentSeparator A black-box separator surrogate: each component is split to the product stream by a fixed recovery, the complement going to the residue. Both products inherit the inlet T and P, and the reported duty `Q` is the enthalpy imbalance needed to hold both at the inlet temperature. Useful as a column stand-in when only the recovery specification is known. ```python rec = {"propane": 0.95, "n_butane": 0.99} # heavies to product sep = ComponentSeparator( ComponentSeparatorParams(recovery_to_product=rec, default_recovery=0.0), thermo, ) residue, product, info = sep(feed) # info: {"Q", "H_in", "H_out"} ``` --- ## Combustion & Gas-Turbine Units These units model a **Brayton cycle** working fluid — air and combustion gas at high temperature and moderate pressure, where the cubic-EOS departure is negligible. They therefore use **ideal-gas** properties with temperature-dependent Cp (`difflow.combustion.IdealGasThermo`), the standard model for gas-turbine cycle analysis, and are named distinctly from the real-gas `Compressor` / `Turboexpander` above (a different thermodynamic model for a different service). All internal temperature and air/fuel-ratio solves are `optimistix` root finds, so shaft work, firing temperature, air/fuel ratio and efficiency are differentiable with respect to feed conditions, pressures, efficiencies and fuel composition. The fuel-hydrocarbon (and N₂/CO₂) ideal-gas Cp come from the database (the same cubics used by the NGL work); the module adds O₂, H₂O-vapor and Ar Cp, air composition, and per-fuel lower heating values and combustion stoichiometry. ```python from difflow import ( Combustor, CombustorParams, GasCompressor, GasCompressorParams, GasTurbine, GasTurbineParams, brayton_cycle, BraytonCycleParams, make_cycle_thermo, ) from difflow.combustion import AIR_COMPOSITION from difflow.streams import make_stream thermo = make_cycle_thermo() # ideal-gas thermo over cycle + fuel species ``` ### GasCompressor Adiabatic ideal-gas compression to `pressure_ratio × P_in` with an isentropic efficiency (the efficiency inflates the enthalpy rise). Work consumed is $W = H_\text{out} - H_\text{in} > 0$. ```python air = make_stream(dict(AIR_COMPOSITION), T=288.15, P=101325.0) comp = GasCompressor(GasCompressorParams(pressure_ratio=18.0, eta_isentropic=0.89), thermo) compressed, info = comp(air) # info: {"W", "T_isen", "T_out", ...} ``` ### Combustor Complete-combustion reactor, $C_xH_y + (x + y/4)\,O_2 \to x\,CO_2 + (y/2)\,H_2O$, with two modes: - **`"adiabatic"`** — both feeds fixed; solves the adiabatic flame temperature. - **`"fixed_T"`** — scales the air stream to the air/fuel ratio that hits a target firing temperature `T_out` (closed-form, affine in the air amount). ```python fuel = make_stream({"methane": 1.0}, T=298.15, P=18 * 101325.0) comb = Combustor(CombustorParams(mode="fixed_T", T_out=1673.15, dp_frac=0.04), thermo) products, info = comb(fuel, compressed) # info: {"T_out", "air_scale", "o2_demand", "Q"} ``` ### GasTurbine Adiabatic ideal-gas expansion to a back-pressure `P_out` with an isentropic efficiency. Work extracted is $W = H_\text{in} - H_\text{out} > 0$. ```python turb = GasTurbine(GasTurbineParams(P_out=101325.0, eta_isentropic=0.90), thermo) exhaust, info = turb(products) ``` ### brayton_cycle Assembles compressor → combustor → turbine into an intensive (per mole of fuel) simple- or combined-cycle solve. Defaults are a modern F-class machine at ISO conditions and reproduce published performance: simple-cycle η ≈ 0.40 (8470 Btu/kWh), combined-cycle η ≈ 0.568 (6009 Btu/kWh). ```python fuel_comp = {"methane": 0.95, "ethane": 0.03, "propane": 0.01, "nitrogen": 0.005, "carbon_dioxide": 0.005} result = brayton_cycle(fuel_comp, BraytonCycleParams(combined_cycle=True)) # result: {"eta_thermal", "eta_gt_only", "work_net", "air_fuel_molar", ...} ``` --- ## Summary Tables ### Reactor Comparison | Reactor | Mixing | Residence Time | Best For | |---------|--------|----------------|----------| | CSTR | Perfect | Distribution | Liquid-phase, uniform T | | PFR | None (axial) | Uniform | Gas-phase, high X | | GasPFR | None | Variable | Gas with ΔP, mole change | | Fed-Batch | Perfect | Variable | Selectivity control | ### Heat Exchanger Comparison | Type | Arrangement | ΔT Driving Force | Max T Approach | Typical Applications | |------|-------------|------------------|----------------|---------------------| | Counter-current | Opposite flow | Maximum | T_c,out → T_h,in | Max efficiency, heat recovery | | Cross-flow (unmixed) | Perpendicular | Good | Intermediate | Car radiators, HVAC, finned-tube HX | | Cross-flow (mixed) | Perpendicular | Better | Intermediate | Compact HX, special geometries | | Co-current | Parallel flow | Moderate | T_c,out ≤ T_h,out | Simple applications, temperature control | **Note**: For the same UA and inlet conditions, effectiveness ranking is: Counter-current > Cross-flow (mixed) > Cross-flow (unmixed) > Co-current ### Separation Method Selection | Method | Basis | Typical Application | |--------|-------|---------------------| | Flash | VLE | Light/heavy split | | Distillation | Boiling point | High purity, sharp split | | LLE | Solubility | Heat-sensitive, azeotropes |