04 — Divertor Heat Flux & Blanket Neutronics

This tutorial demonstrates two nuclear engineering modules:

1. Divertor Thermal Simulation — heat-flux spreading on the divertor target
2. Blanket Neutronics — tritium breeding ratio (TBR) calculation

Both are critical for DEMO-class reactor feasibility.

License: © 1998–2026 Miroslav Šotek. GNU AGPL v3.

---

import numpy as np
import matplotlib.pyplot as plt

Part A: Divertor Thermal Simulation

The divertor receives the exhaust heat from the plasma scrape-off layer. ITER's divertor must handle ~10 MW/m² steady-state, with transients up to 20 MW/m².

from scpn_fusion.core.divertor_thermal_sim import DivertorThermalSim

# ITER-like divertor parameters
sim = DivertorThermalSim(
    target_length_m=1.2,
    target_width_m=0.3,
    n_points=100,
    coolant_temp_K=423.0,  # 150°C water
)

# Sweep over power levels
powers = np.linspace(5.0, 25.0, 20)  # MW/m²
peak_temps = []
for q in powers:
    result = sim.compute_steady_state(heat_flux_mw_m2=q)
    peak_temps.append(result["peak_temp_K"])

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(powers, np.array(peak_temps) - 273.15, "b-o", markersize=4)
ax.axhline(y=1800, color="r", linestyle="--", label="Tungsten recrystallization (1800°C)")
ax.axhline(y=3422, color="darkred", linestyle="--", label="Tungsten melting (3422°C)")
ax.set_xlabel("Heat Flux [MW/m²]")
ax.set_ylabel("Peak Surface Temperature [°C]")
ax.set_title("Divertor Target Temperature vs Heat Flux")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"At 10 MW/m²: {peak_temps[5] - 273.15:.0f}°C")
print(f"At 20 MW/m²: {peak_temps[15] - 273.15:.0f}°C")

Part B: Blanket Neutronics (Tritium Breeding)

A D-T fusion reactor must breed its own tritium. The Tritium Breeding Ratio (TBR) must exceed 1.0 for fuel self-sufficiency. ITER targets TBR ≈ 1.10.

from scpn_fusion.nuclear.blanket_neutronics import BlanketNeutronics

neutronics = BlanketNeutronics()

# Compare blanket concepts
concepts = {
    "HCLL (Li-Pb)": {"li6_enrichment": 0.30, "blanket_thickness_m": 0.8},
    "HCPB (Li-ceramic)": {"li6_enrichment": 0.60, "blanket_thickness_m": 0.7},
    "WCLL (Water-cooled)": {"li6_enrichment": 0.40, "blanket_thickness_m": 0.75},
}

results = {}
for name, params in concepts.items():
    tbr = neutronics.compute_tbr(**params)
    results[name] = tbr
    print(f"{name:25s}: TBR = {tbr:.3f} {'✓' if tbr > 1.0 else '✗'}")

# Plot
fig, ax = plt.subplots(figsize=(8, 4))
bars = ax.barh(
    list(results.keys()), list(results.values()), color=["#2196F3", "#4CAF50", "#FF9800"]
)
ax.axvline(x=1.0, color="r", linestyle="--", linewidth=2, label="TBR = 1.0 (self-sufficiency)")
ax.set_xlabel("Tritium Breeding Ratio")
ax.set_title("Blanket Concept TBR Comparison")
ax.legend()
plt.tight_layout()
plt.show()

Benchmark

%%timeit -n1 -r3
for _name, _params in concepts.items():
    neutronics.compute_tbr(**_params)

Part C: Li-6 Enrichment Sensitivity

TBR depends strongly on Li-6 enrichment. Natural lithium is only 7.5% Li-6.

enrichments = np.linspace(0.075, 0.90, 50)
tbr_values = [
    neutronics.compute_tbr(li6_enrichment=e, blanket_thickness_m=0.8) for e in enrichments
]

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(enrichments * 100, tbr_values, "g-", linewidth=2)
ax.axhline(y=1.0, color="r", linestyle="--", label="TBR = 1.0")
ax.axhline(y=1.1, color="orange", linestyle="--", alpha=0.7, label="TBR = 1.1 (target)")
ax.axvline(x=7.5, color="gray", linestyle=":", label="Natural Li-6 (7.5%)")
ax.set_xlabel("Li-6 Enrichment [%]")
ax.set_ylabel("Tritium Breeding Ratio")
ax.set_title("TBR vs Li-6 Enrichment (0.8m blanket)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Performance Benchmarks

Timing the key computations in this notebook:

1. Single divertor steady-state solve (compute_steady_state)
2. Heat flux sweep (20-point parameter scan)
3. Single TBR calculation (compute_tbr)
4. Li-6 enrichment sweep (50-point parameter scan)

import timeit


# 1. Single divertor steady-state solve
def bench_divertor_single():
    sim.compute_steady_state(heat_flux_mw_m2=10.0)


t_div = timeit.repeat(bench_divertor_single, number=100, repeat=5)
print("compute_steady_state single call (100 calls):")
print(f"  Mean: {np.mean(t_div) * 1000:.1f} ms +/- {np.std(t_div) * 1000:.1f} ms")
print(f"  Per call: {np.mean(t_div) / 100 * 1000:.2f} ms")


# 2. Heat flux sweep (20-point parameter scan)
def bench_divertor_sweep():
    for q in np.linspace(5.0, 25.0, 20):
        sim.compute_steady_state(heat_flux_mw_m2=q)


t_sweep = timeit.repeat(bench_divertor_sweep, number=10, repeat=3)
print("\nDivertor 20-point heat flux sweep (10 runs):")
print(f"  Mean: {np.mean(t_sweep) * 1000:.1f} ms +/- {np.std(t_sweep) * 1000:.1f} ms")
print(f"  Per run: {np.mean(t_sweep) / 10 * 1000:.2f} ms")


# 3. Single TBR calculation
def bench_tbr_single():
    neutronics.compute_tbr(li6_enrichment=0.30, blanket_thickness_m=0.8)


t_tbr = timeit.repeat(bench_tbr_single, number=1000, repeat=5)
print("\ncompute_tbr single call (1000 calls):")
print(f"  Mean: {np.mean(t_tbr) * 1000:.1f} ms +/- {np.std(t_tbr) * 1000:.1f} ms")
print(f"  Per call: {np.mean(t_tbr) / 1000 * 1e6:.2f} us")


# 4. Li-6 enrichment sweep (50-point parameter scan)
def bench_tbr_sweep():
    for e in np.linspace(0.075, 0.90, 50):
        neutronics.compute_tbr(li6_enrichment=e, blanket_thickness_m=0.8)


t_tbr_sweep = timeit.repeat(bench_tbr_sweep, number=10, repeat=3)
print("\nTBR 50-point enrichment sweep (10 runs):")
print(f"  Mean: {np.mean(t_tbr_sweep) * 1000:.1f} ms +/- {np.std(t_tbr_sweep) * 1000:.1f} ms")
print(f"  Per run: {np.mean(t_tbr_sweep) / 10 * 1000:.2f} ms")

Summary

• The divertor model computes 1D/2D temperature profiles under steady-state and transient heat loads
• The neutronics model estimates TBR from Li-6 enrichment, blanket geometry, and neutron multiplication
• Both modules feed into the compact reactor optimizer (see notebook 01) for integrated design
• With Rust acceleration, these calculations can be embedded in real-time digital twin scenarios