CSTR and PFR Basics: Conventional vs Difflow

CSTR and PFR Basics: Conventional vs Difflow#

This notebook demonstrates how difflow solves reactor problems and compares the results to conventional analytical/numerical solutions.

Problem: First-order irreversible reaction A → B

Rate: r = k·C_A
k = 0.5 /s
Inlet: F_A0 = 10 mol/s, F_B0 = 0 mol/s
Volume: V = 2 m³
Volumetric flow: Q = 0.2 m³/s (so τ = V/Q = 10 s)

import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

import numpy as np

# Problem parameters
k = 0.5      # rate constant (1/s)
F_A0 = 10.0  # inlet molar flow of A (mol/s)
F_B0 = 0.0   # inlet molar flow of B (mol/s)
V = 2.0      # reactor volume (m³)
Q = 0.2      # volumetric flow rate (m³/s)
tau = V / Q  # residence time (s)

print(f"Residence time τ = V/Q = {tau:.1f} s")
print(f"Damköhler number Da = k·τ = {k * tau:.1f}")

Residence time τ = V/Q = 10.0 s
Damköhler number Da = k·τ = 5.0

1. CSTR: Continuous Stirred Tank Reactor#

Conventional Solution#

For a CSTR with first-order reaction:

Mole balance: \(F_{A0} = F_A + V·r = F_A + V·k·C_A\)

Since \(C_A = F_A/Q\):

\[F_{A0} = F_A + V \cdot k \cdot \frac{F_A}{Q} = F_A(1 + k\tau)\]

\[F_A = \frac{F_{A0}}{1 + k\tau}\]

\[X = 1 - \frac{F_A}{F_{A0}} = \frac{k\tau}{1 + k\tau}\]

# Conventional CSTR solution
F_A_cstr_conventional = F_A0 / (1 + k * tau)
F_B_cstr_conventional = F_A0 - F_A_cstr_conventional  # by stoichiometry
X_cstr_conventional = (k * tau) / (1 + k * tau)

print("=" * 50)
print("CSTR - Conventional Solution")
print("=" * 50)
print(f"Outlet F_A = {F_A_cstr_conventional:.4f} mol/s")
print(f"Outlet F_B = {F_B_cstr_conventional:.4f} mol/s")
print(f"Conversion X = {X_cstr_conventional*100:.2f}%")

==================================================
CSTR - Conventional Solution
==================================================
Outlet F_A = 1.6667 mol/s
Outlet F_B = 8.3333 mol/s
Conversion X = 83.33%

Difflow Solution#

from difflow import CSTR, CSTRParams, make_stream, get_flows, IdealThermo, SpeciesData

# Define species (minimal properties needed)
species_data = {
    "A": SpeciesData(name="A", MW=100.0, Cp_coeffs=(75.0, 0.0, 0.0, 0.0),
                    Hvap_coeffs=(35000.0, 0.38, 500.0), antoine_coeffs=(10.0, 3000.0, -50.0)),
    "B": SpeciesData(name="B", MW=100.0, Cp_coeffs=(75.0, 0.0, 0.0, 0.0),
                    Hvap_coeffs=(30000.0, 0.38, 450.0), antoine_coeffs=(10.0, 2800.0, -40.0)),
}
thermo = IdealThermo(species_data)

# Rate function: r = k * C_A
def rate_fn(C, T, params):
    return jnp.array([params["k"] * C["A"]])

# Stoichiometry: A → B
stoich = jnp.array([[-1.0], [+1.0]])  # A consumed, B produced

# Create CSTR
cstr_params = CSTRParams(
    V=jnp.array(V),
    rate_fn=rate_fn,
    stoich=stoich,
    rate_params={"k": jnp.array(k)},
    species_order=["A", "B"],
)
cstr = CSTR(cstr_params, thermo=thermo, mode="isothermal")

# Create inlet stream and solve
inlet = make_stream({"A": F_A0, "B": F_B0}, T=300.0, P=101325.0)
outlet, info = cstr(inlet, T_spec=300.0, volumetric_flow=Q)

F_A_cstr_difflow = float(get_flows(outlet)["A"])
F_B_cstr_difflow = float(get_flows(outlet)["B"])
X_cstr_difflow = float(info["conversion"]["A"])

print("=" * 50)
print("CSTR - Difflow Solution")
print("=" * 50)
print(f"Outlet F_A = {F_A_cstr_difflow:.4f} mol/s")
print(f"Outlet F_B = {F_B_cstr_difflow:.4f} mol/s")
print(f"Conversion X = {X_cstr_difflow*100:.2f}%")

WARNING:2026-03-01 07:48:26,933:jax._src.xla_bridge:905: Platform 'METAL' is experimental and not all JAX functionality may be correctly supported!

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
W0000 00:00:1772369306.933781 12631199 mps_client.cc:510] WARNING: JAX Apple GPU support is experimental and not all JAX functionality is correctly supported!
I0000 00:00:1772369306.942727 12631199 service.cc:145] XLA service 0xb42526400 initialized for platform METAL (this does not guarantee that XLA will be used). Devices:
I0000 00:00:1772369306.942736 12631199 service.cc:153]   StreamExecutor device (0): Metal, <undefined>
I0000 00:00:1772369306.943738 12631199 mps_client.cc:406] Using Simple allocator.
I0000 00:00:1772369306.943750 12631199 mps_client.cc:384] XLA backend will use up to 55662313472 bytes on device 0 for SimpleAllocator.
WARNING:2026-03-01 07:48:26,943:jax._src.xla_bridge:905: Platform 'mps' is experimental and not all JAX functionality may be correctly supported!

Metal device set to: Apple M4 Pro

==================================================
CSTR - Difflow Solution
==================================================
Outlet F_A = 1.6667 mol/s
Outlet F_B = 8.3333 mol/s
Conversion X = 83.33%

# Compare results
print("\n" + "=" * 50)
print("CSTR Comparison")
print("=" * 50)
print(f"{'':20} {'Conventional':>15} {'Difflow':>15} {'Diff':>10}")
print("-" * 60)
print(f"{'F_A (mol/s)':20} {F_A_cstr_conventional:>15.4f} {F_A_cstr_difflow:>15.4f} {abs(F_A_cstr_conventional - F_A_cstr_difflow):>10.2e}")
print(f"{'F_B (mol/s)':20} {F_B_cstr_conventional:>15.4f} {F_B_cstr_difflow:>15.4f} {abs(F_B_cstr_conventional - F_B_cstr_difflow):>10.2e}")
print(f"{'Conversion (%)':20} {X_cstr_conventional*100:>15.2f} {X_cstr_difflow*100:>15.2f} {abs(X_cstr_conventional - X_cstr_difflow)*100:>10.2e}")
print("\n✓ Results match!")

==================================================
CSTR Comparison
==================================================
                        Conventional         Difflow       Diff
------------------------------------------------------------
F_A (mol/s)                   1.6667          1.6667   0.00e+00
F_B (mol/s)                   8.3333          8.3333   0.00e+00
Conversion (%)                 83.33           83.33   0.00e+00

✓ Results match!

2. PFR: Plug Flow Reactor#

Conventional Solution#

For a PFR with first-order reaction:

Design equation: \(\frac{dF_A}{dV} = r_A = -k C_A = -k \frac{F_A}{Q}\)

This is a separable ODE:

\[\int_{F_{A0}}^{F_A} \frac{dF_A}{F_A} = -\frac{k}{Q} \int_0^V dV\]

\[\ln\frac{F_A}{F_{A0}} = -\frac{kV}{Q} = -k\tau\]

\[F_A = F_{A0} e^{-k\tau}\]

\[X = 1 - e^{-k\tau}\]

# Conventional PFR solution (analytical)
F_A_pfr_conventional = F_A0 * np.exp(-k * tau)
F_B_pfr_conventional = F_A0 - F_A_pfr_conventional
X_pfr_conventional = 1 - np.exp(-k * tau)

print("=" * 50)
print("PFR - Conventional Solution (Analytical)")
print("=" * 50)
print(f"Outlet F_A = {F_A_pfr_conventional:.4f} mol/s")
print(f"Outlet F_B = {F_B_pfr_conventional:.4f} mol/s")
print(f"Conversion X = {X_pfr_conventional*100:.2f}%")

==================================================
PFR - Conventional Solution (Analytical)
==================================================
Outlet F_A = 0.0674 mol/s
Outlet F_B = 9.9326 mol/s
Conversion X = 99.33%

Difflow Solution#

from difflow import PFR, PFRParams

# Create PFR (uses same rate_fn and stoich as CSTR)
pfr_params = PFRParams(
    V=jnp.array(V),
    rate_fn=rate_fn,
    stoich=stoich,
    rate_params={"k": jnp.array(k)},
    species_order=["A", "B"],
    n_save_points=101,  # Points to save in output profile
)
pfr = PFR(pfr_params, thermo=thermo, mode="isothermal")

# Solve PFR
outlet_pfr, info_pfr = pfr(inlet, T_spec=300.0, volumetric_flow=Q)

F_A_pfr_difflow = float(get_flows(outlet_pfr)["A"])
F_B_pfr_difflow = float(get_flows(outlet_pfr)["B"])
X_pfr_difflow = float(info_pfr["conversion"]["A"])

print("=" * 50)
print("PFR - Difflow Solution")
print("=" * 50)
print(f"Outlet F_A = {F_A_pfr_difflow:.4f} mol/s")
print(f"Outlet F_B = {F_B_pfr_difflow:.4f} mol/s")
print(f"Conversion X = {X_pfr_difflow*100:.2f}%")

==================================================
PFR - Difflow Solution
==================================================
Outlet F_A = 0.0674 mol/s
Outlet F_B = 9.9326 mol/s
Conversion X = 99.33%

# Compare results
print("\n" + "=" * 50)
print("PFR Comparison")
print("=" * 50)
print(f"{'':20} {'Conventional':>15} {'Difflow':>15} {'Diff':>10}")
print("-" * 60)
print(f"{'F_A (mol/s)':20} {F_A_pfr_conventional:>15.4f} {F_A_pfr_difflow:>15.4f} {abs(F_A_pfr_conventional - F_A_pfr_difflow):>10.2e}")
print(f"{'F_B (mol/s)':20} {F_B_pfr_conventional:>15.4f} {F_B_pfr_difflow:>15.4f} {abs(F_B_pfr_conventional - F_B_pfr_difflow):>10.2e}")
print(f"{'Conversion (%)':20} {X_pfr_conventional*100:>15.2f} {X_pfr_difflow*100:>15.2f} {abs(X_pfr_conventional - X_pfr_difflow)*100:>10.2e}")
print("\n✓ Results match!")

==================================================
PFR Comparison
==================================================
                        Conventional         Difflow       Diff
------------------------------------------------------------
F_A (mol/s)                   0.0674          0.0674   6.33e-08
F_B (mol/s)                   9.9326          9.9326   6.33e-08
Conversion (%)                 99.33           99.33   6.33e-07

✓ Results match!

3. CSTR vs PFR Comparison#

For positive-order reactions (like first-order), PFR always gives higher conversion than CSTR at the same volume and flow rate.

print("=" * 50)
print("CSTR vs PFR (same V, Q, k)")
print("=" * 50)
print(f"CSTR conversion: {X_cstr_conventional*100:.2f}%")
print(f"PFR conversion:  {X_pfr_conventional*100:.2f}%")
print(f"\nPFR advantage: {(X_pfr_conventional - X_cstr_conventional)*100:.2f} percentage points")
print(f"\nFor first-order kinetics:")
print(f"  CSTR: X = kτ/(1+kτ) = {k*tau}/(1+{k*tau}) = {X_cstr_conventional:.4f}")
print(f"  PFR:  X = 1 - e^(-kτ) = 1 - e^(-{k*tau}) = {X_pfr_conventional:.4f}")

==================================================
CSTR vs PFR (same V, Q, k)
==================================================
CSTR conversion: 83.33%
PFR conversion:  99.33%

PFR advantage: 15.99 percentage points

For first-order kinetics:
  CSTR: X = kτ/(1+kτ) = 5.0/(1+5.0) = 0.8333
  PFR:  X = 1 - e^(-kτ) = 1 - e^(-5.0) = 0.9933

4. Bonus: Automatic Differentiation#

Unlike conventional solutions, difflow provides gradients automatically!

# Gradient of conversion with respect to volume
def cstr_conversion(V_val):
    params = CSTRParams(V=V_val, rate_fn=rate_fn, stoich=stoich,
                        rate_params={"k": jnp.array(k)}, species_order=["A", "B"])
    reactor = CSTR(params, thermo=thermo, mode="isothermal")
    out, inf = reactor(inlet, T_spec=300.0, volumetric_flow=Q)
    return inf["conversion"]["A"]

def pfr_conversion(V_val):
    params = PFRParams(V=V_val, rate_fn=rate_fn, stoich=stoich,
                       rate_params={"k": jnp.array(k)}, species_order=["A", "B"], n_save_points=101)
    reactor = PFR(params, thermo=thermo, mode="isothermal")
    out, inf = reactor(inlet, T_spec=300.0, volumetric_flow=Q)
    return inf["conversion"]["A"]

# Compute gradients
dX_dV_cstr = jax.grad(cstr_conversion)(jnp.array(V))
dX_dV_pfr = jax.grad(pfr_conversion)(jnp.array(V))

print("=" * 50)
print("Sensitivity Analysis (Automatic Differentiation)")
print("=" * 50)
print(f"dX/dV for CSTR: {float(dX_dV_cstr):.6f} per m³")
print(f"dX/dV for PFR:  {float(dX_dV_pfr):.6f} per m³")
print(f"\nInterpretation: Increasing volume by 1 m³ increases conversion by:")
print(f"  CSTR: {float(dX_dV_cstr)*100:.2f} percentage points")
print(f"  PFR:  {float(dX_dV_pfr)*100:.2f} percentage points")

==================================================
Sensitivity Analysis (Automatic Differentiation)
==================================================
dX/dV for CSTR: 0.069444 per m³
dX/dV for PFR:  0.016845 per m³

Interpretation: Increasing volume by 1 m³ increases conversion by:
  CSTR: 6.94 percentage points
  PFR:  1.68 percentage points

# Verify gradients against analytical derivatives
# CSTR: X = kτ/(1+kτ) = kV/Q / (1 + kV/Q)
# dX/dV = k/Q / (1 + kV/Q)² = k/Q / (1 + kτ)²
dX_dV_cstr_analytical = (k/Q) / (1 + k*tau)**2

# PFR: X = 1 - exp(-kτ) = 1 - exp(-kV/Q)
# dX/dV = (k/Q) * exp(-kV/Q) = (k/Q) * exp(-kτ)
dX_dV_pfr_analytical = (k/Q) * np.exp(-k*tau)

print("\nVerification against analytical gradients:")
print(f"  CSTR: difflow={float(dX_dV_cstr):.6f}, analytical={dX_dV_cstr_analytical:.6f}")
print(f"  PFR:  difflow={float(dX_dV_pfr):.6f}, analytical={dX_dV_pfr_analytical:.6f}")
print("\n✓ Gradients match analytical derivatives!")

Verification against analytical gradients:
  CSTR: difflow=0.069444, analytical=0.069444
  PFR:  difflow=0.016845, analytical=0.016845

✓ Gradients match analytical derivatives!

Summary#

Reactor	Conventional Formula	Difflow Result	Match?
CSTR	X = kτ/(1+kτ) = 83.33%	83.33%	✓
PFR	X = 1 - e^(-kτ) = 99.33%	99.33%	✓

Key advantages of difflow:

Same interface for CSTR and PFR (just change the class)
Automatic gradients via JAX (no manual differentiation)
Works with complex kinetics (just change rate_fn)
Composable into flowsheets with recycles

CSTR and PFR Basics: Conventional vs Difflow

Contents

CSTR and PFR Basics: Conventional vs Difflow#

1. CSTR: Continuous Stirred Tank Reactor#

Conventional Solution#

Difflow Solution#

2. PFR: Plug Flow Reactor#

Conventional Solution#

Difflow Solution#

3. CSTR vs PFR Comparison#

4. Bonus: Automatic Differentiation#

Summary#