Part 6, Chapter 1

Magnetic Confinement

Tokamak equilibrium, safety factor, stability, and the Lawson criterion

1.1 Tokamak Equilibrium and the Grad-Shafranov Equation

A tokamak confines plasma in a toroidal geometry using a combination of toroidal and poloidal magnetic fields. The equilibrium of an axisymmetric toroidal plasma is governed by force balance between the magnetic pressure, magnetic tension, and plasma pressure gradients. In cylindrical coordinates (R, phi, Z), all equilibrium quantities depend only on R and Z, and the magnetic field can be written in terms of a poloidal flux function psi(R,Z).

The toroidal component of the magnetic field is expressed as B_phi = F(psi)/R, where F(psi) encodes the poloidal current. The poloidal field components are B_R = -(1/R) partial psi / partial Z and B_Z = (1/R) partial psi / partial R. Substituting into the MHD force balance J x B = grad p yields the Grad-Shafranov equation:

$$R\frac{\partial}{\partial R}\!\left(\frac{1}{R}\frac{\partial\psi}{\partial R}\right) + \frac{\partial^2\psi}{\partial Z^2} = -\mu_0 R^2 \frac{dp}{d\psi} - \frac{1}{2}\frac{dF^2}{d\psi}$$

This is a nonlinear elliptic PDE. The two free functions p(psi) and F(psi) must be specified to close the system. The left-hand side is often written compactly as Delta* psi, the Grad-Shafranov operator. The equation states that the toroidal current density j_phi = R dp/dpsi + (1/mu_0 R) F dF/dpsi provides the source for the poloidal flux.

The simplest analytic solution is the Solov'ev equilibrium, where p'(psi) = C_1 and FF'(psi) = C_2 are taken as constants. The Grad-Shafranov equation then becomes linear and admits closed-form solutions that capture the essential features of tokamak flux surfaces, including the Shafranov shift (outward displacement of inner flux surfaces due to the toroidal geometry).

1.2 Safety Factor and Magnetic Topology

The safety factor q characterizes how many times a magnetic field line wraps around the torus toroidally for each poloidal circuit. For a large-aspect-ratio circular tokamak, it is given by:

$$q = \frac{r B_\phi}{R B_\theta}$$

More generally, q is defined as a flux-surface average:

$$q(\psi) = \frac{1}{2\pi}\oint \frac{B \cdot \nabla\phi}{B \cdot \nabla\theta}\,d\theta$$

The profile q(r) is central to stability analysis. Rational surfaces where q = m/n (integers) are locations where MHD instabilities can develop because field lines close on themselves after m toroidal and n poloidal transits, allowing perturbations to resonate. The magnetic shear s = (r/q)(dq/dr) measures how rapidly field-line pitch varies, with positive shear generally stabilizing tearing modes.

Typical tokamak profiles have q_0 approximately 1 at the magnetic axis and q_a of 3-5 at the edge. The cylindrical safety factor for a tokamak with plasma current I_p, minor radius a, and toroidal field B_0 is q_a = 2 pi a^2 B_0 / (mu_0 R_0 I_p), which shows the inverse relationship between q and plasma current.

1.3 Kink Stability: The Kruskal-Shafranov Limit

The external kink instability is one of the most dangerous MHD modes in a tokamak. It arises when the plasma column undergoes a helical deformation that grows without bound. The Kruskal-Shafranov criterion provides the fundamental stability condition:

$$q(a) > 1 \quad \text{(Kruskal-Shafranov limit)}$$

If q(a) falls below 1, the m=1, n=1 external kink mode becomes unstable, leading to a violent disruption. In practice, tokamaks operate with q(a) greater than 2 to maintain a margin against internal kink modes and other resistive instabilities. The energy principle for ideal MHD stability gives the change in potential energy:

$$\delta W = \frac{1}{2}\int \left[\frac{|\delta B|^2}{\mu_0} + \gamma p |\nabla \cdot \boldsymbol{\xi}|^2 - \boldsymbol{\xi}_\perp \cdot \nabla p\,(\nabla \cdot \boldsymbol{\xi}_\perp) - j_\parallel (\boldsymbol{\xi}_\perp \times \hat{b}) \cdot \delta B\right] d^3x$$

The last term (the kink drive) shows that parallel current j_parallel drives the instability. When delta W is less than 0 for any trial displacement xi, the equilibrium is unstable. This provides a rigorous necessary and sufficient condition for ideal MHD stability.

1.4 The Lawson Criterion and Ignition

For a magnetically confined fusion plasma to produce net energy, the fusion power must exceed all loss channels. The Lawson criterion provides the minimum confinement parameter for energy breakeven. For D-T fusion at optimal temperature (T approximately 15 keV), the fusion triple product must satisfy:

$$n\,\tau_E\,T > 3 \times 10^{21}\;\text{keV·s/m}^3$$

The derivation begins with power balance. The alpha particle heating power is:

$$P_\alpha = \frac{n_D n_T \langle\sigma v\rangle E_\alpha}{1} = \frac{n^2}{4}\langle\sigma v\rangle E_\alpha$$

where E_alpha = 3.5 MeV is the alpha particle energy and we assume n_D = n_T = n/2. The energy loss rate is P_loss = 3nT / tau_E (thermal energy content divided by confinement time). Setting P_alpha greater than or equal to P_loss gives the ignition condition. The fusion gain factor Q = P_fusion / P_external quantifies performance: Q = 1 is scientific breakeven, Q = 5 means 20% recirculating power, and Q = infinity is ignition (self-sustaining burn with no external heating needed).

ITER is designed to achieve Q = 10, with n approximately 10^20 m^-3, T approximately 15 keV, and tau_E approximately 3.7 s, giving a triple product of about 5.6 x 10^21 keV s/m^3 -- well above the Lawson threshold.

Interactive Simulations

Solov'ev Equilibrium: Tokamak Flux Surfaces

Python

script.py113 lines

import numpy as np

# Solov'ev equilibrium: analytic solution to the Grad-Shafranov equation
# with p'(psi) = C1 (constant) and FF'(psi) = C2 (constant)
#
# The GS equation becomes:
#   Delta* psi = -mu0 R^2 C1 - C2
#
# A particular solution for up-down symmetric equilibrium:
#   psi(R,Z) = A1/8 * R^4 + A2 * (R^2 * Z^2)
#            + C3 * (R^4 - 4*R^2*Z^2) + C4 * (...)
#
# We use a standard parameterization with R0, a, kappa, delta

R0 = 6.2      # Major radius (m) - ITER-like
a  = 2.0      # Minor radius (m)
kappa = 1.7   # Elongation
delta = 0.33  # Triangularity
B0 = 5.3      # Toroidal field on axis (T)
mu0 = 4 * np.pi * 1e-7

# Solov'ev parameters
# psi = psi0 * [x^4/8 + A*(x^2*y^2/2 - x^4/8 - ln(x)*(1+A*y^2)/2)]
# where x = R/R0, y = Z/R0

# Grid in (R, Z)
NR, NZ = 200, 200
R = np.linspace(R0 - 1.3*a, R0 + 1.3*a, NR)
Z = np.linspace(-1.3*kappa*a, 1.3*kappa*a, NZ)
RR, ZZ = np.meshgrid(R, Z)

# Simplified Solov'ev solution
# psi = C * [(R^2 - R0^2)^2 / (8*R0^2) + Z^2 * R^2 / (2*R0^2) * alpha]
# where alpha controls elongation
alpha_param = kappa**2 / (1 + kappa**2)
C_psi = -mu0 * B0 / (R0 * 3.5)  # Scale factor

x = RR / R0
z = ZZ / R0

# Solov'ev flux function (simplified form)
psi = C_psi * R0**2 * (
    x**4 / 8.0
    + alpha_param * (x**2 * z**2 / 2.0 - x**4 / 8.0)
    + (1 - alpha_param) * (-0.5 * np.log(np.maximum(x, 0.01)))
)

# Normalize so psi=0 on axis, psi=1 on boundary
psi_axis = C_psi * R0**2 * (1.0/8.0 + alpha_param*(0 - 1.0/8.0))
psi_edge = C_psi * R0**2 * (
    ((R0+a)/R0)**4 / 8.0
    + alpha_param * (0 - ((R0+a)/R0)**4 / 8.0)
    + (1 - alpha_param) * (-0.5 * np.log((R0+a)/R0))
)
psi_norm = (psi - psi_axis) / (psi_edge - psi_axis)

# Compute q-profile (approximate for circular cross-section limit)
r_minor = np.linspace(0.05, a, 50)
# q(r) ~ r * B0 / (R0 * B_theta(r))
# For Solov'ev, B_theta ~ mu0 * j0 * r / 2 (uniform current)
j0 = 2 * B0 / (mu0 * R0 * 3.0)  # j such that q(a) ~ 3
B_theta = mu0 * j0 * r_minor / 2
q_profile = r_minor * B0 / (R0 * B_theta)

# Safety factor at key locations
print("Solov'ev Equilibrium for ITER-like Tokamak")
print("=" * 55)
print(f"Major radius R0 = {R0} m")
print(f"Minor radius a  = {a} m")
print(f"Elongation kappa = {kappa}")
print(f"Triangularity delta = {delta}")
print(f"Toroidal field B0 = {B0} T")
print(f"Aspect ratio A = {R0/a:.1f}")
print()
print("Safety factor profile:")
print(f"  q(0)   = {q_profile[0]:.3f}")
print(f"  q(a/2) = {q_profile[len(q_profile)//2]:.3f}")
print(f"  q(a)   = {q_profile[-1]:.3f}")
print()

# Compute magnetic shear
dq_dr = np.gradient(q_profile, r_minor)
shear = r_minor / q_profile * dq_dr
print("Magnetic shear s = (r/q)(dq/dr):")
print(f"  s(a/4) = {shear[len(shear)//4]:.3f}")
print(f"  s(a/2) = {shear[len(shear)//2]:.3f}")
print(f"  s(3a/4) = {shear[3*len(shear)//4]:.3f}")
print()

# Compute plasma beta
# beta = 2*mu0*<p> / B0^2
# For Solov'ev, p(psi) = p0*(1 - psi_norm)
p0 = 3e5  # Pa (central pressure)
p_avg = p0 / 2  # average for linear profile
beta = 2 * mu0 * p_avg / B0**2
beta_N = beta / (1e6 / (a * B0 / (mu0 * 2 * np.pi * R0 * j0 * np.pi * a**2 / 2)))
print(f"Central pressure p0 = {p0/1e3:.0f} kPa")
print(f"Volume-average beta = {beta*100:.2f}%")
print()

# Lawson criterion check
n_e = 1.0e20    # m^-3
T_keV = 15.0    # keV
tau_E = 3.7      # s
triple = n_e * T_keV * tau_E
print("Lawson Criterion Check (ITER parameters):")
print(f"  n_e = {n_e:.1e} m^-3")
print(f"  T = {T_keV} keV")
print(f"  tau_E = {tau_E} s")
print(f"  n*T*tau_E = {triple:.2e} keV s/m^3")
print(f"  Threshold  = 3.0e+21 keV s/m^3")
print(f"  Ignition margin: {triple/3e21:.2f}x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

← Previous Next →