Grad–Shafranov equation

The Grad–Shafranov equation (H. Grad and H. Rubin (1958); Vitalii Dmitrievich Shafranov (1966)) is the equilibrium equation in ideal magnetohydrodynamics (MHD) for a two dimensional plasma, for example the axisymmetric toroidal plasma in a tokamak. This equation takes the same form as the Hicks equation from fluid dynamics. This equation is a two-dimensional, nonlinear, elliptic partial differential equation obtained from the reduction of the ideal MHD equations to two dimensions, often for the case of toroidal axisymmetry (the case relevant in a tokamak). Taking

$(r,\theta ,z)$

as the cylindrical coordinates, the flux function

$\psi$

is governed by the equation,

$\frac{{\partial }^2\psi }{\partial r^2}-\frac{1}{r}\frac{\partial \psi }{\partial r}+\frac{{\partial }^2\psi }{\partial z^2}=-{\mu }_0{r}^2\frac{dp}{d\psi }-\frac{1}{2}\frac{d{F}^2}{d\psi },$

where

${\mu }_0$

is the magnetic permeability,

$p(\psi )$

is the pressure,

$F(\psi )=r{B}_{\theta }$

and the magnetic field and current are, respectively, given by

$$\begin{array}{rl}\mathbf{B}& =\frac{1}{r}\mathrm{\nabla }\psi \times {\hat{\mathbf{e}}}_{\theta }+\frac{F}{r}{\hat{\mathbf{e}}}_{\theta },\\ {\mu }_0\mathbf{J}& =\frac{1}{r}\frac{dF}{d\psi }\mathrm{\nabla }\psi \times {\hat{\mathbf{e}}}_{\theta }-[\frac{\partial }{\partial r}\left(\frac{1}{r}\frac{\partial \psi }{\partial r}\right)+\frac{1}{r}\frac{{\partial }^2\psi }{\partial z^2}]{\hat{\mathbf{e}}}_{\theta }.\end{array}$$

The nature of the equilibrium, whether it be a tokamak, reversed field pinch, etc. is largely determined by the choices of the two functions $F(\psi )$ and $p(\psi )$ as well as the boundary conditions.

Derivation (in Cartesian coordinates)

In the following, it is assumed that the system is 2-dimensional with $z$ as the invariant axis, i.e. $\frac{\partial }{\partial z}$ produces 0 for any quantity. Then the magnetic field can be written in cartesian coordinates as $$\mathbf{B}=(\frac{\partial A}{\partial y},-\frac{\partial A}{\partial x},B_z(x,y)),$$ or more compactly, $$\mathbf{B}=\mathrm{\nabla }A\times \hat{\mathbf{z}}+B_z\hat{\mathbf{z}},$$ where $A(x,y)\hat{\mathbf{z}}$ is the vector potential for the in-plane (x and y components) magnetic field. Note that based on this form for B we can see that A is constant along any given magnetic field line, since $\mathrm{\nabla }A$ is everywhere perpendicular to B. (Also note that -A is the flux function $\psi$ mentioned above.) Two dimensional, stationary, magnetic structures are described by the balance of pressure forces and magnetic forces, i.e.: $$\mathrm{\nabla }p=\mathbf{j}\times \mathbf{B},$$ where p is the plasma pressure and j is the electric current. It is known that p is a constant along any field line, (again since $\mathrm{\nabla }p$ is everywhere perpendicular to B). Additionally, the two-dimensional assumption ($\frac{\partial }{\partial z}=0$) means that the z- component of the left hand side must be zero, so the z-component of the magnetic force on the right hand side must also be zero. This means that ${\mathbf{j}}_{\perp }\times {\mathbf{B}}_{\perp }=0$, i.e. ${\mathbf{j}}_{\perp }$ is parallel to ${\mathbf{B}}_{\perp }$. The right hand side of the previous equation can be considered in two parts: $$\mathbf{j}\times \mathbf{B}=j_z(\hat{\mathbf{z}}\times {\mathbf{B}}_{\mathbf{\perp }})+{\mathbf{j}}_{\mathbf{\perp }}\times \hat{\mathbf{z}}{B}_z,$$ where the $\perp$ subscript denotes the component in the plane perpendicular to the $z$-axis. The $z$ component of the current in the above equation can be written in terms of the one-dimensional vector potential as $$j_z=-\frac{1}{{\mu }_0}{\mathrm{\nabla }}^{2}A.$$ The in plane field is $${\mathbf{B}}_{\perp }=\mathrm{\nabla }A\times \hat{\mathbf{z}},$$ and using Maxwell–Ampère's equation, the in plane current is given by $${\mathbf{j}}_{\perp }=\frac{1}{{\mu }_0}\mathrm{\nabla }{B}_z\times \hat{\mathbf{z}}.$$ In order for this vector to be parallel to ${\mathbf{B}}_{\perp }$ as required, the vector $\mathrm{\nabla }{B}_z$ must be perpendicular to ${\mathbf{B}}_{\perp }$, and $B_z$ must therefore, like $p$, be a field-line invariant. Rearranging the cross products above leads to $$\hat{\mathbf{z}}\times {\mathbf{B}}_{\perp }=\mathrm{\nabla }A-(\hat{\mathbf{z}}\cdot \mathrm{\nabla }A)\hat{\mathbf{z}}=\mathrm{\nabla }A,$$ and $${\mathbf{j}}_{\perp }\times B_z\hat{\mathbf{z}}=\frac{B_z}{{\mu }_0}(\hat{\mathbf{z}}\cdot \mathrm{\nabla }{B}_z)\hat{\mathbf{z}}-\frac{1}{{\mu }_0}{B}_z\mathrm{\nabla }{B}_z=-\frac{1}{{\mu }_0}{B}_z\mathrm{\nabla }{B}_z.$$ These results can be substituted into the expression for $\mathrm{\nabla }p$ to yield: $$\mathrm{\nabla }p=-\left[\frac{1}{{\mu }_0}{\mathrm{\nabla }}^{2}A\right]\mathrm{\nabla }A-\frac{1}{{\mu }_0}{B}_z\mathrm{\nabla }{B}_z.$$ Since $p$ and $B_z$ are constants along a field line, and functions only of $A$, hence $\mathrm{\nabla }p=\frac{dp}{dA}\mathrm{\nabla }A$ and $\mathrm{\nabla }{B}_z=\frac{d{B}_z}{dA}\mathrm{\nabla }A$. Thus, factoring out $\mathrm{\nabla }A$ and rearranging terms yields the Grad–Shafranov equation: $${\mathrm{\nabla }}^{2}A=-{\mu }_0\frac{d}{dA}(p+\frac{B_z^2}{2{\mu }_0}).$$

Derivation in contravariant representation

This derivation is only used for Tokamaks, but it can be enlightening. Using the definition of 'The Theory of Toroidally Confined Plasmas 1:3'(Roscoe White), Writing $\overrightarrow{B}$ by contravariant basis $(\mathrm{\nabla }\mathrm{\Psi },\mathrm{\nabla }\varphi ,\mathrm{\nabla }\zeta )$: $$\overrightarrow{B}=\mathrm{\nabla }\mathrm{\Psi }\times \mathrm{\nabla }\varphi +\overline{F}\mathrm{\nabla }\varphi ,$$ we have $\overrightarrow{j}$: $$\begin{gathered} {\mu }_0\overrightarrow{j}=\mathrm{\nabla }\times \overrightarrow{B}=-{\mathrm{\Delta }}^{*}\mathrm{\Psi }\mathrm{\nabla }\varphi +\mathrm{\nabla }\overline{F}\times \mathrm{\nabla }\varphi \text{, where} \\ {\mathrm{\Delta }}^{*}=r{\partial }_r(r^{-1}{\partial }_r)+\underset{\varphi }{\overset{2}{\partial }}\text{;} \end{gathered}$$ then force balance equation: $${\mu }_0\overrightarrow{j}\times \overrightarrow{B}={\mu }_0\mathrm{\nabla }p\text{.}$$ Working out, we have: $$-{\mathrm{\Delta }}^{*}\mathrm{\Psi }=\overline{F}\frac{d\overline{F}}{d\mathrm{\Psi }}+{\mu }_0{R}^2\frac{dp}{d\mathrm{\Psi }}\text{.}$$

Further reading

• Grad, H., and Rubin, H. (1958) Hydromagnetic Equilibria and Force-Free Fields Archived 2023-06-21 at the Wayback Machine. Proceedings of the 2nd UN Conf. on the Peaceful Uses of Atomic Energy, Vol. 31, Geneva: IAEA p. 190.
• Shafranov, V.D. (1966) Plasma equilibrium in a magnetic field, Reviews of Plasma Physics, Vol. 2, New York: Consultants Bureau, p. 103.
• Woods, Leslie C. (2004) Physics of plasmas, Weinheim: WILEY-VCH Verlag GmbH & Co. KGaA, chapter 2.5.4
• Haverkort, J.W. (2009) Axisymmetric Ideal MHD Tokamak Equilibria. Notes about the Grad–Shafranov equation, selected aspects of the equation and its analytical solutions.
• Haverkort, J.W. (2009) Axisymmetric Ideal MHD equilibria with Toroidal Flow. Incorporation of toroidal flow, relation to kinetic and two-fluid models, and discussion of specific analytical solutions.