MHD Instabilities: Sawteeth and NTMs¶
This tutorial covers the two resistive MHD instabilities that most directly affect tokamak operation: sawteeth (core relaxation at the \(q = 1\) surface) and neoclassical tearing modes (NTMs, magnetic islands at \(q = 3/2\) and \(q = 2\)). You will learn how to:
Detect sawtooth crashes using the Porcelli trigger model
Apply the Kadomtsev reconnection model
Evolve NTM island width via the Modified Rutherford Equation
Stabilise NTMs with targeted ECCD
Prerequisites: Current Profile Evolution & Steady State.
Part I: Sawteeth¶
Physics¶
When the safety factor drops below 1 in the plasma core, the internal kink mode (m=1, n=1) becomes unstable. The resulting sawtooth cycle consists of:
Ramp phase (slow, \(\sim 10\)–\(100\) ms): Temperature and density rise in the core as heating exceeds transport.
Crash (fast, \(\sim 100\) \(\mu\)s): Magnetic reconnection at \(q = 1\) expels core energy, flattening profiles out to the mixing radius \(r_\text{mix}\).
Sawteeth are generally benign (they flush helium ash from the core) but large crashes can seed NTMs by generating magnetic perturbations at the \(q = 3/2\) and \(q = 2\) surfaces.
Sawtooth Simulation¶
import numpy as np
from scpn_fusion.core.sawtooth import (
SawtoothCycler,
SawtoothMonitor,
kadomtsev_crash,
)
nr = 100
rho = np.linspace(0, 1, nr)
# Initial profiles
Te = 15.0 * (1 - 0.9 * rho**2) # core-peaked temperature [keV]
q = 0.85 + 1.5 * rho**2 # q(0) = 0.85, so q=1 surface exists
# Sawtooth monitor: checks Porcelli trigger each call
monitor = SawtoothMonitor(s_crit=0.3, rho_star_crit=0.05)
# Check trigger
triggered, info = monitor.check_trigger(rho, q, Te, ne_keV=10.0, B0=5.3)
print(f"Sawtooth triggered: {triggered}")
if triggered:
print(f" q=1 location: rho = {info['rho_q1']:.3f}")
print(f" Mixing radius: rho = {info['rho_mix']:.3f}")
# Apply Kadomtsev crash
Te_post = kadomtsev_crash(rho, Te, info['rho_q1'], info['rho_mix'])
print(f" Te(0): {Te[0]:.1f} -> {Te_post[0]:.1f} keV")
Full Sawtooth Cycling¶
cycler = SawtoothCycler(
rho=rho, R0=6.2, a=2.0, B0=5.3,
s_crit=0.3, rho_star_crit=0.05,
)
# Simulate 50 ms of plasma evolution with sawteeth
dt = 1e-3 # 1 ms timestep
Te_current = Te.copy()
q_current = q.copy()
crash_times = []
for step in range(50):
t = step * dt
# Artificial core heating (ramp phase)
Te_current[:nr//3] += 0.1 # 0.1 keV per ms in the core
crashed, Te_current, info = cycler.step(
Te_current, q_current, ne_keV=10.0, dt=dt,
)
if crashed:
crash_times.append(t)
print(f" Crash at t = {t*1e3:.0f} ms")
print(f"\n{len(crash_times)} sawtooth crashes in 50 ms")
if len(crash_times) > 1:
period = np.diff(crash_times).mean() * 1e3
print(f"Average period: {period:.1f} ms")
Part II: Neoclassical Tearing Modes¶
Physics¶
NTMs are magnetic islands at rational surfaces (\(q = m/n\)) driven by the loss of bootstrap current inside the island. The island creates a flat pressure region (no gradient, no bootstrap current), reducing the total current and further destabilising the island through a positive feedback loop.
The island half-width \(w\) evolves according to the Modified Rutherford Equation (see (2) in the textbook):
Key terms:
\(\Delta' < 0\) (classically stable) — tearing is subcritical
\(\alpha_\text{bs}\) (bootstrap drive) — makes the mode grow
\(\alpha_\text{GGJ}\) (curvature) — stabilises small islands
\(\alpha_\text{pol}\) (polarisation) — provides a seed threshold
\(\alpha_\text{ECCD}\) — targeted ECCD adds stabilising current
NTM Simulation¶
from scpn_fusion.core.ntm_dynamics import (
NTMIslandDynamics,
NTMController,
find_rational_surfaces,
)
# Find rational surfaces in the q-profile
q_profile = 0.85 + 1.5 * rho**2 + 0.5 * rho**4
surfaces = find_rational_surfaces(rho, q_profile, [1.0, 1.5, 2.0])
for q_val, rho_loc in surfaces.items():
if rho_loc is not None:
print(f"q = {q_val} at rho = {rho_loc:.3f}")
# Initialise NTM at q = 2 surface
ntm = NTMIslandDynamics(
m=2, n=1,
r_s=surfaces[2.0] * 2.0, # dimensional radius [m]
a=2.0,
R0=6.2,
B0=5.3,
tau_R=300.0, # resistive time [s]
)
# Seed the island (e.g. from a sawtooth crash)
ntm.w = 0.02 # 2 cm seed island
# Set plasma parameters at the rational surface
ntm.j_bs = 0.5e6 # bootstrap current density [A/m^2]
ntm.j_phi = 1.0e6 # total current density [A/m^2]
# Evolve without ECCD (island grows)
dt = 0.1
for step in range(200):
ntm.step(dt, j_eccd=0.0)
print(f"Island width after 20 s (no ECCD): {ntm.w*100:.1f} cm")
print(f"Saturated: {ntm.is_saturated()}")
ECCD Stabilisation¶
# Reset and apply ECCD
ntm.w = 0.02
controller = NTMController(
target_modes=[(2, 1)],
eccd_efficiency=0.25,
eccd_power_max_mw=5.0,
)
w_history = []
for step in range(500):
# Controller decides ECCD power based on island size
j_eccd = controller.compute_eccd(ntm)
ntm.step(dt, j_eccd=j_eccd)
w_history.append(ntm.w * 100) # cm
print(f"Island width after 50 s (with ECCD): {ntm.w*100:.2f} cm")
Plotting the Result¶
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(np.arange(len(w_history)) * dt, w_history)
ax.set_xlabel("Time [s]")
ax.set_ylabel("Island Width [cm]")
ax.set_title("NTM (2,1) Stabilisation with ECCD")
ax.axhline(2.0, color="red", ls="--", label="Seed width")
ax.legend()
plt.tight_layout()
plt.show()
Design Guidelines¶
Sawtooth control:
Longer sawtooth periods → larger crashes → bigger NTM seeds. ECCD at the \(q = 1\) surface (pacing) can shorten the period.
Ion cyclotron resonance heating (ICRH) with minority ions stabilises sawteeth by enhancing fast-ion pressure — but delays the crash, potentially creating a monster sawtooth.
NTM prevention:
Keep \(q_\text{min} > 1.5\) (no \(q = 3/2\) surface).
Maintain ECCD standby at the \(q = 2\) surface.
Monitor island width via locked-mode detection coils; trigger ECCD above a threshold island size.
ECCD alignment:
The ECCD deposition must be localised within \(\pm 2\) cm of the rational surface. Real-time EFIT (see Real-Time Equilibrium Reconstruction & Shape Control) provides the equilibrium update needed to steer the ECCD mirrors.