Two-dimensional shaping of Solov’ev equilibrium with vacuum using external coils

Jiaxing Liu International Joint Research Laboratory of Magnetic Confinement Fusion and Plasma Physics, State Key Laboratory of Advanced Electromagnetic Engineering and Technology, School of Electrical and Electronic Engineering, Huazhong University of Science and Technology, Wuhan, 430074, China Ping Zhu zhup@hust.edu.cn International Joint Research Laboratory of Magnetic Confinement Fusion and Plasma Physics, State Key Laboratory of Advanced Electromagnetic Engineering and Technology, School of Electrical and Electronic Engineering, Huazhong University of Science and Technology, Wuhan, 430074, China Department of Engineering Physics, University of Wisconsin-Madison, Madison, Wisconsin 53706, United States of America Haolong Li College of Physics and Optoelectronic Engineering, Shenzhen University, Shenzhen 518060, China

Abstract

In this work, we demonstrate a method for constructing the Solov’ev equilibrium with any given 2D shape surrounded by a vacuum region using external poloidal field coils. The computational domain consists of two parts: the plasma region, where the solution is the same as the Solov’ev solution, and the vacuum region, where the magnetic field generated by external coils as well as plasma current is determined using the Green function method through a matching condition near the separatrix. However, the method is not limited to the Solov’ev equilibrium in particular. The accuracy, efficiency, and robustness of such a scheme suggest that this method may be applied to the 2D shaping of tokamak plasma with vacuum region using external coils in general.

For a tokamak plasma, its stability and transport can be significantly affected by the two-dimensional (2D) shaping of equilibrium poloidal cross-section[1, 2, 3, 4]. Thus it is often desirable to construct the tokamak equilibrium with any prescribed poloidal shape for its plasma-vacuum interface using external coils in both theory and experiments.

2D shaping of tokamak equilibrium has been a subject of extensive studies. For certain pressure and current profiles, several analytical solutions of the Grad-Shafranov (G-S) equation allowing arbitrary plasma shape have been derived, where the shaping parameters, including elongation and triangularity, are determined using the expansion coefficients of the homogeneous solutions [5, 6, 7, 8]. In simulation and experiments, 2D shaping studies have been often carried out in the design and optimization of equilibrium configuration feedback control schemes, which have been implemented in most tokamak devices (e.g. EAST[9] and DIII-D[10]). one common such scheme usually involves the measurement and calculation of poloidal flux at specific control points[11, 12, 13]. For the given plasma current and pressure profiles, any deviation of the poloidal shaping from its target at a certain control surface can be corrected by adjusting the external coil current, and the exact adjustment can be determined by the matching condition at the control points. Such a method has been found successful in practice [10, 9]. However, the direct relation between the 2D plasma shape and external coil configurations, which would allow the identification of their overall geometric and physical connections, has yet to be explicitly established.

This paper presents a method for finding the direct and explicit relation between the tokamak shaping and the external poloidal field coil configuration, taking the well-known Solov’ev equilibrium as an example. As one of the analytical solutions of the G-S equation, Solov’ev equilibrium is often used in the benchmark of the numerical G-S equation solver and the analysis of tokamak equilibrium, such as the study on 2D plasma shaping[6, 8]. The conventional Solov’ev solution, however, is derived with toroidal current filling all space, which is in contrast to the tokamak experiment where there is virtually a vacuum out of separatrix. Even though many studies have been done on the relationship between Solov’ev solutions and plasma shape[6, 8], most of these studies are based on the conventional Solov’ev solution without vacuum region. Xu and Fitzpatrick find an overall equilibrium solution with the Solov’ev solution in plasma along with a vacuum region[14] using the Green function and multiple expansion method, which may fail when the external coils are close to the plasma boundary as in realistic experiments. Instead, we directly use the Green function method for both plasma current and external poloidal field coils, which allows us to not only find the matching vacuum solution to a given Solov’ev equilibrium inside the separatrix but also establish the direct and explicit relation between the plasma shape and the external coil configuration. In general, this method can be directly applied to obtain tokamak equilibrium with desired 2D plasma shape in both simulation and experiment.

The rest of this paper is organized as follows: first, the numerical method is presented for finding the relation between plasma shaping and external poloidal field coil configuration. Then, we demonstrate the usage of the method and evaluate the associated numerical error using an example case of Solov’ev equilibrium. Finally, we apply the method to discuss the explicit relations between several typical plasma shaping and their corresponding external poloidal field coil configurations.

Numerical method: For the MHD equilibrium with toroidal axisymmetry, the poloidal magnetic flux function $\psi$ is governed by the G-S equation

R\frac{\partial}{\partial R}\left(\frac{1}{R}\frac{\partial\psi}{\partial R}\right)+\frac{\partial^{2}\psi}{\partial Z^{2}}=-RJ_{\phi}

(1)

with

RJ_{\phi}=R^{2}\frac{dp}{d\psi}-F\frac{dF}{d\psi}

(2)

Here, $\left(R,\phi,Z\right)$ corresponds to the right-hand cylindrical coordinate system. The above equation is written in the dimensionless form, where the magnetic field, poloidal magnetic flux, plasma pressure $p$ , toroidal current density $J_{\phi}$ , and poloidal current are normalized by the constants $R_{0},B_{0},R_{0}^{2}B_{0},B_{0}^{2}/\mu_{0},B_{0}/\left(\mu_{0}R_{0}\right)$ , and $R_{0}B_{0}$ respectively, with $R_{0},B_{0}$ being the major radius and magnetic field strength at the magnetic axis.

The solutions to the G-S equation consist of the special and homogeneous contributions:

\psi=\psi_{0}+\sum_{i=1}^{n}c_{i}\psi_{i}\qquad\left(i=1,2,3...,n\right)

(3)

where $\psi_{0}$ is the special solution specified by the pressure and the current profiles, and $\psi_{i}$ is a set of homogeneous solutions of the G-S equation (1) that satisfy certain boundary conditions. The homogeneous solution used here assumes that the overall solution can be written as a Taylor series, starting from a constant and increasing up to and including $n$ -th order terms[6, 15, 16, 17, 18]. Next, we can solve the G-S equation by finding the corresponding matching solutions with vacuum in the form of Eq. (3) from two sub-regions of the computational domain. The first part is the plasma region with closed magnetic field lines, and the second is the vacuum region with open magnetic field lines.

$\psi$ in the closed field line region: In this region, we first determined the special solution given by the pressure and current profiles. In the case of Solov’ev equilibrium, the pressure and current profiles are as follows:

\begin{array}[]{rl}-\frac{dp}{d\psi}&=a\\ F\frac{dF}{d\psi}&=b\end{array}

(4)

For any given special solution, the 2D plasma shaping is further determined by the homogeneous solution or in terms of its expansion in Eq. 3, a particular set of the coefficients $c_{i}$ . On the other hand, all the desired local and global geometric features of the plasma shaping can be used to construct the corresponding set of constraints on the flux function or its partial derivatives at locations, such as the X points, the magnetic axis, or the surface, which allows the specification of the homogeneous solution to the extent of an equal number of coefficients $c_{i}$ . This completes the equilibrium solution in the plasma region on the closed field lines inside the separatrix.

$\psi$ in the open field line region: The magnetic flux $\psi$ in the open field line region is determined by the plasma current and the external poloidal field coils, which, in general, can be calculated using the Green function method:

\psi=\psi_{p}+\psi_{c}

(5)

\psi_{p}(R,Z)=\iint_{\Omega}G\left(R,Z;R^{\prime},z^{\prime}\right)J_{\phi}\left(R^{\prime},Z^{\prime}\right)dR^{\prime}dZ^{\prime}

(6)

\psi_{c}(R,Z)=\sum_{i=1}^{N_{c}}G\left(R,Z;R_{i}^{c},Z_{i}^{c}\right)I_{i}

(7)

where $\Omega$ represents the plasma region inside the separatrix, $J_{\phi}$ is the toroidal current density, $N_{c}$ is the total number of the external coils, and $I_{i}$ is the current in each of the external coils, the Green function $G\left(R,Z;R^{\prime},z^{\prime}\right)$ is given by [19]

G\left(\mathbf{R};\mathbf{R}^{\prime}\right)=\frac{1}{2\pi}\frac{\sqrt{RR^{\prime}}}{k}\left[\left(2-k^{2}\right)K(k)-2E(k)\right]

(8)

where $K\left(k\right)$ and $E\left(k\right)$ are the complete elliptic integrals of the first and second kind, and the argument is defined by

k^{2}=\frac{4RR^{\prime}}{\left[\left(R+R^{\prime}\right)^{2}+\left(Z-Z^{\prime}\right)^{2}+\epsilon\right]}

(9)

which can be derived through calculating the vector potential generated by an axisymmetric current ring with unit current[20]. And $\epsilon$ is a small quantity added to avoid the singularity at the control surface caused by the complete elliptic integral of the first kind $K$ .

Matching on last closed field surface: To guarantee the continuity of poloidal magnetic flux $\psi$ between the open and the closed field line regions, the values of $\psi$ from both regions could be matched on a surface close to the LCFS at chosen matching points $j$

\psi_{pj}+\psi_{cj}=\left(1-\eta\right)\psi_{X}

(10)

where $\psi_{X}$ is the poloidal magnetic flux at the X point on separatrix, $\psi=\left(1-\eta\right)\psi_{X}$ is the flux on the matching surface, $\eta$ is a small number chosen to avoid singularities caused by discontinuities in the tangential slopes of the magnetic field lines at the X points[14]. For the purpose of plasma 2D shaping control, the extent that the exact configuration setting of the external poloidal field coils needs to be specified may be determined by the set of chosen locations where the matching conditions are required to be satisfied. From the expression of $\psi_{cj}$ , Eq. (7-9), one can see that the dependence of matching conditions on external coil current is linear, whereas their dependence on external coil location is nonlinear. Thus, it would be more straightforward and practical to solve for the set of external coil current at the fixed coil locations in order to satisfy the set of given matching conditions, which is often the case in reality.

In this way, we obtain the complete equilibrium solution to the G-S equation in Eq. (1) where the corresponding configuration of the external poloidal field coils can maintain the prescribed plasma profiles and 2D shape at the separatrix are also determined.

An example case and error evaluation: An up-down symmetric tokamak equilibrium with separatrix, which can be minimally specified with the first four items expansion terms of the homogeneous solutions along with a special solution of the G-S equation (1), is shown here as an example of the solution method described in the previous section:

\begin{array}[]{rl}\psi_{0}&=\frac{a}{8}R^{4}+\frac{b}{2}\left(R^{2}\ln R-\frac{R^{2}}{2}\right)\\ \psi_{1}&=1,\qquad\qquad\qquad\psi_{2}=R^{2}\\ \psi_{3}&=Z^{2}-R^{2}\ln R,\qquad\psi_{4}=R^{4}-4R^{2}Z^{2}\end{array}

(11)

In this case, the four unknown coefficients $c_{1}-c_{4}$ in Eq. (3) can be determined by the locations of the X points $(R_{X},\pm Z_{X})$ and the magnetic axis $\left(R_{0},0\right)$ together with the assumption that the poloidal flux function $\psi=0$ at the magnetic axis. This yields the corresponding constraints:

\begin{array}[]{rl}\frac{\partial\psi}{\partial R}&=0,\quad\frac{\partial\psi}{\partial Z}=0,\qquad\mbox{at X points }(R_{X},\pm Z_{X})\\ \hfill\hfil\\ \frac{\partial\psi}{\partial R}&=0,\quad\psi=0,\qquad\mbox{at magnetic axis }\left(R_{0},0\right)\\ \hfill\hfil\end{array}

(12)

The constraint that the derivative of $\psi$ on $Z$ at the magnetic axis equals zero is naturally satisfied, owing to the symmetry of $\psi$ in $Z$ . The constraints in Eq.(12) lead to a set of linear inhomogeneous algebraic equations for the unknown $c_{i}$ which uniquely specify the 2D shape of the magnetic flux surface in the closed field line region. And the external coil setting designed to maintain the 2D shape of the closed field line region can be determined by the matching conditions in Eq. (10), which would reduce to essentially a set of linear equations for the external coil current $I_{i}$ at fixed coil locations.

For example, the plasma current and pressure are specified by $a=1.2,b=-1.0$ . The shape parameters are chosen as $R_{0}=1.0,R_{X}=0.95,Z_{X}=\pm 0.06$ . Corresponding to the up-down symmetry of this particular equilibrium, a number of external coils at fixed locations with up-down symmetry are also prescribed. In this way, we only need to solve half number of the external coil currents. The exact number of the external coils is determined by the number of match points selected. For the Green function integral in Eq. (6), the two-dimensional trapezoidal rule is adopted on a $100\times 100$ uniform rectangular grid.

Using the method outlined in previous sections, we obtain the Solov’ev equilibrium with a vacuum region outside the separatrix (Fig. 1) and the corresponding external coil configuration. In comparison to the conventional Solov’ev equilibrium without a vacuum region, the contours of $\psi$ are the same in the closed field line region, and the shape of the separatrix is nearly identical subject to the numerical error. The relative error of the vacuum solution on the separatrix can be defined by $\Lambda=\frac{\psi_{X}^{vac}(\theta)-\psi_{X}}{\psi_{X}}$ , where $\psi_{X}$ is the Solov’ev solution on the separatrix and $\psi_{X}^{vac}(\theta)$ is the numerical vacuum solution along the separatrix and dependent on the poloidal angle $\theta$ at the given $\psi_{X}$ on the separatrix. As shown in the upper figure of Fig. 2, the maximum difference is at the X points, which is caused by the discontinuity of the tangential slope of the magnetic field line at the X points. Nonetheless, as seen in the lower figure of Fig. 2, $\psi$ is continuous at separatrix. Outside the separatrix, the contour of $\psi$ differs from the Solov’ev solution, as also indicated by the different external angles of field lines at the X points. The error $\Lambda$ depends on the matching parameter $\eta$ and the grid size. Keeping the grid size fixed, the relative error increases with $\eta$ , as shown in Fig. 3. For several fixed given $\eta$ , the numerical convergence in terms of $\Lambda$ with grid size is also confirmed.

2D shaping by external coils: In the above method, only a set of linear equations for the coil currents need to be solved for a set of prescribed coil locations, which suggests that the method may be used as an efficient 2D shaping control scheme using external coils. To demonstrate this, we solve for example the external coil settings for several Solov’ev equilibria with vacuum regions and the designed features of 2D shape, such as varying degrees of up-down asymmetry, elongation, and triangularity.

Up-down asymmetry: For an up-down asymmetric equilibrium, the positions of the magnetic axis and two different X points can be chosen to determine the plasma shape. As a result of the choice, we will need to select at least three additional expansions terms of the homogeneous solutions of Eq. (1).

\begin{array}[]{rl}\psi_{5}&=Z,\qquad\qquad\qquad\psi_{6}=R^{2}Z,\\ \psi_{7}&=Z^{3}-3ZR^{2}\ln R\end{array}

(13)

Now, seven constraints are required to solve the seven unknown coefficients $c_{i}$ .

\begin{array}[]{rl}\frac{\partial\psi}{\partial R}&=0,\quad\frac{\partial\psi}{\partial Z}=0,\qquad\mbox{two X points}~{}(R_{x1},Z_{x1}),~{}(R_{x2},Z_{x2})\\ \hfill\hfil\\ \frac{\partial\psi}{\partial R}&=0,\quad\frac{\partial\psi}{\partial Z}=0,\quad\mbox{magnetic axis}~{}(R_{0},0)\\ \hfill\hfil\\ \psi&=0,\qquad\mbox{magnetic axis}~{}(R_{0},0)\end{array}

(14)

The constraint on the derivative of $\psi$ on Z at the magnetic axis is required now because it is not satisfied automatically due to the loss of symmetry in Z. Similarly, the currents of all coils, instead of half, must be calculated. For the same external coil locations, and the plasma current and pressure profiles in the up-down symmetric Solov’ev equilibrium specified in Eq. (11), the conventional Solov’ev equilibrium and the Solov’ev equilibrium with a vacuum region with the same plasma shape are shown in Fig. 4 respectively.

Elongation: For an up-down symmetry equilibrium, the position of the X points can be used to change the elongation. Using the same constraints as in Eq. (12) and different positions of X points, an equilibrium with different elongation is obtained which is shown in Fig. 5 along with the classical Solov’ev equilibrium without vacuum for comparison.

Triangularity: It can be shown that if only the first four homogeneous solutions of the G-S equation, stated in Eq. (11), are chosen, the triangularity of Solov’ev equilibrium with X points is equal to 1. To vary the triangularity of equilibrium, another two expansions terms of the homogeneous solutions of the G-S equation are added

\begin{array}[]{rl}\psi_{5}&=2Z^{4}-9Z^{2}R^{2}-12Z^{2}R^{2}\ln R+3R^{4}\ln R\\ \psi_{6}&=R^{6}-12R^{4}Z^{2}+8R^{2}Z^{4}\end{array}

(15)

Two corresponding constraints are required to determine these two unknown coefficients, such as the requirement that $\psi$ at the inner and outer points of LCFS be the same as $\psi$ at the X point

\begin{array}[]{rl}\psi_{inner}&=\psi_{X},\qquad\mbox{inner points }(R_{\text{inner}},0)\\ \psi_{outer}&=\psi_{X},\qquad\mbox{outer points }(Z_{\text{inner}},0)\end{array}

(16)

The Solov’ev equilibrium thus solved is shown in Fig. 6 along with the conventional Solov’ev equilibrium with the same plasma current, pressure profiles, and plasma shape.

In summary, a general method for constructing an equilibrium with any specified 2D shape surrounded by a vacuum region using external poloidal field coils is presented, which allows the determination of the direct relation between the plasma shaping and the corresponding external poloidal field coil setting. And this method is demonstrated to be effective for specifying the elongation, triangularity, and up-down symmetry of Solov’ev equilibrium with a vacuum region by adjusting the current of a set of external coils at fixed locations, for example. However, this method is not limited to the Solov’ev equilibrium. We plan on implementing the method into 2D equilibrium solvers, such as NIMEQ[21, 22] and CHEASE[23], to obtain a desirable equilibrium with designed plasma shape using external coils and extending this method to the calculation of external coil setting for more general and realistic tokamak equilibrium with any prescribed 2D shape, which may be eventually developed into a practical 2D shaping control scheme.

Acknowledgements.

The author Jiaxing Liu thanks Mr. Rui Han (University of Science and Technology of China) for helpful discussions. This work was supported by the National Key Research and Development Program of China Grant No. 2019YFE03050004, the Fundamental Research Funds for the Central Universities at Huazhong University of Science and Technology Grant No. 2019kfyXJJS193, the National Natural Science Foundation of China Grant No. 51821005, and the U.S. Department of Energy Grant Nos. DE-FG02-86ER53218 and DE-SC0018001. The computing work in this paper is supported by the Public Service Platform of High Performance Computing by Network and Computing Center of HUST.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

Wahlberg [2007] C. Wahlberg, “Effect of a noncircular plasma cross section on the ideal (1, 1) mode in low-shear tokamaks,” Physics of Plasmas 14, 042506 (2007).
Aiba et al. [2007] N. Aiba, S. Tokuda, T. Takizuka, G. Kurita, and T. Ozeki, “Effects of ‘sharpness’ of the plasma cross-section on the MHD stability of tokamak edge plasmas,” Nuclear Fusion 47, 297–304 (2007).
Risk [1981] H. M. Risk, “MHD stability of toroidal plasma with modulated curvature planar magnetic axis,” Plasma Physics 23, 385–397 (1981).
Grauer and Rebhan [1984] R. Grauer and E. Rebhan, “Analysis of ideal mhd equilibrium and axisymmetric stability for finite aspect ratio tokamaks with elliptic cross-section and flat current profile,” Journal of Plasma Physics 32, 99–117 (1984).
Solov’ev [1968] L. S. Solov’ev, “The theory of hydromagnetic stability of toroidal plasma configurations,” Soviet Physics JETP 26, 400–407 (1968).
Zheng, Wootton, and Solano [1996] S. B. Zheng, A. J. Wootton, and E. R. Solano, “Analytical tokamak equilibrium for shaped plasmas,” Physics of Plasmas 3, 1176–1178 (1996).
Guazzotto and Freidberg [2007] L. Guazzotto and J. P. Freidberg, “A family of analytic equilibrium solutions for the Grad–Shafranov equation,” Physics of Plasmas 14, 112508 (2007).
Cerfon and Freidberg [2010] A. J. Cerfon and J. P. Freidberg, “"One size fits all" analytic solutions to the grad–shafranov equation,” Physics of Plasmas 17, 032502 (2010).
Yuan et al. [2013] Q. Yuan, B. Xiao, Z. Luo, M. Walker, A. Welander, A. Hyatt, J. Qian, R. Zhang, D. Humphreys, J. Leuer, R. Johnson, B. Penaflor, and D. Mueller, “Plasma current, position and shape feedback control on EAST,” Nuclear Fusion 53, 043009 (2013).
Hyatt et al. [2013] A. W. Hyatt, D. A. Humphreys, A. S. Welander, N. W. Eidietis, J. R. Ferron, J. M. Hanson, R. D. Johnson, E. Kolemen, M. A. Lanctot, D. Moreau, B. G. Penaflor, E. Schuster, F. Turco, M. L. Walker, R. Coon, and J. Qian, “Designing, constructing and using plasma control system algorithms on DIII-D,” in 2013 IEEE 25th Symposium on Fusion Engineering (SOFE) (2013) pp. 1–6.
Hofmann and Jardin [1990] F. Hofmann and S. Jardin, “Plasma shape and position control in highly elongated tokamaks,” Nuclear Fusion 30, 2013–2022 (1990).
Albanese, Coccorese, and Rubinacci [2011] R. Albanese, E. Coccorese, and G. Rubinacci, “Plasma modelling for the control of vertical instabilities in tokamaks,” Nuclear Fusion 29, 1013 (2011).
Pironti and Portone [1998] A. Pironti and A. Portone, “Optimal choice of the geometrical descriptors for tokamak plasma shape control,” Fusion Engineering and Design 43, 115–127 (1998).
Xu and Fitzpatrick [2019] T. Xu and R. Fitzpatrick, “Vacuum solution for Solov’ev’s equilibrium configuration in tokamaks,” Nuclear Fusion 59, 064002 (2019).
Shafranov and Yurchenko [2014] V. Shafranov and E. Yurchenko, “Condition for flute instability of a toroidal-geometry plasma,” Journal of Experimental and Theoretical Physics 26, 682–686 (2014).
Freidberg [2014] J. P. Freidberg, Ideal MHD (Cambridge University Press, 2014).
[17] L. E. Zakharov and V. D. Shafranov, “Equilibrium of a toroidal plasma with noncircular cross section,” Sov. Phys.-Tech. Phys. (Engl. Transl.), v. 18, no. 2, pp. 151-156 .
Reusch and Neilson [1986] M. Reusch and G. Neilson, “Toroidally symmetric polynomial multipole solutions of the vector laplace equation,” Journal of Computational Physics 64, 416–432 (1986).
Jardin [2010] S. C. Jardin, Computational Methods in Plasma Physics (CRC press, 2010).
Jackson [1998] J. D. Jackson, Classical Electrodynamics (Wiley, 1998).
Howell and Sovinec [2014] E. C. Howell and C. R. Sovinec, “Solving the Grad–Shafranov equation with spectral elements,” Computer Physics Communications 185, 1415–1421 (2014).
Li and Zhu [2021] H. Li and P. Zhu, “Solving the grad–shafranov equation using spectral elements for tokamak equilibrium with toroidal rotation,” Computer Physics Communications 260, 107264 (2021).
Lütjens, Bondeson, and Sauter [1996] H. Lütjens, A. Bondeson, and O. Sauter, “The CHEASE code for toroidal MHD equilibria,” Computer Physics Communications 97, 219–260 (1996).

Refer to caption — Figure 1: Contour of $\psi$ in Solov’ev solution with (lower) and without (upper) vacuum region, where $a=1.2$ , $b=-1.0$ , $R_{0}=1.0$ , $R_{X}=0.95$ , $Z_{X}=\pm 0.06$ as in Eq. (11) and Eq. (12)