Categories: Physics, Plasma physics.

Grad-Shafranov equation

Nuclear fusion reactors tend to have a torus shape, in which the plasma is confined by a pinch, i.e. by magnetic fields chosen so that the Lorentz force stops particles escaping. Effectively, we are taking a cylindrical screw pinch and bending it into a torus.

We would like to find the equilibrium state of the plasma in the general case of a reactor with toroidal symmetry. Using ideal magnetohydrodynamics (MHD), we start by assuming that the fluid is stationary, and that the confining field B\vb{B} is fixed:

u=0∂u∂t=0∂B∂t=0E=0\begin{aligned} \vb{u} = 0 \qquad \qquad \pdv{\vb{u}}{t} = 0 \qquad \qquad \pdv{\vb{B}}{t} = 0 \qquad \qquad \vb{E} = 0 \end{aligned}

Notice that E=0\vb{E} = 0 is a result of ideal MHD’s generalized Ohm’s law. Under these assumptions, the relevant equations to be solved are Gauss’ law for magnetism, Ampère’s law, and the momentum equation of MHD, respectively:

0=∇⋅Bμ0J=∇×B∇p=J×B\begin{aligned} 0 = \nabla \cdot \vb{B} \qquad \qquad \mu_0 \vb{J} = \nabla \cross \vb{B} \qquad \qquad \nabla p = \vb{J} \cross \vb{B} \end{aligned}

The idea is to analyze them in this order, exploiting toroidal symmetry along the way, to arrive at a general equilibrium condition. Cylindrical polar coordinates (r,θ,z)(r, \theta, z) are a natural choice, with the zz-axis running through the middle of the torus.

As preparation, it is a good idea to write B\vb{B} as the curl of a magnetic vector potential A\vb{A}, which looks like this in cylindrical polar coordinates:

B=∇×A=[1r∂Az∂θ−∂Aθ∂z∂Ar∂z−∂Az∂r1r(∂(rAθ)∂r−∂Ar∂θ)]=[−∂Aθ∂z∂Ar∂z−∂Az∂r1r∂(rAθ)∂r]\begin{aligned} \vb{B} = \nabla \cross \vb{A} = \begin{bmatrix} \displaystyle \frac{1}{r} \pdv{A_z}{\theta} - \pdv{A_\theta}{z} \\ \displaystyle \pdv{A_r}{z} - \pdv{A_z}{r} \\ \displaystyle \frac{1}{r} \Big( \pdv{(r A_\theta)}{r} - \pdv{A_r}{\theta} \Big) \end{bmatrix} = \begin{bmatrix} \displaystyle - \pdv{A_\theta}{z} \\ \displaystyle \pdv{A_r}{z} - \pdv{A_z}{r} \\ \displaystyle \frac{1}{r} \pdv{(r A_\theta)}{r} \end{bmatrix} \end{aligned}

Here, it is convenient to define the so-called stream function ψ\psi as follows:

ψ≡rAθ\begin{aligned} \boxed{ \psi \equiv r A_\theta } \end{aligned}

Such that B\vb{B} can be written as below, where we will regard BθB_\theta as a given quantity:

B=[−1r∂ψ∂zBθ1r∂ψ∂r]whereBθ=∂Ar∂z−∂Az∂r\begin{aligned} \vb{B} = \begin{bmatrix} \displaystyle -\frac{1}{r} \pdv{\psi}{z} \\ B_\theta \\ \displaystyle \frac{1}{r} \pdv{\psi}{r} \end{bmatrix} \qquad \mathrm{where} \qquad B_\theta = \pdv{A_r}{z} - \pdv{A_z}{r} \end{aligned}

Inserting this into Gauss’ law, we see that it is trivially satisfied, thanks to circular symmetry guaranteeing that ∂Bθ/∂θ=0\ipdv{B_\theta}{\theta} = 0:

0=∇⋅B=−1r∂∂r(rr∂ψ∂z)+1r∂Bθ∂θ+∂∂z(1r∂ψ∂r)=−1r∂2ψ∂r∂z+1r∂2ψ∂z∂r=0\begin{aligned} 0 = \nabla \cdot \vb{B} &= - \frac{1}{r} \pdv{}{r}\bigg( \frac{r}{r} \pdv{\psi}{z} \bigg) + \frac{1}{r} \pdv{B_\theta}{\theta} + \pdv{}{z}\bigg( \frac{1}{r} \pdv{\psi}{r} \bigg) \\ &= - \frac{1}{r} \mpdv{\psi}{r}{z} + \frac{1}{r} \mpdv{\psi}{z}{r} = 0 \end{aligned}

What matters is that we have expressions for the components of B\vb{B}. Moving on, to find the current density J\vb{J}, we use Ampère’s law and symmetry to get:

J=1μ0∇×B=1μ0[1r∂Bz∂θ−∂Bθ∂z∂Br∂z−∂Bz∂r1r(∂(rBθ)∂r−∂Br∂θ)]=1μ0[0∂Br∂z−∂Bz∂r1r∂(rBθ)∂r]\begin{aligned} \vb{J} = \frac{1}{\mu_0} \nabla \cross \vb{B} = \frac{1}{\mu_0} \begin{bmatrix} \displaystyle \frac{1}{r} \pdv{B_z}{\theta} - \pdv{B_\theta}{z} \\ \displaystyle \pdv{B_r}{z} - \pdv{B_z}{r} \\ \displaystyle \frac{1}{r} \Big( \pdv{(r B_\theta)}{r} - \pdv{B_r}{\theta} \Big) \end{bmatrix} = \frac{1}{\mu_0} \begin{bmatrix} \displaystyle 0 \\ \displaystyle \pdv{B_r}{z} - \pdv{B_z}{r} \\ \displaystyle \frac{1}{r} \pdv{(r B_\theta)}{r} \end{bmatrix} \end{aligned}

Where we have assumed that BθB_\theta depends only on rr, not zz or θ\theta. Substituting this into the MHD momentum equation gives the following pressure gradient ∇p\nabla p:

∇p=J×B=[JθBz−JzBθJzBr−JrBzJrBθ−JθBr]=[JθBz−JzBθJzBr−JθBr]\begin{aligned} \nabla p &= \vb{J} \cross \vb{B} = \begin{bmatrix} J_\theta B_z - J_z B_\theta \\ J_z B_r - J_r B_z \\ J_r B_\theta - J_\theta B_r \end{bmatrix} = \begin{bmatrix} J_\theta B_z - J_z B_\theta \\ J_z B_r \\ - J_\theta B_r \end{bmatrix} \end{aligned}

Now, the idea is to focus on this rr-component to get an equation for ψ\psi, whose solution can then be used to calculate the θ\theta and zz-components of ∇p\nabla p. Therefore, we evaluate:

∂p∂r=JθBz−JzBθ=1μ0(∂Br∂z−∂Bz∂r)Bz−1μ0r∂(rBθ)∂rBθ=−1μ0(∂∂z(1r∂ψ∂z)+∂∂r(1r∂ψ∂r))1r∂ψ∂r−1μ0r∂(rBθ)∂rBθ=−1μ0r(1r∂2ψ∂z2+∂∂r(1r∂ψ∂r))∂ψ∂r−1μ0r∂(rBθ)∂rBθ\begin{aligned} \pdv{p}{r} &= J_\theta B_z - J_z B_\theta \\ &= \frac{1}{\mu_0} \bigg( \pdv{B_r}{z} - \pdv{B_z}{r} \bigg) B_z - \frac{1}{\mu_0 r} \pdv{(r B_\theta)}{r} B_\theta \\ &= - \frac{1}{\mu_0} \bigg( \pdv{}{z}\Big(\frac{1}{r} \pdv{\psi}{z}\Big) + \pdv{}{r}\Big(\frac{1}{r} \pdv{\psi}{r}\Big) \bigg) \frac{1}{r} \pdv{\psi}{r} - \frac{1}{\mu_0 r} \pdv{(r B_\theta)}{r} B_\theta \\ &= - \frac{1}{\mu_0 r} \bigg( \frac{1}{r} \pdvn{2}{\psi}{z} + \pdv{}{r}\Big( \frac{1}{r} \pdv{\psi}{r} \Big) \bigg) \pdv{\psi}{r} - \frac{1}{\mu_0 r} \pdv{(r B_\theta)}{r} B_\theta \end{aligned}

By using the chain rule to rewrite ∂/∂r=(∂ψ/∂r)  ∂/∂ψ\ipdv{}{r}= (\ipdv{\psi}{r}) \; \ipdv{}{\psi}, we get ∂ψ/∂r\ipdv{\psi}{r} in each term:

∂ψ∂r∂p∂ψ=−1μ0r(1r∂2ψ∂z2+∂∂r(1r∂ψ∂r))∂ψ∂r−1μ0r∂ψ∂r∂(rBθ)∂ψBθ\begin{aligned} \pdv{\psi}{r} \pdv{p}{\psi} &= - \frac{1}{\mu_0 r} \bigg( \frac{1}{r} \pdvn{2}{\psi}{z} + \pdv{}{r}\Big( \frac{1}{r} \pdv{\psi}{r} \Big) \bigg) \pdv{\psi}{r} - \frac{1}{\mu_0 r} \pdv{\psi}{r} \pdv{(r B_\theta)}{\psi} B_\theta \end{aligned}

Dividing out ∂ψ/∂r\ipdv{\psi}{r} and multiplying by μ0r2\mu_0 r^2 leads us to the Grad-Shafranov equation, which gives the equilibrium condition of a plasma in a toroidal reactor:

∂2ψ∂z2+r∂∂r(1r∂ψ∂r)=−μ0r2∂p∂ψ−r∂(rBθ)∂ψBθ\begin{aligned} \boxed{ \pdvn{2}{\psi}{z} + r \pdv{}{r}\bigg( \frac{1}{r} \pdv{\psi}{r} \bigg) = - \mu_0 r^2 \pdv{p}{\psi} - r \pdv{(r B_\theta)}{\psi} B_\theta } \end{aligned}

Weirdly, ψ\psi appears both as an unknown and as a differentiation variable, but this equation can still be solved analytically by assuming a certain ψ\psi-dependence of pp and rBθr B_\theta.

Suppose that BθB_\theta is induced by a poloidal electrical current IpolI_\mathrm{pol}, i.e. a current around the “tube” of the torus, then, assuming IpolI_\mathrm{pol} only depends on rr, we have:

Bθ=μ0Ipol(r)2πr\begin{aligned} B_\theta = \frac{\mu_0 I_\mathrm{pol}(r)}{2 \pi r} \end{aligned}

Inserting this into the Grad-Shafranov equation yields its following alternative form:

∂2ψ∂z2+r∂∂r(1r∂ψ∂r)=−μ0r2∂p∂ψ−μ028π2∂Ipol2∂ψ\begin{aligned} \boxed{ \pdvn{2}{\psi}{z} + r \pdv{}{r}\bigg( \frac{1}{r} \pdv{\psi}{r} \bigg) = - \mu_0 r^2 \pdv{p}{\psi} - \frac{\mu_0^2}{8 \pi^2} \pdv{I_\mathrm{pol}^2}{\psi} } \end{aligned}

References

  1. M. Salewski, A.H. Nielsen, Plasma physics: lecture notes, 2021, unpublished.