Chemical Unit Operations

Chemical Unit Operations#

This document provides comprehensive documentation for all chemical unit operations available in Difflow.

Reactors#

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#

@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=...)):

Isothermal (mode='isothermal', T_spec given at call time): Outlet temperature is fixed, heat duty calculated
Adiabatic (mode='adiabatic'): Q = 0, outlet temperature calculated
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#

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:

An equation of state – eos=<cubic 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.
A constant – molar_density=<mol/m^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.

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)#

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#

@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#

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)#

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#

@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#

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#

@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#

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:

from difflow.units.fed_batch import SemiBatchReactor, FedBatchParams

reactor = SemiBatchReactor(fb_params, thermo, mode="isothermal")   # fb_params: FedBatchParams

Everything in 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.

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.

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 returns from a Cantera YAML file, so a published mechanism goes straight into a reactor:

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#

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#

@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.

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 (Non-Ideal)#

Uses fugacity coefficients from cubic equations of state (Peng-Robinson or SRK) for non-ideal VLE.

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 (Isenthalpic)#

Performs adiabatic flash at constant pressure and enthalpy. Solves for flash temperature.

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:

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#

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#

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#

@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:

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#

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:

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#

@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:

\[\begin{split}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}\end{split}\]

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:

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, 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

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#

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:

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#

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 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#

@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 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:

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#

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.

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#

@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#

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#

@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#

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:

Counter-current (highest)
Cross-flow with mixed streams
Cross-flow (both unmixed)
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#

@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#

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#

@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).

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#

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#

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#

@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#

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#

@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#

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#

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 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.

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.

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) 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$.

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.

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.

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.

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$.

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).

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$.

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).

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

	`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

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

Chemical Unit Operations

Contents

Chemical Unit Operations#

Reactors#

CSTR (Continuous Stirred-Tank Reactor)#

Process Role#

Parameters#

Inputs#

Outputs#

Operating Modes#

Governing Equations#

Example Usage#

Concentration Basis (Molar Density)#

Design Considerations#

PFR (Plug Flow Reactor)#

Process Role#

Parameters#

Inputs#

Outputs#

Governing Equations#

Example Usage#

GasPFR (Gas-Phase PFR with Pressure Drop)#

Additional Parameters#

Additional Equations#

Outputs#

FedBatchReactor#

Process Role#

Parameters#

Inputs#

Governing Equations#

Utility Functions#

SemiBatchReactor#

Declarative Kinetics#

The rate law#

What it refuses, and why#

Units#

Separators#

Flash Drum#

Process Role#

Parameters#

Inputs#

Outputs#

Governing Equations#

Flash Classes#

Flash (Ideal)#

EOSFlash (Non-Ideal)#

PHFlash (Isenthalpic)#

Bubble and Dew Point Methods#

Mixer#

Governing Equations#

Example Usage#

Splitter#

Governing Equations#

Example Usage#

Distillation#

ShortcutColumn#

Process Role#

Parameters#

Governing Equations#

Outputs#

Example Usage#

Utility Functions#

DistillationColumn (Rigorous)#

Parameters#

Governing Equations (MESH)#

Product temperatures#

Thermodynamics: ideal K-values or a cubic EOS#

Solver Paths#

Outputs#

Example Usage#

Heat Exchangers#

Heater#

Parameters#

Energy models#

Example Usage#

Cooler#

Parameters#

Governing Equations#

CounterCurrentHX#

Process Role#