Chapters 6 and 7 handled equilibrium: what the state of a supercritical fluid is, and how to compute it. This chapter handles the rates: how fast momentum, mass and heat actually move. Viscosity, diffusivity and thermal conductivity are what set pressure drop, extraction time and heat-exchanger area, and near the critical point they behave in ways that equilibrium intuition gets wrong — thermal conductivity diverges while thermal diffusivity collapses, and viscosity sails straight through the critical point without noticing it. We finish by turning the transport picture into a solvent-selection procedure across the full supercritical solvent set, not just CO₂ and water. Every printed output below was produced by running the code.
Learning Objectives
After completing this chapter, you will be able to:
- Distinguish transport properties from thermodynamic properties, and identify which design quantity each one controls
- Decompose viscosity into a dilute-gas term and a density-dependent residual term, and implement the Lucas correlation for the former
- Explain why viscosity shows essentially no critical anomaly while thermal conductivity shows a large one
- Estimate solute diffusivity in a supercritical fluid from Stokes-Einstein and Wilke-Chang, and state the uncertainty honestly
- Compute and interpret the Schmidt, Prandtl, Reynolds and Sherwood numbers for a supercritical process
- Quantify the critical enhancement of thermal conductivity and the critical slowing down of thermal diffusivity
- Size the pressure drop and the mass-transfer coefficient of a packed supercritical extraction bed, and determine which resistance controls
- Compare CO₂, water, ethanol, propane, nitrogen, xenon and fluorinated fluids at matched reduced states
- Apply a constraint-then-rank procedure to select a supercritical solvent for a stated separation
8.1 Why Transport Properties Need Their Own Chapter
Equilibrium Tells You Whether, Rates Tell You How Long
An equation of state answers the question "at this temperature and pressure, what is the density, and how much solute will dissolve?" It says nothing at all about how long the extraction takes, how much pumping power the bed consumes, or how much heat-exchanger area the cycle needs. Those follow from three coefficients, each of which relates a flux to a gradient:
| Coefficient | Constitutive law | Transports | Sets |
|---|---|---|---|
| Dynamic viscosity $\eta$ (Pa·s) | $\tau = -\eta \, \partial u/\partial y$ | Momentum | Pressure drop, pumping power, flow regime |
| Binary diffusivity $D_{12}$ (m²/s) | $J = -D_{12} \, \partial c/\partial y$ | Mass | Extraction time, chromatographic band width |
| Thermal conductivity $\lambda$ (W/(m·K)) | $q = -\lambda \, \partial T/\partial y$ | Heat | Heat-exchanger area, recuperator effectiveness |
Chapters 1 and 2 gave the qualitative headline: a supercritical fluid has liquid-like density with gas-like viscosity and intermediate diffusivity. That headline is correct as an order of magnitude and it is what makes supercritical processing attractive. It is also, as this chapter will show numerically, considerably less dramatic at the high-density operating conditions industry actually uses than the usual "10-100 times" phrasing suggests.
The Three Regimes of a Transport Coefficient
All three coefficients can be written as a sum of contributions with different physical origins. For viscosity and thermal conductivity the standard decomposition used by every reference correlation is
$$ \eta(T,\rho) = \eta_0(T) + \Delta\eta(T,\rho) + \eta_c(T,\rho) $$ $$ \lambda(T,\rho) = \lambda_0(T) + \Delta\lambda(T,\rho) + \lambda_c(T,\rho) $$- The dilute-gas term $\eta_0$, $\lambda_0$: the zero-density limit, governed by binary collisions and predictable from kinetic theory (Chapman-Enskog) or a corresponding-states correlation. It depends on temperature only.
- The residual (excess) term $\Delta\eta$, $\Delta\lambda$: the effect of finite density. For most fluids this is dominated by density and only weakly dependent on temperature — a fact Example 2 verifies to within a few percent.
- The critical enhancement $\eta_c$, $\lambda_c$: an anomaly driven by long-range density fluctuations whose correlation length $\xi$ diverges at the critical point. This term is the reason supercritical transport properties cannot be extrapolated from either the gas or the liquid side.
The asymmetry that surprises everyone
The critical enhancement is not the same size for the three coefficients. In mode-coupling theory the divergences carry very different exponents: thermal conductivity diverges strongly (roughly as $\xi$), viscosity diverges only logarithmically-weakly (an exponent of order 0.04, which amounts to a few percent within millikelvins of $T_c$), and the mutual diffusion coefficient goes to zero. So at the same state point one coefficient blows up, one is unaffected, and one vanishes. Sections 8.2-8.4 demonstrate all three numerically.
8.2 Viscosity
The Dilute-Gas Term and the Lucas Correlation
At zero density the viscosity of a gas follows from kinetic theory. The practical engineering route is the Lucas corresponding-states correlation, which needs only $T_c$, $P_c$, $M$, $Z_c$ and the dipole moment:
$$ \eta_0 \xi = \left[0.807 T_r^{0.618} - 0.357 e^{-0.449 T_r} + 0.340 e^{-4.058 T_r} + 0.018\right] F_P^{\circ} F_Q^{\circ} $$ $$ \xi = 0.176 \left(\frac{T_c}{M^3 P_c^4}\right)^{1/6} $$with $T_c$ in K, $M$ in g/mol, $P_c$ in bar, $\xi$ in $(\mu\text{P})^{-1}$ and $\eta_0$ in micropoise ($1\,\mu\text{P} = 10^{-7}$ Pa·s). $F_P^{\circ}$ corrects for molecular polarity through the reduced dipole moment $\mu_r = 52.46\, \mu^2 P_c / T_c^2$, and $F_Q^{\circ}$ is a quantum correction needed only for helium, hydrogen and deuterium.
Watch the direction of $\xi$. The group $\xi$ is an inverse viscosity scale, so the correlation is $\eta_0 = [\ldots]/\xi$, not $\eta_0 = [\ldots]\times\xi$. Multiplying instead of dividing gives an answer roughly four orders of magnitude too small — and because the result still looks like a small number with the right kind of exponent, the error survives casual inspection. This is exactly the defect found in the archived source material this chapter was built from.
"""Example 1: dilute-gas viscosity from the Lucas correlation, checked against CoolProp."""
import numpy as np
import CoolProp.CoolProp as CP
R = 8.314462618 # J/(mol.K)
def lucas_dilute_gas_viscosity(T, Tc, Pc_bar, M, Zc, mu_debye=0.0):
"""Lucas correlation for the low-pressure (dilute-gas) viscosity.
Poling, Prausnitz & O'Connell, "The Properties of Gases and Liquids",
5th ed., eqs. 9-4.15 to 9-4.18.
Parameters
----------
T : temperature, K
Tc : critical temperature, K
Pc_bar : critical pressure, bar
M : molar mass, g/mol
Zc : critical compressibility factor, -
mu_debye : dipole moment, debye (0 for non-polar fluids)
Returns
-------
eta : dynamic viscosity, Pa.s
"""
Tr = T / Tc
# Inverse viscosity scaling group, units of 1/micropoise
xi = 0.176 * (Tc / (M ** 3 * Pc_bar ** 4)) ** (1.0 / 6.0)
# Reduced dipole moment and the polarity correction factor F_P
mu_r = 52.46 * mu_debye ** 2 * Pc_bar / Tc ** 2
if mu_r < 0.022:
F_P = 1.0
elif mu_r < 0.075:
F_P = 1.0 + 30.55 * max(0.0, 0.292 - Zc) ** 1.72
else:
F_P = 1.0 + 30.55 * max(0.0, 0.292 - Zc) ** 1.72 * abs(0.96 + 0.1 * (Tr - 0.7))
bracket = (0.807 * Tr ** 0.618 - 0.357 * np.exp(-0.449 * Tr)
+ 0.340 * np.exp(-4.058 * Tr) + 0.018)
eta_micropoise = bracket * F_P / xi
return eta_micropoise * 1e-7 # 1 micropoise = 1e-7 Pa.s
# Dipole moments (debye) from the CRC Handbook
FLUIDS = {
'CO2': 0.0,
'Nitrogen': 0.0,
'Propane': 0.084,
'Ethanol': 1.69,
'Water': 1.85,
}
print("=== Lucas dilute-gas viscosity vs CoolProp (evaluated at Tr = 1.05) ===")
print(f"{'Fluid':10s} {'T (K)':>8s} {'Lucas':>12s} {'CoolProp':>12s} {'error':>8s} {'F_P':>5s}")
print(f"{'':10s} {'':>8s} {'(uPa.s)':>12s} {'(uPa.s)':>12s} {'(%)':>8s}")
print("-" * 62)
for fluid, mu in FLUIDS.items():
Tc = CP.PropsSI('Tcrit', fluid)
Pc = CP.PropsSI('pcrit', fluid)
rho_c = CP.PropsSI('rhocrit', fluid)
M = CP.PropsSI('molar_mass', fluid) * 1000.0 # g/mol
Zc = Pc * (M / 1000.0) / (rho_c * R * Tc)
T = 1.05 * Tc
eta_lucas = lucas_dilute_gas_viscosity(T, Tc, Pc / 1e5, M, Zc, mu)
# Dilute-gas limit of the reference correlation: evaluate at ~zero density
eta_ref = CP.PropsSI('V', 'T', T, 'D', 1e-6, fluid)
err = 100.0 * (eta_lucas - eta_ref) / eta_ref
mu_r = 52.46 * mu ** 2 * (Pc / 1e5) / Tc ** 2
F_P = eta_lucas / lucas_dilute_gas_viscosity(T, Tc, Pc / 1e5, M, Zc, 0.0)
print(f"{fluid:10s} {T:8.2f} {eta_lucas*1e6:12.3f} {eta_ref*1e6:12.3f} "
f"{err:+8.2f} {F_P:5.3f}")
print()
print("Reduced dipole moments (mu_r) and the polarity branch taken:")
for fluid, mu in FLUIDS.items():
Tc = CP.PropsSI('Tcrit', fluid)
Pc_bar = CP.PropsSI('pcrit', fluid) / 1e5
mu_r = 52.46 * mu ** 2 * Pc_bar / Tc ** 2
branch = 'non-polar' if mu_r < 0.022 else ('weakly polar' if mu_r < 0.075 else 'polar')
print(f" {fluid:10s} mu = {mu:4.2f} D -> mu_r = {mu_r:8.5f} ({branch})")Three to four percent for a correlation that needs nothing but the critical constants and a dipole moment is a good result, and it degrades in exactly the expected direction: the two hydrogen-bonding fluids are the worst, and both are handled by the polarity factor rather than by the underlying kinetic theory. Note also that all five dilute-gas viscosities land in the narrow band 9-25 µPa·s. At zero density, all gases are much the same.
Density Dependence: the Residual Term Does the Work
Real supercritical operation is at $\rho/\rho_c \approx 1$-$2$, where the residual term dominates. Subtracting the dilute-gas term isolates it, and the result is close to a function of density alone:
"""Example 2: residual viscosity of CO2 and the (absent) critical spike."""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
FLUID = 'CO2'
Tc = CP.PropsSI('Tcrit', FLUID)
Pc = CP.PropsSI('pcrit', FLUID)
rho_c = CP.PropsSI('rhocrit', FLUID)
def eta0(T):
"""Dilute-gas viscosity: the reference correlation at ~zero density."""
return CP.PropsSI('V', 'T', T, 'D', 1e-6, FLUID)
def eta(T, rho):
return CP.PropsSI('V', 'T', T, 'D', rho, FLUID)
print(f"CO2: Tc = {Tc:.2f} K, Pc = {Pc/1e6:.3f} MPa, rho_c = {rho_c:.1f} kg/m3")
print()
print("=== Residual viscosity d_eta = eta(T, rho) - eta0(T) ===")
print("If d_eta collapses onto a single curve in rho, it is density-controlled,")
print("not temperature-controlled.")
print()
header = f"{'rho/rho_c':>10s}" + ''.join(f"{f'{t-273.15:.0f} C':>12s}" for t in
[310.0, 320.0, 340.0, 380.0, 450.0])
print(header)
print(f"{'':>10s}" + ''.join(f"{'(uPa.s)':>12s}" for _ in range(5)))
print("-" * (10 + 12 * 5))
for rho_r in [0.2, 0.5, 1.0, 1.5, 2.0]:
rho = rho_r * rho_c
row = f"{rho_r:10.1f}"
for T in [310.0, 320.0, 340.0, 380.0, 450.0]:
row += f"{(eta(T, rho) - eta0(T)) * 1e6:12.3f}"
print(row)
print()
print("Spread of d_eta across the 310-450 K range, at fixed density:")
for rho_r in [0.2, 0.5, 1.0, 1.5, 2.0]:
rho = rho_r * rho_c
vals = np.array([eta(T, rho) - eta0(T) for T in [310.0, 320.0, 340.0, 380.0, 450.0]])
print(f" rho/rho_c = {rho_r:3.1f}: mean = {vals.mean()*1e6:7.3f} uPa.s, "
f"spread = {(vals.max()-vals.min())/vals.mean()*100:5.1f} % of mean")
print()
print("=== Is there a viscosity spike at the critical point? ===")
print("Walking an isotherm at Tr = 1.001 straight through rho_c:")
T_near = 1.001 * Tc
print(f"{'rho/rho_c':>10s} {'eta (uPa.s)':>14s} {'d(eta)/d(rho) (uPa.s per kg/m3)':>34s}")
print("-" * 60)
rhos = rho_c * np.array([0.90, 0.95, 0.99, 1.00, 1.01, 1.05, 1.10])
etas = np.array([eta(T_near, r) for r in rhos])
grad = np.gradient(etas, rhos)
for r, e, g in zip(rhos, etas, grad):
print(f"{r/rho_c:10.2f} {e*1e6:14.4f} {g*1e6:34.5f}")
print()
print("Compare with the isobaric heat capacity on the same isotherm, which does diverge:")
for rho_r in [0.90, 1.00, 1.10]:
rho = rho_r * rho_c
cp = CP.PropsSI('C', 'T', T_near, 'D', rho, FLUID)
print(f" rho/rho_c = {rho_r:4.2f}: cp = {cp:10.1f} J/(kg.K), "
f"eta = {eta(T_near, rho)*1e6:7.3f} uPa.s")
# Figure: residual viscosity collapse
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
rho_grid = np.linspace(1.0, 2.2 * rho_c, 300)
for T in [310.0, 320.0, 340.0, 380.0, 450.0]:
e0 = eta0(T)
ax1.plot(rho_grid / rho_c, [(eta(T, r)) * 1e6 for r in rho_grid],
label=f'{T - 273.15:.0f} °C')
ax2.plot(rho_grid / rho_c, [(eta(T, r) - e0) * 1e6 for r in rho_grid],
label=f'{T - 273.15:.0f} °C')
ax1.set_xlabel(r'$\rho/\rho_c$'); ax1.set_ylabel(r'$\eta$ (µPa·s)')
ax1.set_title('Total viscosity'); ax1.legend(); ax1.grid(alpha=0.3)
ax2.set_xlabel(r'$\rho/\rho_c$'); ax2.set_ylabel(r'$\eta-\eta_0(T)$ (µPa·s)')
ax2.set_title('Residual viscosity collapses onto a density curve')
ax2.legend(); ax2.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('co2_residual_viscosity.png', dpi=150)
print("\nFigure saved to 'co2_residual_viscosity.png'")What Example 2 establishes
- Residual viscosity is a density function. Over a 140 K span, $\Delta\eta$ at fixed density varies by 1-3% once $\rho \geq \rho_c$. Density, not temperature, is the control variable for viscosity — which is why the whole chapter is organised around density rather than pressure.
- The dilute-gas term is a small correction, not the answer. At $\rho = 2\rho_c$ the residual term is 85 µPa·s against a dilute-gas term near 17 µPa·s. Estimating supercritical viscosity from kinetic theory alone underpredicts by a factor of six.
- Viscosity does not care about the critical point. On an isotherm 0.3 K above $T_c$, $\eta$ passes through $\rho_c$ smoothly, with a derivative that changes by only 20% across the whole $0.9$-$1.1\,\rho_c$ window. On the same isotherm $c_p$ reaches $5.8\times10^{5}$ J/(kg·K) — a divergence of the kind derived in Chapter 6. Whatever critical enhancement viscosity has is buried inside the correlation's own uncertainty.
The practical consequence is convenient: you can build a supercritical hydraulic model with a viscosity correlation that ignores criticality entirely, and be wrong by less than the uncertainty in your bed voidage. The same shortcut applied to thermal conductivity would be a serious error, as Section 8.4 shows.
8.3 Diffusivity
Three Different Diffusion Coefficients
"Diffusivity" is ambiguous, and in the near-critical region the ambiguity matters:
| Coefficient | Physical meaning | Behaviour near $T_c$ |
|---|---|---|
| Self-diffusion $D_s$ | Random walk of a molecule among identical molecules | Smooth; no anomaly |
| Tracer diffusion $D_{12}^{\infty}$ | Dilute solute in the solvent — the extraction-relevant one | Weak anomaly at infinite dilution |
| Mutual diffusion $D_{12}$ | Relaxation of a composition gradient at finite concentration | Vanishes at the mixture critical locus |
For supercritical extraction at the dilute concentrations typical of a solubility-limited process, the tracer coefficient is the right one, and it is well described hydrodynamically.
Stokes-Einstein: Why Low Viscosity Buys Fast Diffusion
The hydrodynamic result for a sphere of radius $r$ in a continuum of viscosity $\eta$ is
$$ D_{12} = \frac{k_B T}{n \pi \eta r}, \qquad n = 6 \ \text{(no slip)}, \quad n = 4 \ \text{(perfect slip)} $$which contains the entire supercritical selling point in one line: $D \propto T/\eta$, and $\eta$ in a supercritical fluid is well below that of any liquid. It also makes the temperature dependence non-obvious, because raising $T$ raises the numerator and lowers $\eta$ through the density drop, so $D$ rises faster than linearly in $T$.
The classical engineering alternative is the Wilke-Chang correlation quoted in Chapter 2,
$$ D_{AB} = \frac{7.4\times10^{-8}\,(\phi M_B)^{0.5}\,T}{\eta\,V_A^{0.6}} \quad [\text{cm}^2/\text{s},\ \eta \text{ in cP}] $$which was fitted to liquid solvents and is being extrapolated when applied to a supercritical fluid. Comparing the two is a cheap and honest uncertainty estimate.
"""Example 3: solute diffusivity in scCO2 - Stokes-Einstein, Wilke-Chang, Schmidt number."""
import numpy as np
import CoolProp.CoolProp as CP
k_B = 1.380649e-23 # J/K
N_A = 6.02214076e23 # 1/mol
# Naphthalene: the canonical solute for scCO2 diffusion measurements.
V_A = 147.6 # Le Bas molar volume, cm3/mol
M_B = 44.01 # CO2 molar mass, g/mol
def hydrodynamic_radius(V_le_bas_cm3):
"""Equivalent hard-sphere radius from a Le Bas molar volume."""
V_m3 = V_le_bas_cm3 * 1e-6 # m3/mol
return (3.0 * V_m3 / (4.0 * np.pi * N_A)) ** (1.0 / 3.0)
def stokes_einstein(T, eta, radius, n=6.0):
"""D = k_B T / (n pi eta r). n = 6 -> no slip, n = 4 -> perfect slip."""
return k_B * T / (n * np.pi * eta * radius)
def wilke_chang(T, eta_Pa_s, phi=1.0):
"""Wilke-Chang correlation, returned in m2/s.
D_AB = 7.4e-8 (phi M_B)^0.5 T / (eta V_A^0.6) [cm2/s, eta in cP]
"""
eta_cP = eta_Pa_s * 1e3
D_cm2_s = 7.4e-8 * np.sqrt(phi * M_B) * T / (eta_cP * V_A ** 0.6)
return D_cm2_s * 1e-4
r_A = hydrodynamic_radius(V_A)
print(f"Naphthalene: Le Bas molar volume {V_A:.1f} cm3/mol "
f"-> equivalent radius {r_A * 1e10:.2f} A")
print()
print("=== Diffusivity of naphthalene in scCO2 at 40 C ===")
print(f"{'P':>6s} {'rho':>9s} {'eta':>10s} {'D (S-E, 6pi)':>14s} "
f"{'D (Wilke-Chang)':>17s} {'Sc':>7s}")
print(f"{'(MPa)':>6s} {'(kg/m3)':>9s} {'(uPa.s)':>10s} {'(1e-8 m2/s)':>14s} "
f"{'(1e-8 m2/s)':>17s} {'(-)':>7s}")
print("-" * 68)
T = 40 + 273.15
for P_MPa in [8, 10, 15, 20, 25, 30]:
P = P_MPa * 1e6
rho = CP.PropsSI('D', 'T', T, 'P', P, 'CO2')
eta = CP.PropsSI('V', 'T', T, 'P', P, 'CO2')
D_se = stokes_einstein(T, eta, r_A)
D_wc = wilke_chang(T, eta)
Sc = (eta / rho) / D_se
print(f"{P_MPa:6d} {rho:9.1f} {eta*1e6:10.2f} {D_se*1e8:14.3f} "
f"{D_wc*1e8:17.3f} {Sc:7.2f}")
print()
print("=== How do you actually raise D? Two routes, at 15 MPa ===")
print(f"{'T (C)':>7s} {'rho (kg/m3)':>12s} {'eta (uPa.s)':>12s} "
f"{'D (1e-8 m2/s)':>15s} {'Sc':>7s}")
print("-" * 56)
for T_C in [35, 40, 50, 60, 80, 100]:
T_k = T_C + 273.15
rho = CP.PropsSI('D', 'T', T_k, 'P', 15e6, 'CO2')
eta = CP.PropsSI('V', 'T', T_k, 'P', 15e6, 'CO2')
D_se = stokes_einstein(T_k, eta, r_A)
Sc = (eta / rho) / D_se
print(f"{T_C:7d} {rho:12.1f} {eta*1e6:12.2f} {D_se*1e8:15.3f} {Sc:7.2f}")
print()
print("=== Comparison with liquid solvents, same solute ===")
print(f"{'Solvent':22s} {'eta (uPa.s)':>12s} {'D (1e-8 m2/s)':>15s} {'Sc':>9s}")
print("-" * 61)
cases = [
('scCO2, 40 C, 8 MPa', 'CO2', 313.15, 8.0e6),
('scCO2, 40 C, 20 MPa', 'CO2', 313.15, 20.0e6),
('n-Hexane, 25 C', 'n-Hexane', 298.15, 101325.0),
('Ethanol, 25 C', 'Ethanol', 298.15, 101325.0),
('Water, 25 C', 'Water', 298.15, 101325.0),
]
store = {}
for label, fluid, T_l, P_l in cases:
rho_l = CP.PropsSI('D', 'T', T_l, 'P', P_l, fluid)
eta_l = CP.PropsSI('V', 'T', T_l, 'P', P_l, fluid)
D_l = stokes_einstein(T_l, eta_l, r_A)
Sc_l = (eta_l / rho_l) / D_l
store[label] = (eta_l, D_l)
print(f"{label:22s} {eta_l*1e6:12.1f} {D_l*1e8:15.4f} {Sc_l:9.1f}")
print()
print("The 'diffusivity is 10-100x that of a liquid' claim, checked:")
eta_hex, D_hex = store['n-Hexane, 25 C']
for label in ['scCO2, 40 C, 8 MPa', 'scCO2, 40 C, 20 MPa']:
eta_s, D_s = store[label]
print(f" {label:22s}: D / D_hexane = {D_s/D_hex:5.1f} x, "
f"eta_hexane / eta = {eta_hex/eta_s:5.1f} x")The Schmidt number is the number that matters
$Sc = \nu/D_{12}$ compares momentum transport with mass transport, and it is the group that enters every mass-transfer correlation. In scCO₂ it is 3-17. In n-hexane it is 241; in ethanol, 2650. Two to three orders of magnitude in the group that controls the concentration boundary layer is the real quantitative statement behind "supercritical mass transfer is fast", and it is a much sharper statement than a diffusivity ratio, because $Sc$ already accounts for the fact that the supercritical fluid is also less dense.
The "10-100× a liquid" claim needs a pressure attached to it. The last block of Example 3 checks it against n-hexane. At 8 MPa and 40 °C, scCO₂ gives 14× the diffusivity — comfortably inside the usual claim. At 20 MPa and 40 °C, the figure is 3.9×. Both are the same fluid at the same temperature; the difference is 560 kg/m³ of density. Since real extractions run at 15-30 MPa precisely because they need the density for solubility, the modest end of the range is the one that applies to most industrial practice. The advantage is real, but it is a factor of a few, bought at the cost of a pressure vessel.
Where the Hydrodynamic Picture Fails
Two caveats bound the numbers above:
- Absolute accuracy. Stokes-Einstein with a Le Bas hard-sphere radius and Wilke-Chang differ by about 30% here, and both should be treated as order-of-magnitude estimates unless calibrated against measurements for the specific solute. Correlations fitted specifically to supercritical CO₂ (He-Yu; Catchpole-King; Funazukuri and co-workers) do better and are the right tool for design work; they are also fluid-specific and need their published coefficients, which is why this chapter uses the two general-purpose expressions instead.
- The immediate near-critical region. Mutual diffusion slows critically: as the mixture critical locus is approached, $D_{12}$ falls toward zero, exactly opposite to the rise a hydrodynamic correlation would predict from the falling viscosity. Any correlation of the $T/\eta$ form is therefore qualitatively wrong within a few percent of the critical locus. Design mass transfer at $T_r \gtrsim 1.05$, where the anomaly has decayed, and treat the immediate near-critical region as a place to be passed through rather than operated in.
8.4 Thermal Conductivity and the Critical Enhancement
The Term the Other Two Coefficients Do Not Have
Thermal conductivity is the coefficient for which the critical enhancement is unmistakable. Physically, heat is carried not only by molecular collisions but by the collective relaxation of the long-lived density fluctuations whose correlation length $\xi$ grows without bound at the critical point. In the mode-coupling result the enhancement scales roughly as
$$ \lambda_c \sim \frac{k_B T \, \rho \, c_p}{6\pi \eta \, \xi} \cdot \xi \propto \frac{k_B T \rho c_p}{6 \pi \eta} $$so it inherits the divergence of $c_p$ — which Chapter 6 derived — and is large wherever $c_p$ is large. Reference correlations for CO₂ and water include this term explicitly, so it is present in every CoolProp thermal-conductivity call.
"""Example 4: thermal conductivity of CO2 and the critical enhancement."""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
FLUID = 'CO2'
Tc = CP.PropsSI('Tcrit', FLUID)
Pc = CP.PropsSI('pcrit', FLUID)
rho_c = CP.PropsSI('rhocrit', FLUID)
def lam_TP(T, P):
return CP.PropsSI('L', 'T', T, 'P', P, FLUID)
def lam_Trho(T, rho):
return CP.PropsSI('L', 'T', T, 'D', rho, FLUID)
def lam0(T):
"""Dilute-gas thermal conductivity: the correlation at ~zero density."""
return CP.PropsSI('L', 'T', T, 'D', 1e-6, FLUID)
print(f"CO2: Tc = {Tc:.2f} K, Pc = {Pc / 1e6:.3f} MPa, "
f"rho_c = {rho_c:.1f} kg/m3")
print()
print("=== Thermal conductivity along isotherms, mW/(m.K) ===")
temps = [32, 35, 40, 50, 80, 150]
press = [7.5, 8, 9, 10, 15, 20, 30]
print(f"{'P (MPa)':>8s}" + ''.join(f"{f'{t} C':>10s}" for t in temps))
print("-" * (8 + 10 * len(temps)))
for P_MPa in press:
row = f"{P_MPa:8.1f}"
for T_C in temps:
row += f"{lam_TP(T_C + 273.15, P_MPa * 1e6) * 1e3:10.2f}"
print(row)
print()
print("=== The local peak on each isotherm (7.0-13 MPa, 0.002 MPa grid) ===")
print(f"{'T (C)':>7s} {'Tr':>7s} {'P_peak':>9s} {'rho_peak':>10s} "
f"{'lam_peak':>10s} {'lam_min':>10s} {'peak/min':>10s}")
print(f"{'':>7s} {'':>7s} {'(MPa)':>9s} {'(kg/m3)':>10s} "
f"{'(mW/m/K)':>10s} {'(mW/m/K)':>10s} {'':>10s}")
print("-" * 66)
P_scan = np.arange(7.0, 13.0 + 1e-9, 0.002)
for T_C in [31.5, 32, 33, 35, 40, 50, 80]:
T_k = T_C + 273.15
lams = np.array([lam_TP(T_k, p * 1e6) for p in P_scan])
i_pk = int(np.argmax(lams))
if i_pk == 0 or i_pk == len(P_scan) - 1:
print(f"{T_C:7.1f} {T_k / Tc:7.4f} {'--':>9s} {'--':>10s} "
f"{'--':>10s} {'--':>10s} {'monotonic':>10s}")
continue
i_min = i_pk + int(np.argmin(lams[i_pk:]))
rho_peak = CP.PropsSI('D', 'T', T_k, 'P', P_scan[i_pk] * 1e6, FLUID)
print(f"{T_C:7.1f} {T_k / Tc:7.4f} {P_scan[i_pk]:9.3f} {rho_peak:10.1f} "
f"{lams[i_pk] * 1e3:10.2f} {lams[i_min] * 1e3:10.2f} "
f"{lams[i_pk] / lams[i_min]:10.2f}")
print()
print("=== The critical enhancement, isolated at rho = rho_c ===")
print("Approaching Tc from above at fixed critical density:")
print(f"{'T/Tc - 1':>10s} {'T (K)':>9s} {'P (MPa)':>9s} {'lam':>10s} "
f"{'lam0':>9s} {'lam/lam0':>9s}")
print(f"{'':>10s} {'':>9s} {'':>9s} {'(mW/m/K)':>10s} {'(mW/m/K)':>9s} {'':>9s}")
print("-" * 58)
for eps in [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 0.1, 0.5, 1.0]:
T_k = Tc * (1 + eps)
l = lam_Trho(T_k, rho_c)
l0 = lam0(T_k)
P = CP.PropsSI('P', 'T', T_k, 'D', rho_c, FLUID)
print(f"{eps:10.4f} {T_k:9.3f} {P / 1e6:9.3f} {l * 1e3:10.2f} "
f"{l0 * 1e3:9.2f} {l / l0:9.2f}")
lam_water = CP.PropsSI('L', 'T', 298.15, 'P', 101325.0, 'Water')
lam_near = lam_Trho(Tc * 1.0001, rho_c)
print()
print(f"For scale: liquid water at 25 C conducts at "
f"{lam_water * 1e3:.1f} mW/(m.K).")
print(f"CO2 at Tc + 0.03 K and rho_c reaches {lam_near * 1e3:.1f} mW/(m.K), "
f"i.e. {lam_near / lam_water:.2f}x water,")
print(f"while a gas at the same temperature and zero density manages only "
f"{lam0(Tc * 1.0001) * 1e3:.1f} mW/(m.K).")
# Figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
P_grid = np.linspace(6.5, 25, 500)
for T_C in [32, 35, 40, 50, 80]:
T_k = T_C + 273.15
ax1.plot(P_grid, [lam_TP(T_k, p * 1e6) * 1e3 for p in P_grid],
label=f'{T_C} °C')
ax1.axvline(Pc / 1e6, color='k', ls=':', lw=1, label=r'$P_c$')
ax1.set_xlabel('Pressure (MPa)')
ax1.set_ylabel(r'$\lambda$ (mW/(m·K))')
ax1.set_title('Isotherms: a local spike near $P_c$')
ax1.legend(); ax1.grid(alpha=0.3)
eps_grid = np.logspace(-4, 0, 60)
ax2.loglog(eps_grid, [lam_Trho(Tc * (1 + e), rho_c) * 1e3 for e in eps_grid],
'o-', ms=3)
ax2.set_xlabel(r'$T/T_c - 1$')
ax2.set_ylabel(r'$\lambda$ at $\rho_c$ (mW/(m·K))')
ax2.set_title('The critical enhancement diverges')
ax2.grid(alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('co2_thermal_conductivity.png', dpi=150)
print("\nFigure saved to 'co2_thermal_conductivity.png'")Reading the numbers
- The enhancement is real and enormous. Held at critical density, CO₂ reaches 632 mW/(m·K) at 0.03 K above $T_c$ — 37 times the dilute-gas value, and slightly better than liquid water. A fluid that is 99.9% "gas-like" in viscosity terms conducts heat like a liquid.
- It is also extremely local. The enhancement factor falls from 37 to 5.6 by $T_r = 1.01$ and to 4.0 by $T_r = 1.03$. On an isotherm-pressure sweep, the local spike is visible up to about $T_r = 1.013$ (35 °C) and gone by 40 °C, where the isotherm is monotonic in pressure over the whole 7-13 MPa window. Practical extraction conditions therefore sit outside the enhancement; sCO₂ power cycles, which deliberately run near the pseudo-critical line, sit inside it.
- The peak follows the critical density, not a fixed pressure. The local maxima at 31.5-35 °C occur at 450-457 kg/m³, i.e. essentially at $\rho_c = 467.6$ kg/m³, while the pressure at which they occur moves from 7.46 to 8.05 MPa. The enhancement locus is a density locus.
8.5 Prandtl Number, Thermal Diffusivity and Critical Slowing Down
Conductivity Is Not the Same Thing as Heat Penetration
A large $\lambda$ is not automatically good news, because what determines how fast a temperature front moves is the thermal diffusivity
$$ a = \frac{\lambda}{\rho c_p}, \qquad Pr = \frac{\nu}{a} = \frac{\eta c_p}{\lambda} $$Both $\lambda$ and $c_p$ diverge at the critical point, but $c_p$ diverges faster. Their ratio therefore goes to zero: heat diffusion stops precisely where conduction is best. This is critical slowing down, the same phenomenon that makes near-critical temperature control so difficult, and it is directly computable.
"""Example 5: Prandtl number, thermal diffusivity, and critical slowing down."""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
FLUID = 'CO2'
Tc = CP.PropsSI('Tcrit', FLUID)
Pc = CP.PropsSI('pcrit', FLUID)
rho_c = CP.PropsSI('rhocrit', FLUID)
def transport_group(T, P):
"""Return rho, eta, lam, cp, nu, a, Pr at (T, P)."""
rho = CP.PropsSI('D', 'T', T, 'P', P, FLUID)
eta = CP.PropsSI('V', 'T', T, 'P', P, FLUID)
lam = CP.PropsSI('L', 'T', T, 'P', P, FLUID)
cp = CP.PropsSI('C', 'T', T, 'P', P, FLUID)
nu = eta / rho # kinematic viscosity, m2/s
a = lam / (rho * cp) # thermal diffusivity, m2/s
Pr = nu / a # = eta cp / lam
return rho, eta, lam, cp, nu, a, Pr
print("=== Transport groups for CO2 at 40 C ===")
print(f"{'P':>6s} {'rho':>8s} {'eta':>9s} {'lam':>9s} {'cp':>10s} "
f"{'nu':>11s} {'a':>11s} {'Pr':>7s}")
print(f"{'(MPa)':>6s} {'(kg/m3)':>8s} {'(uPa.s)':>9s} {'(mW/m/K)':>9s} "
f"{'(J/kg/K)':>10s} {'(1e-8 m2/s)':>11s} {'(1e-8 m2/s)':>11s} {'(-)':>7s}")
print("-" * 76)
for P_MPa in [8, 10, 15, 20, 25, 30]:
rho, eta, lam, cp, nu, a, Pr = transport_group(313.15, P_MPa * 1e6)
print(f"{P_MPa:6d} {rho:8.1f} {eta*1e6:9.2f} {lam*1e3:9.2f} {cp:10.1f} "
f"{nu*1e8:11.2f} {a*1e8:11.2f} {Pr:7.3f}")
print()
print("=== Critical slowing down of heat transport ===")
print("Thermal diffusivity a = lam / (rho cp) at rho = rho_c, approaching Tc:")
print(f"{'T/Tc - 1':>10s} {'lam':>10s} {'cp':>12s} {'a':>13s} {'Pr':>10s}")
print(f"{'':>10s} {'(mW/m/K)':>10s} {'(J/kg/K)':>12s} {'(1e-8 m2/s)':>13s} {'(-)':>10s}")
print("-" * 58)
for eps in [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 0.1, 0.5]:
T_k = Tc * (1 + eps)
lam = CP.PropsSI('L', 'T', T_k, 'D', rho_c, FLUID)
cp = CP.PropsSI('C', 'T', T_k, 'D', rho_c, FLUID)
eta = CP.PropsSI('V', 'T', T_k, 'D', rho_c, FLUID)
a = lam / (rho_c * cp)
Pr = eta * cp / lam
print(f"{eps:10.4f} {lam*1e3:10.2f} {cp:12.1f} {a*1e8:13.4f} {Pr:10.2f}")
print()
print("Both lambda and cp diverge, but cp diverges faster, so a -> 0.")
print("Heat stops moving at exactly the point where the fluid conducts best.")
print()
print("=== What that means for a heat exchanger ===")
print("Time for heat to diffuse across a 1 mm channel, t ~ L^2 / a:")
L = 1e-3 # m
for label, T_k, P in [
('CO2, 40 C, 20 MPa', 313.15, 20e6),
('CO2, 35 C, 8 MPa (near-critical)', 308.15, 8e6),
('CO2, 32 C, 7.5 MPa (very near-critical)', 305.15, 7.5e6),
('Liquid water, 25 C, 0.1 MPa', 298.15, 101325.0)]:
fluid = 'Water' if 'water' in label else FLUID
lam = CP.PropsSI('L', 'T', T_k, 'P', P, fluid)
rho = CP.PropsSI('D', 'T', T_k, 'P', P, fluid)
cp = CP.PropsSI('C', 'T', T_k, 'P', P, fluid)
a = lam / (rho * cp)
print(f" {label:42s} a = {a*1e8:8.2f}e-8 m2/s -> t = {L**2/a:7.2f} s")
print()
print("=== The pseudo-critical (Widom) locus: where cp peaks along an isobar ===")
print(f"{'P (MPa)':>9s} {'T_pc (C)':>10s} {'cp_max (J/kg/K)':>17s} "
f"{'Pr at T_pc':>12s} {'a (1e-8 m2/s)':>15s}")
print("-" * 66)
for P_MPa in [7.5, 8, 9, 10, 12, 15]:
P = P_MPa * 1e6
T_scan = np.arange(300.0, 400.0, 0.01)
cps = np.array([CP.PropsSI('C', 'T', t, 'P', P, FLUID) for t in T_scan])
i = int(np.argmax(cps))
T_pc = T_scan[i]
rho, eta, lam, cp, nu, a, Pr = transport_group(T_pc, P)
print(f"{P_MPa:9.1f} {T_pc - 273.15:10.2f} {cps[i]:17.1f} {Pr:12.3f} "
f"{a*1e8:15.3f}")
# Figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
P_grid = np.linspace(7.5, 25, 300)
for T_C in [35, 40, 50, 80]:
T_k = T_C + 273.15
ax1.plot(P_grid, [transport_group(T_k, p * 1e6)[6] for p in P_grid],
label=f'{T_C} °C')
ax1.set_xlabel('Pressure (MPa)'); ax1.set_ylabel('Pr (-)')
ax1.set_title('Prandtl number of CO$_2$'); ax1.legend(); ax1.grid(alpha=0.3)
eps_grid = np.logspace(-4, -0.3, 50)
ax2.loglog(eps_grid,
[CP.PropsSI('L', 'T', Tc*(1+e), 'D', rho_c, FLUID)
/ (rho_c * CP.PropsSI('C', 'T', Tc*(1+e), 'D', rho_c, FLUID)) * 1e8
for e in eps_grid], 'o-', ms=3)
ax2.set_xlabel(r'$T/T_c - 1$')
ax2.set_ylabel(r'$a$ at $\rho_c$ ($10^{-8}$ m$^2$/s)')
ax2.set_title('Critical slowing down of heat diffusion')
ax2.grid(alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('co2_prandtl_thermal_diffusivity.png', dpi=150)
print("\nFigure saved to 'co2_prandtl_thermal_diffusivity.png'")The three design consequences
- Near-critical Prandtl numbers are not order unity. $Pr$ rises from 1.7 at 30 MPa to 44 on the pseudo-critical line at 7.5 MPa and to 322 within 0.03 K of $T_c$. Any heat-transfer correlation of the form $Nu = f(Re, Pr)$ that was validated at $Pr \approx 1$ is outside its range there. This is the technical origin of the "deteriorated heat transfer" regime reported in supercritical CO₂ and supercritical water heat-transfer literature.
- Thermal diffusivity collapses by two orders of magnitude. $a$ falls from $7.5\times10^{-8}$ m²/s at $T_r = 1.5$ to $2.2\times10^{-10}$ m²/s at $T_r = 1.0001$. The time for heat to cross a 1 mm channel goes from 20 s at ordinary extraction conditions to 155 s near the critical point — and liquid water, at 6.9 s, beats both.
- The pseudo-critical (Widom) line is where everything happens. Following the $c_p$ maximum along isobars traces a locus from 31.7 °C at 7.5 MPa to 64.3 °C at 15 MPa. It is the locus of maximum $Pr$, minimum $a$, and maximum tunability. sCO₂ Brayton cycles are designed to compress just above it, which buys near-liquid density at near-gas compressibility — and inherits every transport anomaly in this chapter as the price.
Property-independent control tuning will not work near the critical point. Between 8 and 15 MPa on the pseudo-critical line, $a$ changes by a factor of six and $Pr$ by a factor of five. A thermal loop tuned at one operating point will be badly detuned two megapascals away. This is the same conclusion Chapter 6 reached from the $c_p$ divergence, arrived at from the transport side.
8.6 Transport Properties in a Working Extraction Bed
Momentum: the Ergun Equation
A supercritical extraction vessel is a packed bed, and its pressure drop follows the Ergun equation — a viscous term (Kozeny-Carman) plus an inertial term (Burke-Plummer):
$$ \frac{\Delta P}{L} = \underbrace{\frac{150\,\eta\,(1-\varepsilon)^2 u}{\varepsilon^3 d_p^2}}_{\text{viscous}} + \underbrace{\frac{1.75\,\rho\,(1-\varepsilon)u^2}{\varepsilon^3 d_p}}_{\text{inertial}} $$with $u$ the superficial velocity and $\varepsilon$ the bed voidage. Note the two competing density effects: at fixed mass flow, a denser fluid moves more slowly, which cuts both terms, but a denser supercritical fluid is also more viscous.
"""Example 6: pressure drop through a packed extraction bed (Ergun equation)."""
import numpy as np
import CoolProp.CoolProp as CP
# Bed geometry: a laboratory-to-pilot scale extraction vessel
d_p = 1.0e-3 # particle diameter, m (ground coffee / milled plant material)
eps = 0.40 # bed voidage, -
L_bed = 1.0 # bed height, m
D_bed = 0.10 # bed diameter, m
A_bed = np.pi * D_bed ** 2 / 4.0
m_dot = 20.0 / 3600.0 # solvent mass flow, kg/s (20 kg/h)
def ergun(rho, eta, u):
"""Ergun pressure gradient, Pa/m. u = superficial velocity, m/s."""
viscous = 150.0 * eta * (1 - eps) ** 2 * u / (eps ** 3 * d_p ** 2)
inertial = 1.75 * rho * (1 - eps) * u ** 2 / (eps ** 3 * d_p)
return viscous, inertial
print(f"Bed: D = {D_bed*1e3:.0f} mm, L = {L_bed:.1f} m, d_p = {d_p*1e3:.1f} mm, "
f"voidage = {eps:.2f}")
print(f"Solvent mass flow: {m_dot*3600:.0f} kg/h "
f"(cross-section {A_bed*1e4:.1f} cm2)")
print()
print("=== scCO2 at 40 C ===")
print(f"{'P':>6s} {'rho':>8s} {'eta':>9s} {'u':>10s} {'Re_p':>8s} "
f"{'dP visc':>10s} {'dP inert':>10s} {'dP total':>10s}")
print(f"{'(MPa)':>6s} {'(kg/m3)':>8s} {'(uPa.s)':>9s} {'(mm/s)':>10s} {'(-)':>8s} "
f"{'(kPa)':>10s} {'(kPa)':>10s} {'(kPa)':>10s}")
print("-" * 76)
for P_MPa in [8, 10, 15, 20, 30]:
rho = CP.PropsSI('D', 'T', 313.15, 'P', P_MPa * 1e6, 'CO2')
eta = CP.PropsSI('V', 'T', 313.15, 'P', P_MPa * 1e6, 'CO2')
u = m_dot / (rho * A_bed)
Re_p = rho * u * d_p / eta
dv, di = ergun(rho, eta, u)
print(f"{P_MPa:6d} {rho:8.1f} {eta*1e6:9.2f} {u*1e3:10.3f} {Re_p:8.2f} "
f"{dv*L_bed/1e3:10.4f} {di*L_bed/1e3:10.4f} "
f"{(dv+di)*L_bed/1e3:10.4f}")
print()
print("=== Same bed, same mass flow, conventional liquid solvents at 25 C ===")
print(f"{'Solvent':12s} {'rho':>8s} {'eta':>9s} {'u':>10s} {'Re_p':>8s} "
f"{'dP total':>10s} {'ratio to':>10s}")
print(f"{'':12s} {'(kg/m3)':>8s} {'(uPa.s)':>9s} {'(mm/s)':>10s} {'(-)':>8s} "
f"{'(kPa)':>10s} {'scCO2':>10s}")
print("-" * 71)
rho_ref = CP.PropsSI('D', 'T', 313.15, 'P', 20e6, 'CO2')
eta_ref = CP.PropsSI('V', 'T', 313.15, 'P', 20e6, 'CO2')
u_ref = m_dot / (rho_ref * A_bed)
dP_ref = sum(ergun(rho_ref, eta_ref, u_ref)) * L_bed
for label, fluid in [('n-Hexane', 'n-Hexane'), ('Ethanol', 'Ethanol'),
('Water', 'Water')]:
rho = CP.PropsSI('D', 'T', 298.15, 'P', 101325.0, fluid)
eta = CP.PropsSI('V', 'T', 298.15, 'P', 101325.0, fluid)
u = m_dot / (rho * A_bed)
Re_p = rho * u * d_p / eta
dP = sum(ergun(rho, eta, u)) * L_bed
print(f"{label:12s} {rho:8.1f} {eta*1e6:9.1f} {u*1e3:10.3f} {Re_p:8.2f} "
f"{dP/1e3:10.4f} {dP/dP_ref:10.2f}")
print()
print(f"Reference: scCO2 at 40 C / 20 MPa gives dP = {dP_ref/1e3:.4f} kPa over "
f"{L_bed:.1f} m of bed.")
print()
print("=== The pumping-power view ===")
print("Hydraulic power = volumetric flow x pressure drop.")
print(f"{'Case':26s} {'Q (L/h)':>10s} {'dP (kPa)':>10s} {'P_hyd (mW)':>12s}")
print("-" * 60)
cases = [('scCO2, 40 C, 10 MPa', 'CO2', 313.15, 10e6),
('scCO2, 40 C, 20 MPa', 'CO2', 313.15, 20e6),
('n-Hexane, 25 C', 'n-Hexane', 298.15, 101325.0),
('Water, 25 C', 'Water', 298.15, 101325.0)]
for label, fluid, T, P in cases:
rho = CP.PropsSI('D', 'T', T, 'P', P, fluid)
eta = CP.PropsSI('V', 'T', T, 'P', P, fluid)
u = m_dot / (rho * A_bed)
Q = m_dot / rho # m3/s
dP = sum(ergun(rho, eta, u)) * L_bed
print(f"{label:26s} {Q*3.6e6:10.2f} {dP/1e3:10.4f} {Q*dP*1e3:12.4f}")
print()
print("=== Where the bed itself stops being the limitation ===")
print("Fine particles improve mass transfer but the viscous term scales as 1/d_p^2:")
print(f"{'d_p (um)':>10s} {'dP total (kPa)':>16s} {'viscous share':>15s}")
print("-" * 43)
for d_p_um in [2000, 1000, 500, 200, 100, 50]:
d_p = d_p_um * 1e-6
rho, eta = rho_ref, eta_ref
u = m_dot / (rho * A_bed)
dv, di = ergun(rho, eta, u)
print(f"{d_p_um:10d} {(dv+di)*L_bed/1e3:16.3f} "
f"{dv/(dv+di)*100:14.1f}%")Two things stand out. First, the pressure drop is flat in pressure: 0.058-0.077 kPa/m across the whole 8-30 MPa range, because rising viscosity and falling velocity nearly cancel. Bed hydraulics is not a reason to pick an operating pressure. Second, the ratio to a liquid solvent at the same mass flow is 4.3× for hexane and 12.6× for ethanol — real, but four times smaller than a naive comparison of viscosities alone would give, again because the liquid also moves more slowly. And the last block shows where the trade-off bites: milling to 50 µm to speed up intraparticle diffusion raises $\Delta P$ by a factor of 344 and pushes the viscous term to 99% of the total.
Mass: Sherwood Numbers and Which Resistance Controls
External (film) mass transfer in a packed bed is correlated through the Sherwood number. The Wakao-Kaguei correlation is the standard choice:
$$ Sh = \frac{k_f d_p}{D_{12}} = 2 + 1.1\, Sc^{1/3} Re^{0.6}, \qquad Re = \frac{\rho u d_p}{\eta} $$The constant 2 is the stagnant-film limit for an isolated sphere. Comparing the resulting film time constant with the intraparticle diffusion time constant then answers the question that decides how to run the plant.
"""Example 7: external mass transfer in a packed scCO2 extraction bed."""
import numpy as np
import CoolProp.CoolProp as CP
k_B = 1.380649e-23
N_A = 6.02214076e23
# Solute: naphthalene (Le Bas volume 147.6 cm3/mol) -> equivalent radius
r_A = (3.0 * 147.6e-6 / (4.0 * np.pi * N_A)) ** (1.0 / 3.0)
# Bed and particle geometry
d_p = 1.0e-3
eps = 0.40
D_bed = 0.10
A_bed = np.pi * D_bed ** 2 / 4.0
m_dot = 20.0 / 3600.0 # kg/s
tortuosity = 3.0 # typical for a milled plant matrix
porosity_p = 0.5 # intraparticle porosity
def sherwood_wakao_kaguei(Re, Sc):
"""Wakao-Kaguei packed-bed correlation: Sh = 2 + 1.1 Sc^(1/3) Re^0.6.
Valid roughly for 3 < Re < 10000; the constant 2 is the stagnant-film limit.
"""
return 2.0 + 1.1 * Sc ** (1.0 / 3.0) * Re ** 0.6
print("Solute: naphthalene, equivalent radius "
f"{r_A*1e10:.2f} A; bed d_p = {d_p*1e3:.1f} mm, voidage {eps:.2f}")
print("Diffusivity from Stokes-Einstein, D = kT / (6 pi eta r)")
print()
print("=== scCO2 at 40 C, 20 kg/h through a 100 mm bed ===")
print(f"{'P':>6s} {'D':>13s} {'Sc':>7s} {'Re_p':>7s} {'Sh':>7s} "
f"{'k_f':>12s} {'k_f a_v':>11s} {'tau_film':>10s}")
print(f"{'(MPa)':>6s} {'(1e-8 m2/s)':>13s} {'(-)':>7s} {'(-)':>7s} {'(-)':>7s} "
f"{'(1e-5 m/s)':>12s} {'(1/s)':>11s} {'(s)':>10s}")
print("-" * 78)
a_v = 6.0 * (1 - eps) / d_p # interfacial area per unit bed volume, 1/m
rows = []
for P_MPa in [8, 10, 15, 20, 30]:
P = P_MPa * 1e6
rho = CP.PropsSI('D', 'T', 313.15, 'P', P, 'CO2')
eta = CP.PropsSI('V', 'T', 313.15, 'P', P, 'CO2')
D12 = k_B * 313.15 / (6 * np.pi * eta * r_A)
nu = eta / rho
Sc = nu / D12
u = m_dot / (rho * A_bed)
Re = rho * u * d_p / eta
Sh = sherwood_wakao_kaguei(Re, Sc)
k_f = Sh * D12 / d_p
tau_film = 1.0 / (k_f * a_v)
rows.append((P_MPa, D12, Sc, Re, Sh, k_f, tau_film))
print(f"{P_MPa:6d} {D12*1e8:13.3f} {Sc:7.2f} {Re:7.2f} {Sh:7.2f} "
f"{k_f*1e5:12.3f} {k_f*a_v:11.4f} {tau_film:10.1f}")
print()
print(f"Interfacial area a_v = 6(1-eps)/d_p = {a_v:.0f} m2 per m3 of bed")
print()
print("=== Which resistance controls: film or intraparticle? ===")
print("Intraparticle effective diffusivity D_eff = eps_p D / tau,")
print("internal time constant tau_int ~ (d_p/2)^2 / (15 D_eff) for a sphere.")
print(f"{'P (MPa)':>8s} {'D_eff':>13s} {'tau_film':>10s} {'tau_int':>10s} "
f"{'Bi_m':>8s} {'controlling':>14s}")
print(f"{'':>8s} {'(1e-9 m2/s)':>13s} {'(s)':>10s} {'(s)':>10s} {'(-)':>8s}")
print("-" * 68)
for P_MPa, D12, Sc, Re, Sh, k_f, tau_film in rows:
D_eff = porosity_p * D12 / tortuosity
tau_int = (d_p / 2) ** 2 / (15 * D_eff)
Bi_m = k_f * (d_p / 2) / D_eff
which = 'intraparticle' if tau_int > tau_film else 'film'
print(f"{P_MPa:8d} {D_eff*1e9:13.3f} {tau_film:10.1f} {tau_int:10.1f} "
f"{Bi_m:8.1f} {which:>14s}")
print()
print("=== Same bed with a liquid solvent, for scale ===")
print(f"{'Solvent':14s} {'D (1e-8 m2/s)':>15s} {'Sc':>8s} {'Re_p':>7s} "
f"{'Sh':>7s} {'k_f (1e-5 m/s)':>16s}")
print("-" * 70)
for label, fluid, T, P in [('scCO2 (20 MPa)', 'CO2', 313.15, 20e6),
('n-Hexane', 'n-Hexane', 298.15, 101325.0),
('Ethanol', 'Ethanol', 298.15, 101325.0)]:
rho = CP.PropsSI('D', 'T', T, 'P', P, fluid)
eta = CP.PropsSI('V', 'T', T, 'P', P, fluid)
D12 = k_B * T / (6 * np.pi * eta * r_A)
Sc = (eta / rho) / D12
u = m_dot / (rho * A_bed)
Re = rho * u * d_p / eta
Sh = sherwood_wakao_kaguei(Re, Sc)
k_f = Sh * D12 / d_p
print(f"{label:14s} {D12*1e8:15.4f} {Sc:8.1f} {Re:7.2f} {Sh:7.2f} "
f"{k_f*1e5:16.3f}")
print()
print("Note: Re_p here is 0.6-2.4 for the liquids and 8.9-32 for scCO2,")
print("so the Wakao-Kaguei correlation is used slightly below its stated")
print("Re > 3 range for the liquid cases. Treat those two rows as indicative.")Intraparticle diffusion controls, and that changes what you optimise
The mass Biot number is 34-44 across the whole operating range, and the intraparticle time constant exceeds the film time constant by a factor of 4-5 everywhere. In this bed the solvent film is not the bottleneck. The consequences are concrete: raising the CO₂ flow rate mostly wastes solvent, while halving the particle size cuts the internal time constant fourfold. That is why industrial supercritical extraction specifies particle size and matrix pretreatment far more carefully than flow rate — and why Example 6's pressure-drop penalty for fine grinding is the constraint that actually sets the particle size.
Where these numbers stop being trustworthy. $D_{\mathrm{eff}} = \varepsilon_p D_{12}/\tau$ with an assumed intraparticle porosity of 0.5 and tortuosity of 3 is a placeholder for a measurement; real plant matrices vary by an order of magnitude and are frequently anisotropic. The Wakao-Kaguei correlation is also used slightly below its stated $Re > 3$ validity for the liquid comparison rows. The conclusion — intraparticle control, by a comfortable margin — survives these uncertainties, since it would take a 5× error in $D_{\mathrm{eff}}$ to reverse it. The individual time constants do not.
8.7 The Supercritical Solvent Set Beyond CO₂ and Water
Chapters 2 and 3 covered scCO₂ and supercritical water in depth. The remaining candidates get named in comparison tables and then dropped, which understates a useful point: at matched reduced conditions their transport properties are remarkably similar, and the reasons they are not used are almost never transport reasons.
Supercritical Ethanol
$T_c = 241.6$ °C, $P_c = 6.27$ MPa, $\rho_c = 273$ kg/m³. Ethanol occupies the polarity gap between CO₂ and water, retains partial hydrogen-bonding capability, and is renewable and food-grade. Its critical pressure is the lowest of the polar options.
- Polar natural products: polyphenols, glycosides and carotenoids that pure CO₂ cannot reach without a modifier
- Biodiesel: supercritical transesterification at 250-350 °C and 10-20 MPa runs without an alkaline catalyst, in a single phase, and tolerates high free-fatty-acid feedstocks — completing in minutes rather than the hour a catalysed liquid-phase reaction needs
- Biomass delignification and pharmaceutical particle formation
The limitation is flammability at 267 °C: the vessel needs inerting, and the safety case is materially more expensive than for CO₂.
Supercritical Propane
$T_c = 96.7$ °C, $P_c = 4.25$ MPa — the lowest critical pressure of any practical candidate. Propane is a distinctly better lipid solvent than CO₂ at comparable conditions, which is why it survives in vegetable-oil extraction and petroleum deasphalting despite being flammable and explosive. Example 8 shows it has the lowest viscosity and the highest solute diffusivity of the practical set at matched reduced conditions.
Supercritical Nitrogen and Xenon
Nitrogen ($T_c = -147$ °C, $P_c = 3.40$ MPa) is chemically ideal — completely inert, cheap, non-toxic — and practically unusable for extraction, because reaching the supercritical state requires cryogenic plant. It appears in inert-atmosphere processing rather than as a solvent.
Xenon ($T_c = 16.6$ °C, $P_c = 5.84$ MPa) has the most attractive critical temperature of anything on this list, is inert and non-toxic, and is used in specialised pharmaceutical and radiopharmaceutical work. Its cost — thousands of dollars per kilogram — confines it to applications where nothing else will do. It also, as Example 8 demonstrates, has no reference transport model in CoolProp at all, which is itself a lesson: the availability of reference-quality property data is a real engineering constraint on solvent choice, not an afterthought.
Fluorinated Fluids
R-134a ($T_c = 101$ °C, $P_c = 4.06$ MPa) and SF₆ ($T_c = 45.6$ °C, $P_c = 3.76$ MPa) have convenient critical points, zero ozone-depletion potential, well-characterised reference equations, and mature handling practice inherited from refrigeration. They also have global warming potentials of roughly 1430 and 24 300 respectively, which under current phase-down regulation rules them out of any new process where an alternative exists. Their density at matched reduced conditions is the highest of the set, so on transport grounds alone they would be attractive — a clean illustration that solvent selection is not a transport-properties problem.
All of Them at the Same Reduced State
The fair comparison is at matched reduced conditions, because that is what the principle of corresponding states says should collapse the differences.
"""Example 8: transport properties across the supercritical solvent set."""
import numpy as np
import CoolProp.CoolProp as CP
k_B = 1.380649e-23
N_A = 6.02214076e23
r_A = (3.0 * 147.6e-6 / (4.0 * np.pi * N_A)) ** (1.0 / 3.0) # naphthalene
FLUIDS = ['CO2', 'Ethanol', 'Water', 'Propane', 'Nitrogen', 'Xenon',
'R134a', 'SF6']
Tr_TARGET = 1.05
Pr_TARGET = 2.00
def matched_state(fluid, Tr=Tr_TARGET, Pr=Pr_TARGET):
Tc = CP.PropsSI('Tcrit', fluid)
Pc = CP.PropsSI('pcrit', fluid)
return Tr * Tc, Pr * Pc, Tc, Pc
print(f"=== All fluids at the same reduced state: Tr = {Tr_TARGET}, "
f"Pr = {Pr_TARGET} ===")
print(f"{'Fluid':10s} {'Tc':>8s} {'Pc':>7s} {'T':>8s} {'P':>7s} {'rho':>8s} "
f"{'eta':>9s} {'lam':>9s} {'nu':>11s}")
print(f"{'':10s} {'(C)':>8s} {'(MPa)':>7s} {'(C)':>8s} {'(MPa)':>7s} "
f"{'(kg/m3)':>8s} {'(uPa.s)':>9s} {'(mW/m/K)':>9s} {'(1e-8 m2/s)':>11s}")
print("-" * 82)
data = {}
for fluid in FLUIDS:
T, P, Tc, Pc = matched_state(fluid)
rho = CP.PropsSI('D', 'T', T, 'P', P, fluid)
try:
eta = CP.PropsSI('V', 'T', T, 'P', P, fluid)
lam = CP.PropsSI('L', 'T', T, 'P', P, fluid)
except ValueError:
print(f"{fluid:10s} {Tc-273.15:8.2f} {Pc/1e6:7.3f} {T-273.15:8.2f} "
f"{P/1e6:7.3f} {rho:8.1f} "
f"{'no model':>9s} {'no model':>9s} {'--':>11s}")
data[fluid] = None
continue
nu = eta / rho
data[fluid] = (T, P, rho, eta, lam, nu)
print(f"{fluid:10s} {Tc-273.15:8.2f} {Pc/1e6:7.3f} {T-273.15:8.2f} "
f"{P/1e6:7.3f} {rho:8.1f} {eta*1e6:9.2f} {lam*1e3:9.2f} {nu*1e8:11.2f}")
print()
print("=== Derived transport groups at the same reduced state ===")
print(f"{'Fluid':10s} {'cp':>10s} {'a':>12s} {'Pr':>8s} {'D (S-E)':>12s} {'Sc':>8s}")
print(f"{'':10s} {'(J/kg/K)':>10s} {'(1e-8 m2/s)':>12s} {'(-)':>8s} "
f"{'(1e-8 m2/s)':>12s} {'(-)':>8s}")
print("-" * 64)
for fluid in FLUIDS:
if data[fluid] is None:
print(f"{fluid:10s} {'--':>10s} {'--':>12s} {'--':>8s} {'--':>12s} {'--':>8s}")
continue
T, P, rho, eta, lam, nu = data[fluid]
cp = CP.PropsSI('C', 'T', T, 'P', P, fluid)
a = lam / (rho * cp)
Pr = nu / a
D12 = k_B * T / (6 * np.pi * eta * r_A)
Sc = nu / D12
print(f"{fluid:10s} {cp:10.1f} {a*1e8:12.3f} {Pr:8.3f} {D12*1e8:12.3f} "
f"{Sc:8.2f}")
print()
print("=== Ranking (best first) at Tr = 1.05, Pr = 2.00 ===")
avail = [f for f in FLUIDS if data[f] is not None]
by_visc = sorted(avail, key=lambda f: data[f][3])
by_lam = sorted(avail, key=lambda f: -data[f][4])
by_rho = sorted(avail, key=lambda f: -data[f][2])
print("Lowest viscosity (cheapest to pump): " + ", ".join(by_visc))
print("Highest thermal conductivity (heat duty): " + ", ".join(by_lam))
print("Highest density (most solvent power): " + ", ".join(by_rho))
print()
print("=== The catch: reduced states are not equal operating states ===")
print(f"{'Fluid':10s} {'T (C)':>8s} {'P (MPa)':>9s} {'comment'}")
print("-" * 60)
comments = {
'CO2': 'ambient-temperature vessel',
'Ethanol': 'flammable, needs inerting',
'Water': 'corrosion and alloy cost dominate',
'Propane': 'flammable, lowest pressure',
'Nitrogen': 'cryogenic plant required',
'Xenon': 'no reference transport model, and the cost',
'R134a': 'high-GWP refrigerant, phase-down',
'SF6': 'GWP 24 300, restricted use',
}
for fluid in FLUIDS:
T, P, Tc, Pc = matched_state(fluid)
print(f"{fluid:10s} {T-273.15:8.1f} {P/1e6:9.3f} {comments[fluid]}")Corresponding states works, and that is the point
- Viscosity spans a factor of 2.3 across seven chemically unrelated fluids (35-82 µPa·s), and kinematic viscosity a factor of 1.8 (7.1-12.8 $\times 10^{-8}$ m²/s). Prandtl number spans 1.2-2.6 and Schmidt number 5.7-10.0. At matched $T_r$ and $P_r$, all supercritical fluids transport momentum, heat and mass in much the same way.
- The exception is thermal conductivity of water, at 413 mW/(m·K) — five times CO₂. Hydrogen bonding survives into the supercritical state well enough to keep water an outlier in heat transport even where corresponding states has flattened everything else.
- Therefore transport properties are not the discriminator. The absolute temperature and pressure needed to reach the same reduced state, and the safety, regulatory and cost consequences of getting there, differ enormously: 46 °C at 14.8 MPa for CO₂ versus 406 °C at 44.1 MPa for water. That is the real selection problem, and it is what Section 8.8 formalises.
8.8 Choosing a Supercritical Solvent
The Decision Tree
Solute polarity first, then thermal stability, then safety. Transport properties enter last, as a tie-breaker.
The Comparison Table
| Fluid | $T_c$ (°C) | $P_c$ (MPa) | Polarity | Usable density (kg/m³) | Safety | Cost | Principal use |
|---|---|---|---|---|---|---|---|
| CO₂ | 31.0 | 7.38 | Non- to weakly polar (tunable with a modifier) | 200-900 | Excellent: non-toxic, non-flammable, GRAS | Low | Extraction, cleaning, particle formation, power cycles |
| Water | 374.0 | 22.06 | Polar → non-polar above $T_c$ | 50-600 | Good, but severe corrosion and burn hazards | Very low | Oxidation (SCWO), hydrothermal synthesis, gasification |
| Ethanol | 241.6 | 6.27 | Intermediate | 100-500 | Good, but flammable at operating temperature | Medium | Biodiesel, polar natural products, delignification |
| Propane | 96.7 | 4.25 | Non-polar | 150-500 | Flammable and explosive | Low | Vegetable oil extraction, deasphalting |
| Nitrogen | -147.0 | 3.40 | Non-polar | 300-800 | Excellent (inert); asphyxiation risk | Low | Inert-atmosphere processing, not extraction |
| Xenon | 16.6 | 5.84 | Non-polar | 1000-2000 | Excellent (inert) | Very high | Specialised pharmaceutical work; no reference transport model |
| R-134a | 101.1 | 4.06 | Weakly polar | 300-900 | Good; non-flammable | High | Specialised cleaning; restricted by GWP 1430 |
Co-Solvents: Buying Polarity
Adding a few mole percent of a polar modifier is the standard way to move CO₂ into a polarity range it cannot reach alone. Chapter 2 covers the solubility side; the transport side matters too, because the modifier raises viscosity and therefore lowers diffusivity.
| Modifier | Typical loading | Effect | Target solutes |
|---|---|---|---|
| Ethanol | 1-10 vol% (5-15 mol%) | Raises polarity and adds hydrogen bonding | Polyphenols, alkaloids |
| Methanol | 1-5 vol% | Stronger polarity increase than ethanol | Sugars, amino acids |
| Water | 1-10 vol% | Large hydrophilicity increase | Peptides, proteins |
| Acetic acid | 0.1-3 vol% | Acidifies the medium | Basic compounds |
The Hildebrand parameter of pure CO₂ at 40 °C and 20 MPa is around 12 MPa1/2, comparable to hexane; with 10 mol% ethanol it rises to roughly 15-16 MPa1/2, comparable to toluene. The costs are that the modifier must be separated from the product, that selectivity usually falls because unwanted compounds are co-extracted, and that reference mixture property data are thin — CoolProp has no binary interaction data for the CO₂ + ethanol pair, as Chapter 7 demonstrates by raising an exception on it.
Safety and Regulatory Constraints
| Fluid | Principal hazard | Required controls | Regulatory position |
|---|---|---|---|
| CO₂ | High pressure; asphyxiation in poorly ventilated spaces | Relief valves, pressure interlocks, ventilation and CO₂ monitoring | GRAS; widely approved for food and pharmaceutical use |
| Water | Very high pressure; thermal burns; corrosion and salt fouling | Nickel-alloy or lined vessels, insulation, corrosion monitoring | Approved; effluent discharge regulated |
| Ethanol | Flammability at operating temperature | Inert gas blanketing, explosion-proof electrical equipment | Food-grade grades available; flammable-liquid rules apply |
| Propane | Flammability and explosion | Hazardous-area classification, gas detection, automatic isolation | Restricted; flammable-gas regulation |
| R-134a / SF₆ | Asphyxiation; environmental release | Leak detection and recovery, no venting | Restricted or being phased down on GWP grounds |
A Screening Procedure You Can Defend
The right structure separates the hard constraints, which are pass/fail and mostly non-computable, from the ranking, which should use computed properties only. Scoring schemes that add up arbitrary points for temperature, polarity, safety and cost produce a single number that looks objective and cannot be audited; a filter followed by a transparent transport ranking can be.
"""Example 9: transport-aware supercritical solvent screening."""
import numpy as np
import CoolProp.CoolProp as CP
k_B = 1.380649e-23
N_A = 6.02214076e23
r_A = (3.0 * 147.6e-6 / (4.0 * np.pi * N_A)) ** (1.0 / 3.0)
# Non-computable attributes have to be tabulated. Keep them separate from
# anything CoolProp can supply, so it stays obvious which numbers are
# judgement calls and which are reference-quality property data.
CANDIDATES = {
'CO2': dict(polarity='non-polar', flammable=False, cost='low',
cosolvent=True,
notes='GRAS, non-flammable, easy depressurisation'),
'Ethanol': dict(polarity='polar', flammable=True, cost='medium',
cosolvent=False,
notes='renewable, food-grade, needs inerting'),
'Water': dict(polarity='polar', flammable=False, cost='very low',
cosolvent=False,
notes='corrosive above Tc, high-alloy vessels'),
'Propane': dict(polarity='non-polar', flammable=True, cost='low',
cosolvent=False,
notes='better lipid solvent than CO2, explosion risk'),
'Nitrogen': dict(polarity='non-polar', flammable=False, cost='low',
cosolvent=False,
notes='inert, but Tc = -147 C'),
'R134a': dict(polarity='weak', flammable=False, cost='high',
cosolvent=False,
notes='GWP 1430, under phase-down'),
}
POLARITY_OK = {
'non-polar': {'non-polar', 'weak'},
'polar': {'polar'},
}
def screen(max_T_C, required_polarity, allow_flammable,
allow_cosolvent=True, Tr=1.05, Pr=2.0):
"""Hard-filter the candidates, then rank the survivors on transport.
Constraints are pass/fail; the ranking uses computed properties only.
"""
passed, rejected = [], []
for fluid, meta in CANDIDATES.items():
Tc = CP.PropsSI('Tcrit', fluid)
Pc = CP.PropsSI('pcrit', fluid)
T_op, P_op = Tr * Tc, Pr * Pc
reasons, flags = [], []
if Tc < 273.15:
reasons.append('cryogenic critical temperature')
if T_op - 273.15 > max_T_C:
reasons.append(f'operating T {T_op-273.15:.0f} C > limit {max_T_C} C')
if meta['flammable'] and not allow_flammable:
reasons.append('flammable')
if meta['polarity'] not in POLARITY_OK[required_polarity]:
if required_polarity == 'polar' and meta['cosolvent'] and allow_cosolvent:
flags.append('needs a polar co-solvent (5-10 mol% ethanol)')
else:
reasons.append(f"{meta['polarity']} solvent for a "
f"{required_polarity} solute")
if reasons:
rejected.append((fluid, reasons))
continue
try:
rho = CP.PropsSI('D', 'T', T_op, 'P', P_op, fluid)
eta = CP.PropsSI('V', 'T', T_op, 'P', P_op, fluid)
except ValueError:
rejected.append((fluid, ['no reference transport model available']))
continue
D12 = k_B * T_op / (6 * np.pi * eta * r_A)
passed.append(dict(fluid=fluid, T=T_op - 273.15, P=P_op / 1e6, rho=rho,
eta=eta, D12=D12, Sc=(eta / rho) / D12,
cost=meta['cost'], flags=flags))
passed.sort(key=lambda r: -r['D12'])
print(f" Constraints: T <= {max_T_C} C, {required_polarity} solute, "
f"flammable {'allowed' if allow_flammable else 'not allowed'}, "
f"co-solvent {'allowed' if allow_cosolvent else 'not allowed'}")
print(f" Screening state: Tr = {Tr}, Pr = {Pr}")
if passed:
print(f" {'rank':>4s} {'fluid':10s} {'T (C)':>7s} {'P (MPa)':>8s} "
f"{'rho':>8s} {'eta':>8s} {'D (1e-8)':>9s} {'Sc':>6s} {'cost':>9s}")
for i, r in enumerate(passed, 1):
print(f" {i:4d} {r['fluid']:10s} {r['T']:7.1f} {r['P']:8.2f} "
f"{r['rho']:8.1f} {r['eta']*1e6:8.2f} {r['D12']*1e8:9.3f} "
f"{r['Sc']:6.2f} {r['cost']:>9s}")
for f in r['flags']:
print(f" -> {f}")
else:
print(" No candidate satisfies the constraints.")
for fluid, reasons in rejected:
print(f" rejected {fluid:10s} -- {'; '.join(reasons)}")
return passed, rejected
print("=== Case 1: heat-sensitive natural product, non-polar, no flammables ===")
screen(max_T_C=80, required_polarity='non-polar', allow_flammable=False)
print()
print("=== Case 2: polyphenol extraction, polar solute, no flammables ===")
screen(max_T_C=120, required_polarity='polar', allow_flammable=False)
print()
print("=== Case 3: lipid extraction, flammables acceptable ===")
screen(max_T_C=150, required_polarity='non-polar', allow_flammable=True)
print()
print("=== Case 4: hydrothermal oxidation, polar, high temperature allowed ===")
screen(max_T_C=500, required_polarity='polar', allow_flammable=False,
allow_cosolvent=False)
print()
print("Case 1 leaves exactly one survivor, which is the honest answer: scCO2")
print("dominates industrial practice because the constraints that matter")
print("commercially eliminate everything else before any transport comparison")
print("is made. Transport properties decide the outcome only once two or more")
print("candidates clear the constraints, as in Case 3.")What the screening procedure actually shows
- Case 1 has exactly one survivor. For a heat-sensitive non-polar target with no flammables allowed, only CO₂ clears the constraints. This is not a defect of the tool; it is the reason scCO₂ holds the overwhelming majority of industrial supercritical applications.
- Case 2 surfaces the co-solvent route. A polar solute under an 120 °C ceiling eliminates ethanol and water on temperature, so the answer is not another fluid but CO₂ with a polar modifier — which is exactly what commercial polyphenol extraction does.
- Case 3 is where transport decides. With flammables permitted, three candidates clear the constraints and the ranking is a genuine transport ranking: propane first on diffusivity, at 8.5 MPa rather than CO₂'s 14.8 MPa.
- Constraints dominate. Across all four cases, no candidate is ever eliminated by transport properties. They are eliminated by temperature, flammability, polarity or missing data. Transport properties tell you how well a feasible process will run; they do not tell you which processes are feasible.
Summary
Key Takeaways
1. Structure of a transport coefficient
- $\eta = \eta_0(T) + \Delta\eta(\rho) + \eta_c$; the same decomposition holds for $\lambda$.
- The Lucas correlation gives $\eta_0$ within 1-4% from critical constants and a dipole moment.
- The residual term is a function of density to within 1-3% over a 140 K span, and dominates at operating conditions.
2. The critical anomalies are not alike
- Viscosity: no usable anomaly. It passes through the critical point smoothly.
- Thermal conductivity: 37× enhancement at $\rho_c$, 0.03 K above $T_c$; down to 4× by $T_r = 1.03$, and the local spike on a pressure sweep has disappeared by 40 °C.
- Thermal diffusivity: collapses to $2\times10^{-10}$ m²/s — critical slowing down. $Pr$ reaches 322.
- Mutual diffusion: vanishes at the mixture critical locus, so $T/\eta$ correlations are qualitatively wrong there.
3. Diffusivity and mass transfer
- $D \propto T/\eta$; raising temperature at fixed pressure raises $D$ through both terms.
- Schmidt number in scCO₂ is 3-17 against 241 for hexane and 2650 for ethanol. That ratio, not the diffusivity ratio, is the quantitative case for supercritical mass transfer.
- The "10-100× a liquid" claim holds at low density (14× at 8 MPa) and shrinks to 3.9× at 20 MPa.
4. In a real bed
- Ergun pressure drop is nearly flat in pressure (0.058-0.077 kPa/m from 8 to 30 MPa) and 4-13× lower than with a liquid solvent.
- Intraparticle diffusion controls ($Bi_m$ = 34-44), so particle size matters far more than flow rate.
- Grinding finer costs pressure drop as $1/d_p^2$; that trade-off, not mass transfer, sets the particle size.
5. Solvent selection
- At matched $T_r$ and $P_r$, seven unrelated fluids agree within a factor of 2.3 in viscosity and 2.2 in $Pr$. Corresponding states works.
- Water is the outlier in thermal conductivity (5× CO₂), because hydrogen bonding survives above $T_c$.
- Selection is decided by absolute temperature and pressure, safety, regulation, cost and data availability — never by transport properties. Transport properties decide how well a feasible choice performs.
With Chapter 6 supplying the equilibrium state, Chapter 7 the computational tools, and this chapter the rates, the quantitative half of the series is complete: you can now compute where a supercritical process sits, how fast it will run, and which fluid to run it in.
Exercises
Exercise 1: Extend the Lucas Correlation to High Pressure
The Lucas method has a high-pressure extension that multiplies the dilute-gas result by a correction factor built from $T_r$ and $P_r$. Look it up in Poling et al., implement it, and compare against CoolProp for CO₂ at 40 °C over 8-30 MPa. Where does it break down, and does the failure correlate with $\rho/\rho_c$ or with $P_r$?
Exercise 2: Residual Viscosity for a Polar Fluid
Repeat Example 2 for water and for ethanol. Does the residual viscosity collapse onto a density curve as cleanly as it does for CO₂? Quantify the spread at $\rho = \rho_c$ and $\rho = 2\rho_c$, and explain any difference in terms of hydrogen bonding.
Exercise 3: The Cost of a Co-Solvent
Adding 10 mol% ethanol to CO₂ raises the mixture viscosity. Using the CoolProp HEOS backend where it works (and documenting where it does not — see Chapter 7), estimate the viscosity penalty and the resulting reduction in Stokes-Einstein diffusivity. Express the result as the percentage increase in extraction time needed to offset it, assuming intraparticle control.
Exercise 4: Locate the Enhancement Boundary
Example 4 finds a local thermal-conductivity peak up to 35 °C and none at 40 °C. Bisect on temperature to find, to 0.1 °C, the isotherm at which the local maximum disappears. Convert that temperature to $T_r$ and compare it with the $T_r$ at which the enhancement factor at $\rho_c$ falls below 4.
Exercise 5: Deteriorated Heat Transfer
Compute $Nu$ from the Dittus-Boelter correlation $Nu = 0.023 Re^{0.8} Pr^{0.4}$ for CO₂ in a 2 mm tube at 8 MPa, sweeping the bulk temperature through the pseudo-critical point at a fixed mass flux of 500 kg/(m²·s). Plot the heat-transfer coefficient against temperature. Where does the correlation predict a maximum, and why should you not trust it there?
Exercise 6: Optimal Particle Size
Combine Examples 6 and 7 into a single optimisation: minimise total extraction time subject to a pumping-power limit, with particle diameter as the decision variable. Use the intraparticle time constant as the objective and the Ergun pressure drop as the constraint. How sensitive is the optimum to the assumed tortuosity?
Exercise 7: Add a Fluid to the Survey
Extend Example 8 with ethane, ammonia and methanol. Do they fall inside the factor-of-2.3 viscosity band found for the original seven? Ammonia is strongly hydrogen-bonding — does it join water as a thermal-conductivity outlier, and by how much?
References and Further Reading
Transport Property Correlations
- Poling, Prausnitz & O'Connell, The Properties of Gases and Liquids, 5th ed. - Chapters 9-11 cover viscosity, thermal conductivity and diffusion, including the Lucas and Wilke-Chang methods used here
- Bird, Stewart & Lightfoot, Transport Phenomena, 2nd ed. - the constitutive laws and the dimensionless-group framework
- Wakao & Kaguei, Heat and Mass Transfer in Packed Beds - the packed-bed Sherwood correlation and its validity range
Supercritical Transport Properties Specifically
- Sengers & Sengers, Thermodynamic Behavior of Fluids near the Critical Point - the critical enhancement formalism behind $\lambda_c$
- Vesovic et al. and Huber et al. - the reference correlations for CO₂ viscosity and thermal conductivity that CoolProp implements
- Funazukuri, Kong & Kagei - experimental binary diffusion coefficients in supercritical CO₂
- Catchpole & King; He & Yu - diffusivity correlations fitted specifically to supercritical CO₂
Tools and Data
- CoolProp - reference transport properties via the
VandLoutputs - NIST Chemistry WebBook, Thermophysical Properties of Fluid Systems
- NIST REFPROP - the reference implementation these correlations come from