Part 3, Chapter 3

MHD Equilibrium

Force balance and magnetic confinement configurations

3.1 Static Equilibrium Condition

In MHD equilibrium, all time derivatives vanish and flow velocity u = 0. The momentum equation reduces to:

$$\mathbf{J} \times \mathbf{B} = \nabla p$$

Force Balance: Lorentz force balances pressure gradient

Consequence 1

$$\mathbf{B} \cdot \nabla p = 0$$

Pressure is constant along field lines

Consequence 2

$$\mathbf{J} \cdot \nabla p = 0$$

Pressure is constant along current lines

3.2 Magnetic Surfaces

Since both B and J are perpendicular to ∇p, they lie on surfaces of constant pressure called magnetic surfaces or flux surfaces.

Properties of Magnetic Surfaces

Nested toroidal surfaces in toroidal geometry
Field lines wrap helically around torus
Characterized by rotational transform ι or safety factor q
Rational surfaces where field lines close on themselves

Safety Factor

$$q = \frac{\text{toroidal turns}}{\text{poloidal turns}} = \frac{r B_\phi}{R B_\theta}$$

Measures pitch of field line helix; q > 1 typically needed for stability

3.3 Grad-Shafranov Equation

For axisymmetric toroidal equilibria, the force balance reduces to a single scalar equation:

$$R^2 \nabla \cdot \left(\frac{\nabla \psi}{R^2}\right) = -\mu_0 R^2 \frac{dp}{d\psi} - F\frac{dF}{d\psi}$$

Grad-Shafranov Equation

Poloidal Flux

$$\psi = \text{flux function}$$

Labels magnetic surfaces

Poloidal Current

$$F(\psi) = RB_\phi$$

Toroidal field function

3.4 Plasma Beta

The ratio of plasma pressure to magnetic pressure is a key dimensionless parameter:

$$\beta = \frac{p}{B^2/(2\mu_0)} = \frac{\text{thermal pressure}}{\text{magnetic pressure}}$$

Low β

β ≪ 1: Magnetic pressure dominates (tokamaks typically β ~ few %)

High β

β ~ 1: Comparable pressures (spheromaks, RFPs)

β Limit

Maximum β before instabilities (Troyon limit)

3.5 Simple Equilibria

θ-Pinch

$$p + \frac{B_z^2}{2\mu_0} = \frac{B_0^2}{2\mu_0}$$

Axial field confines plasma radially; no end confinement

Z-Pinch

$$\frac{d}{dr}\left(p + \frac{B_\theta^2}{2\mu_0}\right) = -\frac{B_\theta^2}{\mu_0 r}$$

Azimuthal field from axial current; inherently unstable

Screw Pinch

$$\frac{d}{dr}\left(p + \frac{B_\theta^2 + B_z^2}{2\mu_0}\right) = -\frac{B_\theta^2}{\mu_0 r}$$

Combined θ and Z-pinch; basic tokamak geometry

3.6 Bennett Relation

For a uniform pressure Z-pinch in equilibrium:

$$I^2 = \frac{8\pi N k_B (T_e + T_i)}{\mu_0}$$

Bennett Pinch Condition

Where N is the line density (particles per unit length) and I is the total current. This gives the minimum current needed to confine the plasma.

Key Takeaways

✓ Static equilibrium: J × B = ∇p
✓ B and J lie on magnetic flux surfaces of constant p
✓ Grad-Shafranov equation describes axisymmetric equilibria
✓ β = p/(B²/2μ₀) measures confinement efficiency
✓ Pinch configurations: θ-pinch, Z-pinch, screw pinch

Interactive Simulations

Grad-Shafranov: Tokamak Equilibrium Profiles

Python

script.py60 lines

import numpy as np

# Simplified tokamak equilibrium profiles
# Safety factor q(r) = r*B_phi / (R*B_theta)
# Pressure balance: dp/dr = -j_z * B_theta + ...

mu0 = 4 * np.pi * 1e-7

# Tokamak parameters (ITER-like)
R0 = 6.2        # major radius (m)
a = 2.0          # minor radius (m)
B0 = 5.3         # toroidal field on axis (T)
Ip = 15e6        # plasma current (A)
beta_0 = 0.05    # central beta
n0 = 1e20        # central density (m^-3)

print("Tokamak Equilibrium Profiles (ITER-like)")
print("=" * 58)
print(f"R0 = {R0} m, a = {a} m, B0 = {B0} T")
print(f"Ip = {Ip/1e6:.0f} MA, beta_0 = {beta_0}")
print()

N = 20
r = np.linspace(0.01, a, N)
rho = r / a  # normalized radius

# Profiles (parabolic model)
# q(r) = q0 * (1 + (r/a)^2 * (q_a/q0 - 1))
q0 = 1.0   # central safety factor
q_a = 3.5  # edge safety factor

q = q0 * (1 + rho**2 * (q_a/q0 - 1))

# Pressure profile: p(r) = p0 * (1 - (r/a)^2)^2
P_mag = B0**2 / (2 * mu0)
p0 = beta_0 * P_mag
p = p0 * (1 - rho**2)**2

# Poloidal field from q profile
B_theta = r * B0 / (R0 * q)

# Current density (from Ampere)
# j_phi ~ (1/mu0*r) * d(r*B_theta)/dr
j_phi = np.zeros_like(r)
for i in range(1, len(r)-1):
    dr = r[i+1] - r[i-1]
    j_phi[i] = (1/(mu0*r[i])) * (r[i+1]*B_theta[i+1] - r[i-1]*B_theta[i-1]) / dr

print(f"{'r/a':>6} {'q':>8} {'P(kPa)':>10} {'B_pol(T)':>10} {'j(MA/m2)':>10} {'beta_pol':>10}")
print("-" * 56)

for i in range(N):
    beta_pol = 2 * mu0 * p[i] / B_theta[i]**2 if B_theta[i] > 0 else 0
    print(f"{rho[i]:6.2f} {q[i]:8.2f} {p[i]/1000:10.1f} {B_theta[i]:10.3f} {j_phi[i]/1e6:10.3f} {beta_pol:10.2f}")

print(f"\nShafranov shift ~ a * beta_pol * a/(2*R0)")
beta_pol_avg = 2 * mu0 * np.mean(p) / (B_theta[-1]**2)
Delta = a * beta_pol_avg * a / (2 * R0)
print(f"Estimated Shafranov shift: {Delta*100:.1f} cm")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Screw Pinch Equilibrium and Suydam Criterion

Fortran

program.f9074 lines

program screw_pinch
  implicit none
  ! Screw pinch: cylindrical equilibrium with B_z and B_theta
  ! Radial force balance: dP/dr + d(B^2/(2*mu0))/dr + B_theta^2/(mu0*r) = 0
  ! Suydam stability: r*q'^2/q^2 * B_z^2/(mu0) > -8*dP/dr
  real(8) :: mu0, pi
  real(8) :: r, dr, a, B0, B_theta_a
  real(8) :: P, B_z, B_theta, q_val, j_z, j_theta
  real(8) :: P_prev, B_z_prev, B_theta_prev
  real(8) :: dPdr, dq, suydam_lhs, suydam_rhs
  integer :: i, N

mu0 = 4.0d0 * 3.14159265d0 * 1.0d-7
  pi = 4.0d0 * atan(1.0d0)

a = 0.1d0        ! pinch radius (m)
  B0 = 1.0d0       ! axial field on axis (T)
  B_theta_a = 0.5d0  ! poloidal field at edge (T)
  N = 20

dr = a / dble(N)

write(*,'(A)') "Screw Pinch Equilibrium"
  write(*,'(A)') "========================"
  write(*,'(A,F6.2,A)') "Radius a = ", a*100, " cm"
  write(*,'(A,F6.2,A)') "B_z(0) = ", B0, " T"
  write(*,'(A,F6.2,A)') "B_theta(a) = ", B_theta_a, " T"
  write(*,*)

write(*,'(A6,A10,A10,A10,A12,A10)') &
    "r/a", "B_z(T)", "B_th(T)", "P(kPa)", "q", "Suydam"
  write(*,'(A)') repeat("-", 60)

do i = 1, N
    r = dble(i) * dr

! Model profiles
    B_z = B0 * (1.0d0 - 0.3d0 * (r/a)**2)
    B_theta = B_theta_a * (r/a)

! Pressure from force balance (simplified)
    P = (B0**2 - B_z**2) / (2.0d0*mu0) + &
        B_theta_a**2 * (1.0d0 - (r/a)**2) / (4.0d0*mu0)
    P = max(P, 0.0d0)

! Safety factor
    q_val = r * B_z / (a * B_theta) * (a/r)
    if (r > 0.01d0*a) then
      q_val = r * B_z / (r * B_theta) ! simplified
    else
      q_val = B0 * a / (B_theta_a * r)
    end if

! Suydam criterion (simplified check)
    dPdr = -B0**2 * 0.6d0 * r / (a**2 * mu0) &
           - B_theta_a**2 * r / (2.0d0*mu0*a**2)

if (r > 0.05d0*a) then
      suydam_lhs = B_z**2 / (mu0 * q_val**2)
      suydam_rhs = -8.0d0 * dPdr * r
      if (suydam_lhs > suydam_rhs) then
        write(*,'(F6.2,F10.4,F10.4,F10.2,F12.3,A10)') &
          r/a, B_z, B_theta, P/1000, q_val, " Stable"
      else
        write(*,'(F6.2,F10.4,F10.4,F10.2,F12.3,A10)') &
          r/a, B_z, B_theta, P/1000, q_val, " Unstable"
      end if
    else
      write(*,'(F6.2,F10.4,F10.4,F10.2,F12.3,A10)') &
        r/a, B_z, B_theta, P/1000, q_val, " --"
    end if
  end do
end program screw_pinch

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← MHD Equations Next: MHD Waves →