Exact Solov'ev equilibrium with an arbitrary boundary

TL;DR

This paper presents a method to compute exact Solov'ev plasma equilibria with arbitrary boundaries, allowing high accuracy and flexibility in boundary shape, including X-points, using a constrained least-squares approach.

Contribution

It introduces a novel constrained least-squares technique for calculating exact plasma equilibria with arbitrary boundaries, enhancing accuracy and boundary shape flexibility.

Findings

Successfully computed highly-shaped equilibria

Flexible boundary representation with many constraints

Demonstrated accuracy with arbitrary boundary shapes

Abstract

Exact Solov'ev equilibria for arbitrary plasma cross-sections are calculated using a constrained least-squares method. The boundary, with or without $X$-points, can be specified with an arbitrarily large number of constraints to ensure an accurate representation. Thus, the order of the polynomial basis functions in the homogeneous solution of the Grad-Shafranov equation becomes an independent parameter determined only by the accuracy requirements of the overall solution. Examples of exact, highly-shaped equilibria are presented.

Equations

$$-{\Delta }^{\ast }\psi ={\mu }_0{R}^2{p}^{\prime }(\psi )+F{F}^{\prime }(\psi ),$$

$${\Delta }^{\ast }\psi \equiv R^2\nabla \cdot \frac{1}{R^2}\nabla \psi =R\frac{\partial }{\partial R}\frac{1}{R}\frac{\partial \psi }{\partial R}+\frac{{\partial }^2\psi }{\partial Z^2}.$$

$${\Delta }^{\ast }\psi =A{R}_0^2+C{R}^2,$$

$$\psi =\frac{1}{2}(A{R}_0^2+a_0{R}^2)Z^2+\frac{1}{8}(C-a_0)(R^2-R_0^2{)}^2,$$

$${\Delta }^{\ast }\psi =0,$$

$${\Delta }^{\ast }\psi =A+C{R}^2,$$

$${\psi }_p=\frac{C}{8}{R}^4+\frac{A}{2}{R}^2\ln R.$$

$${\psi }_h=\sum _{i=0}{f}_i(R)Z^i.$$

$$R\frac{\partial }{\partial R}\frac{1}{R}\frac{\partial }{\partial R}{f}_{i-2}(R)=-i(i-1)f_i(R),i\ge 2.$$

$${\psi }_h(R,Z)=\sum _{i=0}^{I-1}{c}_i{P}_i(R,Z)+\sum _{i=1}^I{d}_i{Q}_i(R,Z),$$

$$\{c_0,c_1,\dots ,c_{I-1},d_1,d_2,\dots ,d_I\},$$

$$R(\alpha )$$

$$Z(\alpha )$$

$$\delta (\alpha )$$

$$\kappa (\alpha )$$

$$\psi (R_j,Z_j)$$

$$(R_j,Z_j)$$

$$S$$

$$\sum _{i=0}^{I-1}{c}_i⟨P_i,P_k⟩+\sum _{i=1}^I{d}_i⟨Q_i,P_k⟩$$

$$\sum _{i=0}^{I-1}{c}_i⟨P_i,Q_l⟩+\sum _{i=1}^I{d}_i⟨Q_i,Q_l⟩$$

$$⟨p,q⟩\equiv \sum _{j=0}^{M-1}w(R_j,Z_j)p(R_j,Z_j)q(R_j,Z_j).$$

$$\psi (R_X,Z_X)$$

$${\frac{\partial \psi }{\partial R}|}_{(R_X,Z_X)}$$

$$\mathcal{L}(\mathbf{c},\mathbf{d},\mathbf{\lambda })=S+{\lambda }_1{\frac{\partial \psi }{\partial R}|}_{(R_X,Z_X)}+{\lambda }_2{\frac{\partial \psi }{\partial Z}|}_{(R_X,Z_X)}+{\lambda }_3\psi (R_X,Z_X),$$

$$\frac{\partial \mathcal{L}}{\partial c_k}$$

$$\frac{\partial \mathcal{L}}{\partial d_l}$$

$$\frac{\partial \mathcal{L}}{\partial {\lambda }_1}$$

$$\frac{\partial \mathcal{L}}{\partial {\lambda }_2}$$

$$\frac{\partial \mathcal{L}}{\partial {\lambda }_3}$$

$$w(R_j,Z_j)=1-e^{-\{[(R_j-R_X{)}^2+(Z_j-Z_X{)}^2]/\Delta w^2\}},$$

$$P_0$$

$$P_1$$

$$P_2$$

$$P_3$$

$$P_4$$

$$P_5$$

$$P_6$$

$$P_7$$

Full text

I Introduction

An axisymmetric equilibrium, starting point for many linear stability and nonlinear initial-value calculations for toroidal devices, is described by the Grad-Shafranov equation Grad and Rubin (1958); Shafranov (1957),

[TABLE]

where primes denote a derivative with respect to $\psi$. In cylindrical $(R,Z,\zeta)$ coordinates, where the toroidal angle $\zeta$ is in the symmetry direction, the Laplacian-like operator $\Delta^{*}$ has the form

[TABLE]

The equilibrium is determined by specifying the two source functions on the right, $p(\psi),F(\psi)$, for the plasma pressure and poloidal current, respectively.

Although the original references are not easily accessible anymore, derivation of Eq. 1 can be found in many textbooks (see, for example JardinJardin (2010), FreidbergFreidberg (2014)). The review article by Takeda and TokudaTakeda and Tokuda (1991), although somewhat dated, is also still useful. In “fixed-boundary” equilibria, a boundary curve $C_{b}$ in the $(R,Z)$ (poloidal) plane is also specified a priori. Then Eq. 1 is solved with a Dirichlet boundary condition on $C_{b}$, typically $\psi=0.$ In “free-boundary” equilibria, a distribution of external currents in the poloidal field (PF) coils is given and the plasma boundary is determined as part of the solution. Here, effectively a boundary condition at infinity needs to be imposed, $\psi\rightarrow 0$. Usually this problem in the infinite domain is tackled using an artificial computational boundary, $C_{cb}$, a finite distance away from the plasma. Recently we have introduced a new method that directly and efficiently solves the free-boundary problem in the infinite domainHan et al. (2019) without a need for $C_{cb}$. Note that besides $p(\psi)$ and $F(\psi)$, a more general pair of source functions can also be used as inputLütjens et al. (1996). In fact, widely-used EFIT-series of equilibrium codesLao et al. (2005) can calculate equilibria consistent with a wide range of experimental data from kinetic and magnetic diagnostics.

A realistic equilibrium that satisfies a number of constraints, for example, for the plasma pressure and safety factor profiles, generally requires $p$ and $F$ to be nonlinear functions of their argument. Therefore, Eq. 1 is typically a nonlinear equation for the flux function $\psi$, and it is solved numerically. Non-numerical solutions, however, are quite useful in analytic equilibrium and stability studies; they are also very useful in checking and benchmarking of numerical codes. For these reasons, exact analytic solutions of the Grad-Shafranov equation have been sought by a number of workers over the yearsMcCarthy (1999); Atanasiu et al. (2004); Guazzotto and Freidberg (2007). The most well-known and useful of these, however, is a class of exact solutions that have come to be known as the Solov’ev equilibriaSolov’ev (1967), which assume $p$ and $F^{2}$ are linear functions of $\psi$. In Eq. 1 letting $\mu_{0}p^{\prime}=-C,~{}FF^{\prime}=-AR_{0}^{2},$ where $A,C$ are constants and $R_{0}$ is the major radius of the torus, leads to the inhomogenous but linear equation

[TABLE]

with the exact Solov’ev solution

[TABLE]

where $a_{0}$ is another constant. The three arbitrary constants, $a_{0},A,C$, can be adjusted to change, for instance, the separatrix geometrySolov’ev (1967). However, there are not enough degrees of freedom in the solution to specify basic equilibrium parameters such as the plasma $\beta=2\mu_{0}p/B^{2}$, safety factor $q$, and at the same time ensure that the solution conforms to a particular plasma boundary.

Solutions of the homogeneous equation,

[TABLE]

can provide the needed extra degrees of freedom, since a general solution to Eq. 3 can be written as $\psi=\psi_{p}+\psi_{h}$, where $\psi_{p}$ is a particular solution of Eq. 3, and $\psi_{h}$ satisfies Eq. 5. A complete set of multipole solutions of Eq. 5 were found by Reusch and Nielson and applied to PDX vacuum field calculationsReusch and Nielson (1986). A similar polynomial expansion in powers of $R,Z$ was used by Zheng et al.Zheng et al. (1996) to find shaped but up-down symmetric, exact analytic solutions to Eq. 3. These were later generalized by Cerfon and FreidbergCerfon and Freidberg (2010) to geometries without up-down symetry, including diverted configurations with one or more $X$-points.

In both ZhengZheng et al. (1996) and CerfonCerfon and Freidberg (2010), the number of constraints, $M$, used to specify the boundary was limited by the degrees of freedom, $N$, in the homogeneous solution. Thus, a more detailed boundary definition required a larger $N$, which translated into a higher order polynomial expansion for the homogeneous solution $\psi_{h}.$ In this work we show how the requirements of the boundary definition can be decoupled from the polynomial expansion of $\psi_{h}$ using constrained least-squares. Now the boundary, both “fixed” and “free,” can be defined essentially arbitrarily while still maintaining a fixed and small number of terms in the expansion of $\psi_{h}.$ The method is explained in Sect. II with various examples of exact tokamak equilibria. Sect. III summarizes our results.

II Exact solutions with an arbitrary boundary

It is useful to put Eq. 3 in nondimensional form by transforming the variables. After letting $\psi=(B_{0}R_{0}^{2}){\widetilde{\psi}}$, $R=R_{0}{\widetilde{R}}$, $Z=R_{0}{\widetilde{Z}}$, the equation becomes

[TABLE]

where we dropped the tildes in normalized variables. Typically $B_{0}$ is the toroidal field strength at the major radius $R_{0}.$ A $Z$-independent particular solution is given by

[TABLE]

The constants $A,C$ can be used to specify, for example, the plasma $\beta$ and the kink safety factor $q^{*}=\epsilon B_{0}/B_{p},$ where $\epsilon\equiv a/R_{0}$ is the inverse aspect ratio, and $B_{p}$ is the poloidal field at the edge. As stated earlier, to satisfy the constraints that help define the boundary, we will need the degrees of freedom provided by the homogeneous solution.

II.1 Homogeneous solution

Following Zheng et al.Zheng et al. (1996) we write the homogeneous solution in the form

[TABLE]

Substituting in Eq. 5 leads to the recursion relation

[TABLE]

Clearly the even and odd terms decouple. Letting $f_{i}(R)=0$ for $i\geq I$ and integrating yields a solution with $N=2I$ constants of integration (we will assume $I$ is even). Then the homogeneous solution can be written as a sum of polynomial basis functions:

[TABLE]

where $P_{i}$ and $Q_{i}$ have even and odd order, respectively. For up-down symmetric boundaries, the second sum is not needed. The resulting polynomials for $I=10$, obtained using MathematicaWolfram Research, Inc. , are shown in the Appendix. The set of arbitrary coefficients,

[TABLE]

will be determined by the boundary constraints to be explained in the next section. Since the coefficients are arbitrary, each polynomial basis function has to satisfy the homogeneous equation identically, which can be shown to be true by direct substitution. Below we first consider “limited” or fixed-boundary tokamak equilibria without $X$-points.

II.2 Exact fixed-boundary equilibria

A common parametrization of shaped plasma boundaries has the form

[TABLE]

where $(R,Z)$ have been normalized by $R_{0}$. The triangularity $\delta$ and elongation $\kappa$ are defined in terms of their upper and lower values as:

[TABLE]

The geometric meanings of various terms are shown in Fig. 1. Note that the parameter $\alpha$ differs from the angle $\theta=\arctan(Z/(R-R_{0})).$ The actual triangularities, which measure the shift of the highest and lowest points of the boundary with respect to the geometric center, are $\Delta_{U}=\sin\delta_{U}\simeq\delta_{U}$, and $\Delta_{L}=\sin\delta_{L}\simeq\delta_{L}.$

A useful set of boundary constraints uses an equally spaced points on the boundary curve, as shown in Fig. 1, and take the following form:

[TABLE]

where $\psi_{b}$ is the boundary value; we will set $\psi_{b}=0.$

If the number of constraints exactly matches the degrees of freedom in the homogeneous solution, i.e., if $M=N$, then the set of equations in Eq. 14 can be solved easily for the $N=2I$ coefficients ($c_{i},d_{i})$ of Eq. 10, a method adopted by the earlier worksZheng et al. (1996); Cerfon and Freidberg (2010). However, there are obvious advantageous to having $M\gg N$. A highly-shaped boundary may require many more points to specify accurately than the degrees of freedom available in the homogeneous solution. But with $M>N$, the system is overdetermined and the linear system resulting from Eq. 14 has no solution.

II.3 Least squares solution

We can still find a solution for $M>N$ if we do not demand that all points be located exactly on the boundary, i.e., we stop requiring that Eq. 14 be satisfied exactly. Instead we seek the “best” solution that minimizes, for example, a positive-definite error term. In the least-squares method we minimize the sum of the “residuals,” defined by

[TABLE]

where $(R_{j},Z_{j}),~{}j=0,1,\ldots,M-1$ are points on the boundary parametrically determined by Eq. 12, and $w_{j}\equiv w(R_{j},Z_{j})$ is a weight function that can be used to emphasize/de-emphasize some subset of the boundary points. It is not needed for fixed-boundary equilibria and is set to unity.

Letting $\psi_{b}=0$ and seeking a minimum by setting $\partial S/\partial c_{k}=0,~{}\partial S/\partial d_{l}=0$ lead to

[TABLE]

where we defined

[TABLE]

Together Eqs. 16, 17 form a set of $N=2I$ linear equations for the $N$ coefficients $c_{i},d_{i}$ that can be easily inverted. Again, if the boundary curve is symmetric about the mid-plane, we can set $d_{i}=0$ above and solve only for the $N=I$ coefficients $c_{i}$. In this case the parameter $\alpha_{j}$ in Eq. 14 can be redefined to cover, for instance, only the lower half of the boundary.

Examples of highly-shaped exact equilibria obtained with this method are shown in Fig. 2. For each one, $M$ number of boundary constraints are specified at the indicated points (blue dots in the figures) using Eq. 14 and the parametric curve of Eq. 12. For the asymmetric case in (a), both even and odd basis functions are used in the expansion of $\psi_{h}$; thus, there are $N=2I$ degrees of freedom in the solution, where $I=6$ is the degree of the highest-order even basis function in Eq. 10. For the up-down symmetric cases ((b) and (c)) only points in the lower half-plane are used. For these, the odd basis functions are not needed; thus they both have $N=I=6$. For all equilibria we use $M=4N$, so the resulting linear system of equations for the coefficients would be overdetermined without the least squares approach. In the figures, the actual calculated boundary is shown with a dashed (black) curve. Clearly the specified boundary (dots) and the actual boundary curve agree very well, although the least squares solution to the boundary constraints is not exact but only minimizes the sum of residuals in Eq. 15. Note that, although not shown here, all solutions with $I=4$ (maximum order of even polynomials) also produce acceptable results for these highly-shaped fixed-boundary equilibria.

II.4 Exact equilibria with $X$-points

An $X$-point, where the poloidal field vanishes, introduces three additional constraints that the solution $\psi=\psi_{p}+\psi_{h}$ has to satisfy:

[TABLE]

where $(R_{X},Z_{X})$ is the $X$-point location. Now the separatrix becomes the plasma boundary; for simplicity we still choose $\psi_{b}\equiv\psi_{sep}=0.$

These three constraints, however, are of a different nature than the boundary constraints in Eq. 14; they have to be satisfied exactly (otherwise there is no $X$-point), whereas those in Eq. 14 can be satisfied in a “least-squares” sense–we can tolerate the computational boundary agreeing with the specified boundary only approximately to accommodate the $X$-point(s). Thus, we now have a “constrained least-squares” problem, which we solve by introducing a “Lagrangian ${\cal L}$” and the Lagrange multipliers $\lambda_{j},~{}j=1,2,3$:

[TABLE]

where the sum of residuals $S$ is still given by Eq. 15. The unknown coefficient vectors $\bm{c}=\{c_{0},c_{1},\ldots,c_{I-1}\}$, $\bm{d}=\{d_{1},d_{2},\ldots,d_{I}\}$ and $\bm{\lambda}=\{\lambda_{1},\lambda_{2},\lambda_{3}\}$ form the new set of $N+3$ unknowns. They are determined by an augmented set of equations:

[TABLE]

These equations guarantee that, among the set of solutions $\{{\bm{c}},{\bm{d}}\}$ that satisfy the $X$-point constraints in Eqs. 23-25 exactly, we choose the one that minimizes the residual sum $S$. Recall that $N=2I$ for up-down asymmetric equilibria. For symmetric ones, the coefficients $\bm{d}$ for the odd basis functions and Eq. 22 are not needed, and $N=I$ (see Eq. 10).

We find that de-emphasizing the points in the sum of residuals $S$ near an $X$-point leads to better results. Thus, if there is a field-null, we choose a Gaussian weight function in Eq. 15 of the form

[TABLE]

where $\Delta w$ is an appropriately-chosen width and $(R_{X},Z_{X})$ is the location of the null-point. Note that $w\rightarrow 0$ near the $X$-point so that those points are less-constrained to be on the parametrically-defined boundary curve of Eq. 12, which is necessary to accommodate the $X$-point.

Examples of exact equilibria with $X$-points–“free-boundary” versions of the fixed-boundary equilibria in Fig. 2–are shown in Fig. 3. We note that in all three cases, the computed boundary (dashed black curve) conforms to the specified boundary (represented by blue dots) very closely, except near the $X$-points where it necessarily deviates in order to satisfy the $X$-point constraints. Thus, even with one or more $X$-points, an exact equilibrium with a parametrically-specified boundary can be calculated using constrained least-squares method. However, all three examples require higher-order polynomial expansions than their fixed-boundary counterparts in Fig. 2. The lower single null (LSN) configuration in (a) used an expansion with 16 degrees of freedom, or equivalently, the highest-order even basis function had order $I=8.$ Similarly, $I=8$ for (b), the elongated double-null (DN). The DN configuration with the negative triangularity in (c) required $I=12.$

III Summary and discussion

We presented a method to calculate exact Solov’ev equilibria with or without $X$-points for arbitrary plasma cross-sections. Using a least-squares approach to enforce the boundary constraints essentially decouples the accuracy requirements of the boundary specification from the order of the polynomial expansion of the homogeneous solution for the Grad-Shafranov equation. Thus, a much larger number of boundary constraints can be used than the number of degrees of freedom available in the homogeneous solution, ensuring that the numerically-obtained boundary conforms to the specified boundary with high fidelity. Usefulness of the method was illustrated with examples of fixed-boundary equilibria with highly-shaped geometries and their generalization to configurations with one or more $X$-points.

These exact solutions, which require nothing more than the inversion of a relatively small linear system of equations, can be easily performed with a symbolic/numerical algebra system such as Mathematica or Matlab. They should be highly useful, for example, in benchmarking new equilibrium codes. They may also prove useful in semi-analytic equilibrium and stability calculations.

Acknowledgements

This work was supported by MSIP, the Korean Ministry of Science, ICT and Future Planning, through the KSTAR project. Most of the symbolic and numerical calculations were performed with Mathematica.

Appendix

The polynomials $P_{i}(R,Z),Q_{i}(R,Z)$ of Eq. 10 for $I=10$, obtained with MathematicaWolfram Research, Inc. :

[TABLE]

No attempt has been made to simplify or normalize the polynomials. They are presented in the form found by Mathematica. Note that $P_{I-1}$ has order $I$, while $Q_{I}$ has order $I+1$.