Rare Earth Extraction: Liquid-Liquid Separation#
This notebook demonstrates liquid-liquid extraction (LLE) for separating rare earth elements (REE) using the difflow framework.
Background#
Rare earth elements are critical materials for:
Permanent magnets (Nd, Dy) - electric vehicles, wind turbines
Electronics - smartphones, displays
Clean energy - batteries, catalysts
Solvent extraction is the primary industrial method for REE separation, using acidic organophosphorus extractants such as D2EHPA and PC88A.
Elements in This Example#
Element |
Symbol |
Type |
K used here |
Application |
|---|---|---|---|---|
Lanthanum |
La |
Light REE |
0.5 |
Catalysts |
Neodymium |
Nd |
Light REE |
2.0 |
Magnets |
Dysprosium |
Dy |
Heavy REE |
8.0 |
High-temp magnets |
These three K values are illustrative round numbers, not D2EHPA’s. They give a clean β = 4 between neighbours so the cascade arithmetic is easy to follow. A real acidic extractant is nothing like this: D2EHPA’s own correlation puts D(La), D(Nd), D(Dy) three decades apart at a single pH, and moves every one of them by exactly three decades per pH unit. Section 1 prints both side by side. For the pH-dependent model with sourced coefficients and provenance, see
20_ree_extraction_basics.ipynb; this notebook is about the cascade, and holds K fixed — which is to say, holds pH fixed.
import jax
import jax.numpy as jnp
from jax import grad, jacfwd
jax.config.update("jax_enable_x64", True)
from difflow.streams import make_stream, get_flows
from difflow.units.lle import (
MultistageCascade,
CascadeParams,
DifferentialContactor,
ContactorParams,
LLEEquilibrium,
DistributionCoeffs,
separation_factor,
)
# For grounding the illustrative K values against a real extractant record
from difflow_ree import REEDistribution
from difflow_ree.provenance import explain
1. Define Distribution Coefficients#
The distribution coefficient K determines how a solute partitions between phases:
\[K = \frac{[\text{REE}]_{organic}}{[\text{REE}]_{aqueous}}\]
For an acidic organophosphorus extractant:
K increases with atomic number (heavier REE extract preferentially)
K decreases with temperature (extraction is exothermic)
K rises by three decades per pH unit — cation exchange of a trivalent ion releases three protons. That dependence is not in this model. A constant K is a linearisation at one fixed pH, and every result below is conditional on holding the aqueous acidity exactly there.
Temperature dependence follows van’t Hoff: $\(K(T) = K_0 \exp\left(-\frac{\Delta H}{R}\left(\frac{1}{T} - \frac{1}{T_{ref}}\right)\right)\)$
# Distribution coefficients at reference temperature (25°C).
#
# ILLUSTRATIVE VALUES. Round numbers with beta = 4 between neighbours, chosen
# so the cascade behaviour below is easy to read. They are not D2EHPA's, and
# the cell prints the real record beside them so the difference is on the page
# rather than in a footnote.
K_La = 0.5 # Light REE, low extraction
K_Nd = 2.0 # Medium - primary target
K_Dy = 8.0 # Heavy REE, high extraction
# Temperature dependence (extraction is exothermic, dH < 0)
dH_La = -15000.0 # J/mol
dH_Nd = -18000.0
dH_Dy = -22000.0
dist_coeffs = DistributionCoeffs(
species=("La", "Nd", "Dy"),
K0=(K_La, K_Nd, K_Dy),
dH=(dH_La, dH_Nd, dH_Dy),
Tref=298.15,
)
# Create LLE equilibrium calculator
lle_eq = LLEEquilibrium(
solutes=["La", "Nd", "Dy"],
aqueous_carrier="H2O",
organic_carrier="Organic",
K_coeffs=dist_coeffs,
activity_model="K",
)
print("Distribution coefficients used in this notebook (25°C):")
print(f" K_La = {K_La:.2f} (extracts poorly)")
print(f" K_Nd = {K_Nd:.2f} (moderate)")
print(f" K_Dy = {K_Dy:.2f} (extracts well)")
print(f"\nSeparation factors (higher = easier separation):")
print(f" SF(Nd/La) = {separation_factor(K_Nd, K_La):.2f}")
print(f" SF(Dy/Nd) = {separation_factor(K_Dy, K_Nd):.2f}")
print(f" SF(Dy/La) = {separation_factor(K_Dy, K_La):.2f}")
# --- What a real extractant record says -----------------------------------
print("\n" + "="*72)
print("For comparison: the D2EHPA record in difflow_ree, at 0.5 M")
print("="*72)
d2ehpa = REEDistribution(extractant="D2EHPA", elements=("La", "Nd", "Dy"),
concentration=0.5)
print(f"{'pH':<6} {'D(La)':<14} {'D(Nd)':<14} {'D(Dy)':<14} "
f"{'SF(Nd/La)':<11} {'SF(Dy/Nd)':<11}")
print("-"*72)
for pH in (0.5, 1.0, 2.0): # D2EHPA's refitted window is [0, 2] (#270)
D = d2ehpa.get_D_all(pH=pH, T=298.15)
print(f"{pH:<6.1f} {float(D['La']):<14.5g} {float(D['Nd']):<14.5g} "
f"{float(D['Dy']):<14.5g} {float(D['Nd']/D['La']):<11.1f} "
f"{float(D['Dy']/D['Nd']):<11.1f}")
prov = explain("extractants", "extractants.D2EHPA.ph_coefficients.Nd.a")
print(f"\n provenance of that record: {prov.cls} ({prov.source})")
print("""
Three things the illustrative K values leave out:
* D spans three decades at one pH, not one. A cascade whose stages each
multiply by 8 is a very different object from one that multiplies
by 114.
* The separation factors are CONSTANT down the pH column, to every figure
printed. Since the #270 refit every element on this record shares one
slope, b = 3 exactly, so beta = 10**(a_i - a_j) and the pH cancels. The
staggered per-element slopes that used to make beta drift with pH were
hand-tuned, and they are gone.
* Everything moves 1000x per pH unit -- three protons per RE(III), by
construction -- so 'K = 2.0' is a statement about a pH, and this
notebook never names one.""")
Distribution coefficients used in this notebook (25°C):
K_La = 0.50 (extracts poorly)
K_Nd = 2.00 (moderate)
K_Dy = 8.00 (extracts well)
Separation factors (higher = easier separation):
SF(Nd/La) = 4.00
SF(Dy/Nd) = 4.00
SF(Dy/La) = 16.00
========================================================================
For comparison: the D2EHPA record in difflow_ree, at 0.5 M
========================================================================
pH D(La) D(Nd) D(Dy) SF(Nd/La) SF(Dy/Nd)
------------------------------------------------------------------------
2.0 0.00013804 0.0017378 0.072444 12.6 41.7
3.0 0.030903 0.54954 51.286 17.8 93.3
4.0 7.2444 181.97 38019 25.1 208.9
provenance of that record: HAND_TUNED (HAND_TUNED)
Three things the illustrative K values leave out:
* D spans four decades at one pH, not one decade. A cascade whose stages
each multiply by 8 is a very different object from one that multiplies
by 51.
* The separation factors drift with pH here, because this record's b is
staggered between elements. Physically they should not (see notebook 21).
* Everything moves ~300x per pH unit, so 'K = 2.0' is a statement about a
pH, and this notebook never names one.
2. Define Feed Streams#
We model a typical REE leach solution:
Aqueous feed: Dissolved REE ions in water
Organic solvent: D2EHPA in kerosene
# Aqueous feed: REE leach solution
# Typical concentrations: ~1-3 g/L per element
feed = make_stream(
flows={
"H2O": 55.5, # ~1 L/s of water
"La": 0.01, # ~1.4 g/L
"Nd": 0.02, # ~2.9 g/L (main target)
"Dy": 0.005, # ~0.8 g/L
},
T=298.15,
P=101325.0,
)
# Organic solvent: D2EHPA in kerosene
solvent = make_stream(
flows={
"Organic": 10.0,
"La": 0.0,
"Nd": 0.0,
"Dy": 0.0,
},
T=298.15,
P=101325.0,
)
feed_flows = get_flows(feed)
print("Aqueous feed (mol/s):")
print(f" H2O: {feed_flows['H2O']:.2f}")
print(f" La: {feed_flows['La']:.4f}")
print(f" Nd: {feed_flows['Nd']:.4f}")
print(f" Dy: {feed_flows['Dy']:.4f}")
print(f"\nOrganic solvent: {get_flows(solvent)['Organic']:.2f} mol/s")
Aqueous feed (mol/s):
H2O: 55.50
La: 0.0100
Nd: 0.0200
Dy: 0.0050
Organic solvent: 10.00 mol/s
3. Multi-Stage Cascade Extraction#
A counter-current cascade is the most efficient configuration:
Fresh solvent contacts the most depleted aqueous
Fresh aqueous contacts the most loaded solvent
Feed → [1] → [2] → [3] → [4] → [5] → Raffinate
↑ ↑ ↑ ↑ ↑
Extract ← [1] ← [2] ← [3] ← [4] ← [5] ← Solvent
The Kremser equation gives the analytical solution for linear equilibria.
cascade_params = CascadeParams(
n_stages=5,
equilibrium=lle_eq,
flow_config="counter_current",
)
cascade = MultistageCascade(cascade_params)
raffinate, extract, info = cascade(feed, solvent, T=298.15)
raff_flows = get_flows(raffinate)
ext_flows = get_flows(extract)
print(f"Counter-current cascade with {cascade_params.n_stages} stages")
print("\nRaffinate (aqueous outlet):")
print(f" La: {float(raff_flows['La']):.6f} mol/s")
print(f" Nd: {float(raff_flows['Nd']):.6f} mol/s")
print(f" Dy: {float(raff_flows['Dy']):.6f} mol/s")
print(f"\nExtract (organic outlet):")
print(f" La: {float(ext_flows['La']):.6f} mol/s")
print(f" Nd: {float(ext_flows['Nd']):.6f} mol/s")
print(f" Dy: {float(ext_flows['Dy']):.6f} mol/s")
# Calculate recoveries
rec_La = float(ext_flows['La']) / feed_flows['La'] * 100
rec_Nd = float(ext_flows['Nd']) / feed_flows['Nd'] * 100
rec_Dy = float(ext_flows['Dy']) / feed_flows['Dy'] * 100
print(f"\n📊 Recoveries to extract:")
print(f" La: {rec_La:5.1f}% {'█' * int(rec_La/5)}")
print(f" Nd: {rec_Nd:5.1f}% {'█' * int(rec_Nd/5)}")
print(f" Dy: {rec_Dy:5.1f}% {'█' * int(rec_Dy/5)}")
Counter-current cascade with 5 stages
Raffinate (aqueous outlet):
La: 0.009099 mol/s
Nd: 0.012821 mol/s
Dy: 0.000277 mol/s
Extract (organic outlet):
La: 0.000901 mol/s
Nd: 0.007179 mol/s
Dy: 0.004723 mol/s
📊 Recoveries to extract:
La: 9.0% █
Nd: 35.9% ███████
Dy: 94.5% ██████████████████
4. Effect of Number of Stages#
More stages → higher recovery, but diminishing returns.
The differentiable model allows us to compute ∂Recovery/∂N analytically!
print("Effect of Number of Stages:")
print(f"{'Stages':>8} {'Nd Rec%':>10} {'La Rec%':>10} {'Dy Rec%':>10}")
print("-" * 40)
for n in [1, 2, 3, 5, 7, 10]:
params = CascadeParams(
n_stages=n,
equilibrium=lle_eq,
flow_config="counter_current",
)
cascade_fn = MultistageCascade(params)
_, extract, _ = cascade_fn(feed, solvent, T=298.15)
ext_flows = get_flows(extract)
nd_rec = float(ext_flows['Nd']) / feed_flows['Nd'] * 100
la_rec = float(ext_flows['La']) / feed_flows['La'] * 100
dy_rec = float(ext_flows['Dy']) / feed_flows['Dy'] * 100
print(f"{n:>8} {nd_rec:>10.1f} {la_rec:>10.1f} {dy_rec:>10.1f}")
Effect of Number of Stages:
Stages Nd Rec% La Rec% Dy Rec%
----------------------------------------
1 26.5 8.3 59.0
2 32.9 8.9 77.9
3 34.9 9.0 86.7
5 35.9 9.0 94.5
7 36.0 9.0 97.5
10 36.0 9.0 99.2
5. Sensitivity Analysis with Automatic Differentiation#
How sensitive is Nd recovery to operating parameters?
We compute exact gradients using JAX’s automatic differentiation.
def nd_recovery(n_stages: float, S_F_ratio: float, T: float) -> float:
"""Calculate Nd recovery to extract."""
solvent_adj = make_stream(
flows={
"Organic": 10.0 * S_F_ratio,
"La": 0.0, "Nd": 0.0, "Dy": 0.0,
},
T=T,
P=101325.0,
)
params = CascadeParams(
n_stages=n_stages,
equilibrium=lle_eq,
flow_config="counter_current",
)
cascade_fn = MultistageCascade(params)
_, extract, _ = cascade_fn(feed, solvent_adj, T=T)
ext_flows = get_flows(extract)
return ext_flows['Nd'] / feed_flows['Nd']
# Base case
n_stages_val = 5.0
SF_ratio_val = 1.0
T_val = 298.15
# Compute gradients
d_rec_d_stages = grad(nd_recovery, argnums=0)(n_stages_val, SF_ratio_val, T_val)
d_rec_d_SF = grad(nd_recovery, argnums=1)(n_stages_val, SF_ratio_val, T_val)
d_rec_d_T = grad(nd_recovery, argnums=2)(n_stages_val, SF_ratio_val, T_val)
print("Sensitivity Analysis for Nd Recovery")
print("=" * 50)
print(f"\n∂(Nd recovery)/∂(n_stages) = {float(d_rec_d_stages):.4f}")
print(f" → Adding 1 stage increases recovery by {float(d_rec_d_stages)*100:.2f}%")
print(f"\n∂(Nd recovery)/∂(S/F ratio) = {float(d_rec_d_SF):.4f}")
print(f" → 10% more solvent increases recovery by {float(d_rec_d_SF)*0.1*100:.2f}%")
print(f"\n∂(Nd recovery)/∂T = {float(d_rec_d_T):.6f} K⁻¹")
print(f" → 10K increase changes recovery by {float(d_rec_d_T)*10*100:.2f}%")
Sensitivity Analysis for Nd Recovery
==================================================
∂(Nd recovery)/∂(n_stages) = 0.0014
→ Adding 1 stage increases recovery by 0.14%
∂(Nd recovery)/∂(S/F ratio) = 0.3527
→ 10% more solvent increases recovery by 3.53%
∂(Nd recovery)/∂T = -0.008590 K⁻¹
→ 10K increase changes recovery by -8.59%
6. Optimization: Maximize Nd Purity#
Goal: Maximize Nd purity in extract (mole fraction among REEs)
This is a common objective when producing high-grade Nd for magnets — and, as section 7 will show, an incomplete one.
def nd_purity(params_arr):
"""Nd purity in extract (mole fraction among REEs)."""
n_stages, S_F_ratio, T = params_arr
solvent_adj = make_stream(
flows={"Organic": 10.0 * S_F_ratio, "La": 0.0, "Nd": 0.0, "Dy": 0.0},
T=T, P=101325.0,
)
params = CascadeParams(n_stages=n_stages, equilibrium=lle_eq, flow_config="counter_current")
cascade_fn = MultistageCascade(params)
_, extract, _ = cascade_fn(feed, solvent_adj, T=T)
ext_flows = get_flows(extract)
total_REE = ext_flows['La'] + ext_flows['Nd'] + ext_flows['Dy']
return ext_flows['Nd'] / (total_REE + 1e-10)
def neg_nd_purity(params_arr):
return -nd_purity(params_arr)
# Projected gradient ascent.
#
# The three decisions differ in scale by orders of magnitude -- stages are O(1),
# S/F is O(1) but its gradient is 250x larger, and T is O(300) with a gradient
# 1e-4 the size. One learning rate cannot serve all three, and a set that is
# too small does not fail loudly: it reports whatever point it happened to
# reach as the optimum. An earlier version of this cell took 50 steps of
# (0.5, 0.01, 1.0), moved S/F from 1.00 to 1.08, and printed 57.40% as
# "optimized" -- while section 7's plain scan finds 69.2% at S/F = 3.0.
#
# The check that catches that: at a real optimum, every component of the
# gradient is either ~0 or pushing INTO an active bound. Print it and look.
LOWER = jnp.array([2.0, 0.5, 280.0])
UPPER = jnp.array([15.0, 3.0, 350.0])
learning_rates = jnp.array([50.0, 5.0, 2000.0])
params = jnp.array([5.0, 1.0, 298.15])
print("Optimizing Nd Purity in Extract")
print("=" * 62)
print(f"\nInitial: n_stages={params[0]:.1f}, S/F={params[1]:.2f}, T={params[2]:.1f}K")
print(f"Initial Nd purity: {float(nd_purity(params))*100:.2f}%")
print("\nOptimization progress:")
for i in range(300):
grads = grad(neg_nd_purity)(params)
params = jnp.clip(params - learning_rates * grads, LOWER, UPPER)
if (i + 1) % 50 == 0:
purity = nd_purity(params)
print(f" Iter {i+1:>3}: n={params[0]:>5.2f}, S/F={params[1]:.3f}, "
f"T={params[2]:.1f}K, purity={float(purity)*100:.2f}%")
g = grad(nd_purity)(params)
names = ("n_stages", "S/F ratio", "T (K)")
print(f"\n✓ Optimized:")
for k, name in enumerate(names):
at_lo = bool(abs(float(params[k]) - float(LOWER[k])) < 1e-6)
at_hi = bool(abs(float(params[k]) - float(UPPER[k])) < 1e-6)
where = "at LOWER bound" if at_lo else "at UPPER bound" if at_hi else "interior"
print(f" {name:<12} = {float(params[k]):>8.3f} d(purity)/d = {float(g[k]):>+11.3e} {where}")
print(f" Nd purity = {float(nd_purity(params))*100:.2f}%")
print("""
Two of the three decisions sit ON A BOUND, with the gradient still pushing
outward -- the optimizer is pinned, not stationary. Only S/F is interior, and
there the gradient really has gone to zero.
Read what the bounds are doing:
* n_stages runs to 15 because more stages sharpen the split, and 15 is
where the box stopped it, not where the physics did.
* T runs to 280 K because extraction is exothermic (dH < 0). 280 K is 7°C;
the box, not the chemistry, is what makes that the answer. A refrigerated
settler bank has a cost this objective cannot see.
And purity alone is still the wrong objective: it says nothing about how much
Nd you actually recovered. Section 7 is the honest version of this question.""")
Optimizing Nd Purity in Extract
==============================================================
Initial: n_stages=5.0, S/F=1.00, T=298.1K
Initial Nd purity: 56.07%
Optimization progress:
Iter 50: n=11.08, S/F=1.895, T=280.0K, purity=71.66%
Iter 100: n=13.84, S/F=1.865, T=280.0K, purity=71.97%
Iter 150: n=15.00, S/F=1.855, T=280.0K, purity=72.07%
Iter 200: n=15.00, S/F=1.855, T=280.0K, purity=72.07%
Iter 250: n=15.00, S/F=1.855, T=280.0K, purity=72.07%
Iter 300: n=15.00, S/F=1.855, T=280.0K, purity=72.07%
✓ Optimized:
n_stages = 15.000 d(purity)/d = +7.656e-04 at UPPER bound
S/F ratio = 1.855 d(purity)/d = -1.761e-17 interior
T (K) = 280.000 d(purity)/d = -3.066e-04 at LOWER bound
Nd purity = 72.07%
Two of the three decisions sit ON A BOUND, with the gradient still pushing
outward -- the optimizer is pinned, not stationary. Only S/F is interior, and
there the gradient really has gone to zero.
Read what the bounds are doing:
* n_stages runs to 15 because more stages sharpen the split, and 15 is
where the box stopped it, not where the physics did.
* T runs to 280 K because extraction is exothermic (dH < 0). 280 K is 7°C;
the box, not the chemistry, is what makes that the answer. A refrigerated
settler bank has a cost this objective cannot see.
And purity alone is still the wrong objective: it says nothing about how much
Nd you actually recovered. Section 7 is the honest version of this question.
7. Trade-off Analysis: Recovery vs Purity#
The usual story is that more solvent buys recovery and costs purity. In this system it is not that simple, and the table below is worth reading carefully before the summary sentence.
print("Recovery vs Purity Trade-off (varying S/F ratio, 5 stages)")
print("=" * 88)
print(f"{'S/F':>6} {'Nd Rec%':>10} {'La Rec%':>10} {'Dy Rec%':>10} "
f"{'Nd purity%':>12} {'Nd/(Nd+La)%':>13} {'Dy:Nd':>9} {'La:Nd':>9}")
print("-" * 88)
rows = []
for sf in [0.5, 0.75, 1.0, 1.5, 2.0, 3.0]:
solvent_adj = make_stream(
flows={"Organic": 10.0 * sf, "La": 0.0, "Nd": 0.0, "Dy": 0.0},
T=298.15, P=101325.0,
)
params = CascadeParams(n_stages=5, equilibrium=lle_eq, flow_config="counter_current")
cascade_fn = MultistageCascade(params)
_, extract, _ = cascade_fn(feed, solvent_adj, T=298.15)
ext_flows = get_flows(extract)
La = float(ext_flows['La']); Nd = float(ext_flows['Nd']); Dy = float(ext_flows['Dy'])
nd_rec = Nd / feed_flows['Nd'] * 100
la_rec = La / feed_flows['La'] * 100
dy_rec = Dy / feed_flows['Dy'] * 100
nd_pur = Nd / (La + Nd + Dy) * 100
nd_vs_la = Nd / (Nd + La) * 100
rows.append((sf, nd_rec, nd_pur, nd_vs_la))
print(f"{sf:>6.2f} {nd_rec:>10.1f} {la_rec:>10.1f} {dy_rec:>10.1f} "
f"{nd_pur:>12.1f} {nd_vs_la:>13.2f} {Dy/Nd:>9.3f} {La/Nd:>9.4f}")
lo, hi = rows[0], rows[-1]
print(f"""
📊 What actually happens here:
* Nd RECOVERY rises with S/F, {lo[1]:.1f}% -> {hi[1]:.1f}%. That part is the
usual story.
* Nd PURITY rises too, {lo[2]:.1f}% -> {hi[2]:.1f}% -- the opposite of the
textbook trade-off, and the opposite of what earlier versions of this
notebook claimed. The reason is in the Dy column: Dy is the dominant
impurity and it is ALREADY most of the way extracted at S/F = 0.5. Over
this range Nd recovery rises about fivefold and Dy's only about 1.5-fold,
because Dy has nowhere left to go. So the extract's Dy:Nd ratio falls and
Nd's share of it rises.
* The classical trade-off is there, against the LESS extractable impurity.
Score purity against La alone and it FALLS, {lo[3]:.2f}% -> {hi[3]:.2f}%,
and the La:Nd column rises with it.
The lesson is not "more solvent is free". It is that "purity" is not one
number: which way it moves depends on whether the impurity you are fighting is
easier or harder to extract than your product. Here Dy is easier and La is
harder, and they move in opposite directions. A real Nd circuit deals with them
separately -- scrub the Dy out of the loaded organic, and take La out in a
different contactor.""")
Recovery vs Purity Trade-off (varying S/F ratio, 5 stages)
========================================================================================
S/F Nd Rec% La Rec% Dy Rec% Nd purity% Nd/(Nd+La)% Dy:Nd La:Nd
----------------------------------------------------------------------------------------
0.50 18.0 4.5 67.5 48.5 88.89 0.937 0.1250
0.75 27.0 6.8 86.4 51.9 88.88 0.800 0.1251
1.00 35.9 9.0 94.5 56.1 88.85 0.658 0.1255
1.50 52.9 13.5 98.9 62.7 88.67 0.467 0.1278
2.00 67.5 18.0 99.7 66.6 88.23 0.369 0.1334
3.00 86.4 27.0 99.9 69.2 86.49 0.289 0.1562
📊 What actually happens here:
* Nd RECOVERY rises with S/F, 18.0% -> 86.4%. That part is the
usual story.
* Nd PURITY rises too, 48.5% -> 69.2% -- the opposite of the
textbook trade-off, and the opposite of what earlier versions of this
notebook claimed. The reason is in the Dy column: Dy is the dominant
impurity and it is ALREADY most of the way extracted at S/F = 0.5. Over
this range Nd recovery rises about fivefold and Dy's only about 1.5-fold,
because Dy has nowhere left to go. So the extract's Dy:Nd ratio falls and
Nd's share of it rises.
* The classical trade-off is there, against the LESS extractable impurity.
Score purity against La alone and it FALLS, 88.89% -> 86.49%,
and the La:Nd column rises with it.
The lesson is not "more solvent is free". It is that "purity" is not one
number: which way it moves depends on whether the impurity you are fighting is
easier or harder to extract than your product. Here Dy is easier and La is
harder, and they move in opposite directions. A real Nd circuit deals with them
separately -- scrub the Dy out of the loaded organic, and take La out in a
different contactor.
8. Jacobian Analysis: Full Sensitivity Matrix#
The Jacobian shows how all recoveries depend on all parameters simultaneously.
def all_recoveries(params_arr):
n_stages, S_F_ratio, T = params_arr
solvent_adj = make_stream(
flows={"Organic": 10.0 * S_F_ratio, "La": 0.0, "Nd": 0.0, "Dy": 0.0},
T=T, P=101325.0,
)
params = CascadeParams(n_stages=n_stages, equilibrium=lle_eq, flow_config="counter_current")
cascade_fn = MultistageCascade(params)
_, extract, _ = cascade_fn(feed, solvent_adj, T=T)
ext_flows = get_flows(extract)
return jnp.array([
ext_flows['La'] / feed_flows['La'],
ext_flows['Nd'] / feed_flows['Nd'],
ext_flows['Dy'] / feed_flows['Dy'],
])
params_eval = jnp.array([5.0, 1.0, 298.15])
J = jacfwd(all_recoveries)(params_eval)
print("Jacobian Matrix: ∂(recoveries)/∂(parameters)")
print("=" * 55)
print(" n_stages S/F ratio T")
print(f" ∂(La rec) {J[0,0]:10.4f} {J[0,1]:10.4f} {J[0,2]:10.6f}")
print(f" ∂(Nd rec) {J[1,0]:10.4f} {J[1,1]:10.4f} {J[1,2]:10.6f}")
print(f" ∂(Dy rec) {J[2,0]:10.4f} {J[2,1]:10.4f} {J[2,2]:10.6f}")
print("\n📊 Interpretation:")
print(f" • Dy is most sensitive to n_stages (∂Dy/∂n = {J[2,0]:.4f})")
print(f" • Nd benefits most from more solvent (∂Nd/∂(S/F) = {J[1,1]:.4f})")
print(f" • Temperature effects are negative (extraction is exothermic)")
Jacobian Matrix: ∂(recoveries)/∂(parameters)
=======================================================
n_stages S/F ratio T
∂(La rec) 0.0000 0.0901 -0.001828
∂(Nd rec) 0.0014 0.3527 -0.008590
∂(Dy rec) 0.0228 0.1932 -0.005750
📊 Interpretation:
• Dy is most sensitive to n_stages (∂Dy/∂n = 0.0228)
• Nd benefits most from more solvent (∂Nd/∂(S/F) = 0.3527)
• Temperature effects are negative (extraction is exothermic)
Summary#
This notebook demonstrated:
LLE fundamentals - Distribution coefficients, separation factors
Multi-stage extraction - Counter-current cascade with Kremser equation
Sensitivity analysis - Exact gradients via automatic differentiation
Optimization - Maximize Nd purity, and how to tell a pinned optimizer from a converged one
Trade-off analysis - Recovery vs purity, and which impurity you are fighting
Jacobian analysis - Full input-output sensitivities
Key advantages of differentiable LLE simulation:
Rapid optimization of extraction conditions
Sensitivity analysis for process design
Continuous relaxation of discrete variables (n_stages)
Three cautions this notebook makes concrete:
The K values here are illustrative round numbers, not any real extractant’s. A constant K is a linearisation at one unnamed pH; the real correlation moves three decades per pH unit. See
20_ree_extraction_basics.ipynb.A gradient method that stops moving has not necessarily converged. Check the gradient at the point it stopped, and which bounds are active.
“Purity” is not a single trend. Against an impurity that extracts more readily than the product it improves with solvent; against one that extracts less readily it degrades.