---
title: "Grad-Shafranov equation"
sort_title: "Grad-Shafranov equation"
date: 2022-03-06
categories:
- Physics
- Plasma physics
layout: "concept"
---

Nuclear fusion reactors tend to have a torus shape,
in which the plasma is confined by a **pinch**,
i.e. by [magnetic fields](/know/concept/magnetic-field/)
chosen so that the [Lorentz force](/know/concept/lorentz-force/)
stops particles escaping.
Effectively, we are taking a cylindrical [screw pinch](/know/concept/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](/know/concept/magnetohydrodynamics/) (MHD),
we start by assuming that the fluid is stationary,
and that the confining field $$\vb{B}$$ is fixed:

$$\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 $$\vb{E} = 0$$ is a result of the ideal generalized Ohm's law.
Under these assumptions, the relevant MHD equations to be solved are
Gauss' law for magnetism, Ampère's law, and the MHD momentum equation, respectively:

$$\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 goal is to analyze them in this order,
exploiting toroidal symmetry along the way,
to arrive at a general equilibrium condition.
[Cylindrical polar coordinates](/know/concept/cylindrical-polar-coordinates/) $$(r, \theta, z)$$
are a natural choice, with the $$z$$-axis running through the middle of the torus.

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

$$\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:

$$\begin{aligned}
    \boxed{
        \psi
        \equiv r A_\theta
    }
\end{aligned}$$

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

$$\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 $$\ipdv{B_\theta}{\theta} = 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 $$\vb{B}$$.
Moving on, to find the current density $$\vb{J}$$,
we use Ampère's law and symmetry to get:


$$\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_\theta$$ depends only on $$r$$, not $$z$$ or $$\theta$$.
Substituting this into the MHD momentum equation
gives the following pressure gradient $$\nabla p$$:

$$\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 $$r$$-component to get an equation for $$\psi$$,
whose solution can then be used to calculate the $$\theta$$ and $$z$$-components of $$\nabla p$$.
Therefore, we evaluate:

$$\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 $$\ipdv{}{r}= (\ipdv{\psi}{r}) \; \ipdv{}{\psi}$$,
we get $$\ipdv{\psi}{r}$$ in each term:

$$\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 $$\ipdv{\psi}{r}$$ and multiplying by $$\mu_0 r^2$$
leads us to the **Grad-Shafranov equation**,
which gives the equilibrium condition of a plasma in a toroidal reactor:

$$\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 $$p$$ and $$r B_\theta$$.

Suppose that $$B_\theta$$ is induced by a poloidal electrical current $$I_\mathrm{pol}$$,
i.e. a current around the "tube" of the torus,
then, assuming $$I_\mathrm{pol}$$ only depends on $$r$$, we have:

$$\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:

$$\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.