Simple, general, realistic, robust, analytic tokamak equilibria. Part 1. Limiter and divertor tokamaks

Abstract

Tokamak equilibria have been derived that are analytic solutions to the Grad–Shafranov equation. This paper, Part 1, describes a wide range of such equilibria including smooth limiter surfaces, double- and single-null divertor surfaces, arbitrary aspect ratio, elongation, triangularity and beta. Part 2 generalizes the analysis to include edge pedestals and toroidal flow.

1. Introduction

The title of the paper tells the whole story. We have derived simple, general, realistic, robust, analytic solutions to the Grad–Shafranov equation describing tokamak equilibria. What is meant by all of these adjectives? ‘Simple’ refers to the fact that the equilibria require only a few, intuitively simple-to-visualize, terms. ‘General’ indicates that the equilibria are valid for a wide range of configurations, including smooth limiter surfaces, double-null divertor surfaces, single-null divertor surfaces, finite aspect ratio including spherical tokamaks, finite elongation, finite positive and negative triangularity and small, medium and large beta. ‘Realistic’ implies that our profiles are continuous and monotonic with the pressure, pressure gradient and toroidal current density smoothly vanishing at the plasma edge. Finite edge pedestals, an edge contribution due to the bootstrap current, and toroidal flow are generalizations included in Part 2. ‘Robust’ refers to the fact that only four (or six for asymmetric configurations) input parameters are required to obtain a solution for the flux surfaces – three (or five) geometrical parameters and one simple physically intuitive parameter. The model leads to well-behaved solutions every time – no delicate choosing of any parameters. Lastly, ‘analytic’ indicates that our solutions satisfy the exact, unexpanded Grad–Shafranov equation, and are expressed in terms of known functions. The final results are simple analytic expressions for the flux function and importantly its first and second derivatives.

To put the present work in perspective, it is useful to compare with earlier related research. Solov'ev calculated the first purely analytic equilibrium (Solov'ev 1968). The well-known ‘Solov'ev profile’ is characterized by choosing $\textrm{d}p/\textrm{d}\psi$ and $F\,\textrm{d}F/\textrm{d}\psi $, the free functions on the right-hand side of the Grad–Shafranov equation, as constants. The resulting solutions capture many useful features of tokamak equilibria, but are unrealistic in that important physical profiles have large jumps at the plasma edge. Specifically, in physical space, the toroidal and poloidal current densities, plus the plasma pressure gradient have edge jumps. The pressure itself is parabolic. Even so, it is an analytic equilibrium and over the years several generalizations have been developed (Zheng, Wooten & Solano 1996; Weening 2000; Shi 2005; Cerfon & Freidberg 2010) which allowed application to a much wider range of equilibrium plasma shapes, although still with the same jumps at the plasma edge.

The next step of progress (see e.g. Herrnegger 1972; Maschke 1972; De Menna 1977; McCarthy 1999; Atanasiu, Lackner & Miron 2004; Guazzotto & Freidberg 2007; Shi 2009) allows the right-hand-side functions, $\textrm{d}p/\textrm{d}\psi $ and $F\,\textrm{d}F/\textrm{d}\psi $, to be linear (or linear plus constant) in $\psi $ rather than purely constant. This has the important benefit of producing smooth, monotonic profiles for which the current densities, pressure and pressure gradient vanish at the plasma edge. The resulting solutions more closely model experimental profiles than the Solov'ev solutions.

In general, the solutions associated with the linear-in-$\psi $ free functions $\textrm{d}p/\textrm{d}\psi $ and $F\,\textrm{d}F/\textrm{d}\psi $ are complicated, involving a sum of terms obtained by separation of variables. The radial separation functions consist of unintuitive hypergeometric functions, Coulomb wave functions or Whittaker functions of imaginary order and imaginary argument. Furthermore, the number of terms to be maintained plus the choice of separation constants is not robust, and in some cases delicate choices have to be made to obtain physically realistic solutions. Consequently, even though purely analytic solutions should be helpful in experimental comparisons, reactor design, diagnostic development, magnetohydrodynamic stability determination and benchmarking more general Grad–Shafranov solvers, the complexity and sensitivity of the solutions have hindered implementation from widely occurring.

This assessment applies to the references cited above, from Herrnegger (1972) to Shi (2009), including that of the present authors (Guazzotto & Freidberg 2007). While the corresponding research traces a continuing path of progress, describing increasingly realistic plasma shapes and experimental comparisons, there is still more work to be done.

In this context, the present work also focuses on analytic solutions in which the right-hand side of the Grad–Shafranov equation is a combination of linear and constant terms in $\psi $. However, an underutilized mathematic insight allows us to overcome the difficulties just described. In detail, we also use separation of variables and the resulting solutions can again be formally expressed in terms of Whittaker functions. Even so, ‘Whittaker’ is just a ‘name’ for the functions. Our solutions are rewritten in a much more convenient form in which each radial separation function is approximately a sine or cosine. This enables a sharp determination of (i) the number of terms (7 for symmetric and 12 for asymmetric systems) required for the solution as well as (ii) the choice plus number of separation constants (four for both cases). The details are described in the main text. The end result is a robust and general analytic solution to the Grad–Shafranov equation. A computer code has been developed, primarily to plot flux surfaces and carry out post-processing evaluation of various plasma parameters of interest. Typically, the eigenvalue is calculated to six-figure accuracy (although this can easily be increased), while the values of the flux and its first and second derivatives are accurate to 10–15 significant figures.

The code, written in MATLAB, is freely available by contacting the first author (L.G. at lzg0022@auburn.edu). On a laptop computer, it takes approximately 1 s to calculate and plot the flux surfaces and an additional 1 s for standard post-processing evaluations. The post-processing can take substantially longer, of the order of a minute, if a large number of very-high-resolution $q(\psi )$ surfaces are required. In its present form the code, which has not been optimized for speed, is not practical for real-time reconstruction and control. However, optimization could resolve this problem although the solutions would still be limited by the choice of free functions required for analytic solutions.

In terms of organization, the research is separated into two papers because of the large number of examples to be presented. The present paper, Part 1, focuses on standard tokamaks with smooth, double-null divertor or single-null divertor plasma surfaces. The plasma edge is idealized, with the pressure, pressure gradient and current density all vanishing on this surface. Part 2 treats the same plasma shapes, but generalizes the edge treatment. Edge jumps, which are a simple approximation to narrow smooth pedestals, in the pressure, pressure gradient and current density are allowed. Also, the edge-localized portion of the bootstrap current can be included. Toroidal flow is a final feature that is added.

In summary, we believe the work presented here represents a very useful step of progress along the path of obtaining realistic, analytic solutions to the Grad–Shafranov equation. Our hope is that the simplicity, accuracy and speed of the evaluation code will encourage researchers to make use of the results.

2. Formulation of the problem

2.1. Starting equations and boundary conditions

The tokamak equilibria of interest are described by the familiar Grad–Shafranov equation (Lüst & Schlüter 1957; Grad & Rubin 1958; Shafranov 1966):

(2.1)\begin{equation}R\frac{\partial }{{\partial R}}\left( {\frac{1}{R}\frac{{\partial \varPsi }}{{\partial R}}} \right) + \frac{{{\partial ^2}\varPsi }}{{\partial {Z^2}}} ={-} {\mu _0}{R^2}\frac{{\textrm{d}p}}{{\textrm{d}\varPsi }} - \frac{1}{2}\frac{{\textrm{d}{F^2}}}{{\textrm{d}\varPsi }},\end{equation}

with $\varPsi $ the poloidal flux function and $p(\varPsi )$, $F(\varPsi ) = R{B_\phi }$ the two free functions. As discussed, we shall assume in Part 1 that the free functions are chosen to be quadratic in $\varPsi $:

(2.2)\begin{equation}\left. \begin{gathered} {p(\varPsi ) = {p_0}{{\left( {\dfrac{\varPsi }{{{\varPsi_0}}}} \right)}^2},}\\ {{F^2}(\varPsi ) = R_0^2B_0^2\left[ {1 + \dfrac{{2\delta B}}{{{B_0}}}{{\left( {\dfrac{\varPsi }{{{\varPsi_0}}}} \right)}^2}} \right].} \end{gathered}\right\}\end{equation}

Here, ${\varPsi _0}$ is the flux on the magnetic axis, ${R_0}$ is the major radius, ${B_0}$ is the vacuum toroidal field at $R = {R_0}$, ${p_0}$ is the pressure on axis and $\delta B/{B_0}$ is a measure of the toroidal field diamagnetism on axis. Mathematically and physically, we require that each of these parameters except ${\varPsi _0}$ be specified as inputs. The value of ${\varPsi _0}$ is determined as an output of the calculation. Also, based on our definition of realistic profiles, we require that the plasma pressure, plasma pressure gradient and toroidal current density vanish on the plasma surface. This implies that the boundary condition on $\varPsi $ is

(2.3)\begin{equation}\varPsi (\textrm{surf}) = 0.\end{equation}

With the choice of free functions given above, the Grad–Shafranov equation becomes a linear partial differential equation in $\varPsi $:

(2.4)\begin{equation}R\frac{\partial }{{\partial R}}\left( {\frac{1}{R}\frac{{\partial \varPsi }}{{\partial R}}} \right) + \frac{{{\partial ^2}\varPsi }}{{\partial {Z^2}}} + \left( {\frac{{2{\mu_0}{p_0}}}{{\varPsi_0^2}}{R^2} + \frac{{2R_0^2{B_0}\delta B}}{{\varPsi_0^2}}} \right)\varPsi = 0.\end{equation}

2.2. The solution procedure

The procedure used in earlier studies to obtain solutions is to introduce an appropriate set of normalized variables and then solve the resulting equation by separation of variables. For example, define $R = K{u^{1/2}}$, $Z = Kv$, $K = {(\varPsi _0^2/2{\mu _0}{p_0})^{1/4}}$. The Grad–Shafranov equation reduces to

(2.5)\begin{equation}4u\frac{{{\partial ^2}\varPsi }}{{\partial {u^2}}} + \frac{{{\partial ^2}\varPsi }}{{\partial {v^2}}} + (u + \zeta )\varPsi = 0,\quad \zeta = \frac{{2R_0^2{B_0}\delta B{K^2}}}{{\varPsi _0^2}}.\end{equation}

Solutions to this equation can be found by summing a finite or infinite set of terms, each obtained by separation of variables:

(2.6)\begin{equation}\left. \begin{gathered} {\varPsi (u,v) = \sum {{U_n}(u){V_n}(v)} ,}\\ {\dfrac{{{\textrm{d}^2}{V_n}}}{{\textrm{d}{v^2}}} + {\zeta_n}{V_n} = 0,}\\ {\dfrac{{{\textrm{d}^2}{U_n}}}{{\textrm{d}{u^2}}} + \left( {\dfrac{1}{4} + \dfrac{{\zeta - {\zeta_n}}}{{4u}}} \right){U_n} = 0.} \end{gathered} \right\}\end{equation}

Here, ${\zeta _n}$ is the nth separation constant, which can be positive or negative. The solution is given by

(2.7)\begin{equation}\left. \begin{gathered} {{V_n}(v) = {c_1}\,\textrm{exp}(\textrm{i}\zeta_n^{1/2}v) + {c_2}\,\textrm{exp}( - \textrm{i}\zeta_n^{1/2}v),}\\ {{U_n}(u) = {c_3}{W_{\textrm{i}\kappa ,1/2}}(\textrm{i}u) + {c_4}{M_{\textrm{i}\kappa ,1/2}}(\textrm{i}u),}\end{gathered} \right\}\end{equation}

where $\kappa = (\zeta - {\zeta _n})/4$ and ${W_{\textrm{i}\kappa ,1/2}}(\textrm{i}u),\;{M_{\textrm{i}\kappa ,1/2}}(\textrm{i}u)$ are Whittaker functions.

The Whittaker solutions are a mixed blessing. On the positive side, the solutions are known analytic functions. This is important because many basic mathematical properties have been derived and are easily accessible in the literature. Also, they can be called in well-known mathematical programs such as MATLAB and Mathematica or solved for by direct numerical computation.

On the negative side, there is no guarantee that a sum of such terms will converge to a desired solution. The reason is that the natural convergence boundaries for the $R,\;Z$ separated coordinates are rectangles while plasma boundaries are quite different, typically an elongated D-shape. Thus, the radius of convergence of the series may not include the entire plasma shape. It is somewhat like trying to fit a round peg into a square hole. Another point is that the solutions are not very intuitive as they consist of Whittaker functions of imaginary order with imaginary argument. This lack of intuition makes it difficult to choose values for the separation constants, and, equally important, to decide how many terms to keep in the series. Similarly, the normalized scale length is non-intuitively connected to the real geometry. Lastly, while Whittaker functions are available in standard mathematical programs, their evaluation can be time-consuming.

Even with these negative factors, authors have been successful in making the solutions work to a certain extent. However, because of their complexity, they have not been as widely used as might be expected.

2.3. A new approach

After analysing the difficulties with the Whittaker function procedure, we have developed a new approach to solving the linear partial differential equation. The approach is also based on separation of variables, but is much more mathematically intuitive, thereby eliminating many of the difficulties just discussed. Specifically, the number of terms in the series, the choice of separation constants and the behaviour of the solutions are much easier to understand and obtain. This progress is the result of exploiting an underutilized mathematical insight into the nature of the separated solutions.

Several steps are required to reduce the problem to the desired form. First, we introduce a normalized flux and simple normalized coordinates into (2.4). The normalizations are defined by

(2.8)\begin{equation}\left. {\begin{array}{ll@{}} {\varPsi = {\varPsi_0}\psi }&{0 \le \psi \le 1,}\\ {R = {R_0}{r^{1/2}} = {R_0}{{(1 + {\varepsilon^2} + 2\varepsilon x)}^{1/2}}}&{ - 1 \le x \le 1,}\\ {Z = ay}&{ - \kappa \le y \le \kappa .} \end{array}} \right\}\end{equation}

Here, a is the plasma minor radius, $\varepsilon = a/{R_0}$ is the inverse aspect ratio and $\kappa $ is the plasma elongation, all assumed to be known inputs. The Grad–Shafranov equation and boundary condition reduce to

(2.9)\begin{equation}\left. \begin{gathered} (1 + \hat{\varepsilon }x)\dfrac{{{\partial^2}\psi }}{{\partial {x^2}}} + \dfrac{1}{{1 + {\varepsilon^2}}}\dfrac{{{\partial^2}\psi }}{{\partial {y^2}}} ={-} {\alpha^2}(1 + \hat{\varepsilon }\nu x)\psi \quad \psi (\textrm{surf)} = 0,\\ \hat{\varepsilon } = \dfrac{{2\varepsilon }}{{1 + {\varepsilon^2}}} < 1,\\ {{\alpha^2} = \dfrac{{2R_0^2{a^2}}}{{\varPsi_0^2}}\left( {{\mu_0}{p_0} + \dfrac{{{B_0}\delta B}}{{1 + {\varepsilon^2}}}} \right) \ge 0,}\\ \nu = \dfrac{{{\mu_0}{p_0}}}{{{\mu_0}{p_0} + {B_0}\delta B/(1 + {\varepsilon^2})}} \ge 0. \end{gathered} \right\}\end{equation}

There are several points worth discussing. First, note that $\psi = 0$ satisfies the equation and boundary condition. Our problem is thus an eigenvalue problem as has previously been pointed out by several authors (Goedbloed 1984; Takeda & Tokuda 1991; LoDestro & Pearlstein 1994; Pataki et al. 2013).

Non-trivial solutions exist only for special choices of the eigenvalue ${\alpha ^2}$. Only positive eigenvalues ${\alpha ^2} > 0$ need be considered, to ensure that the physical requirement of non-reversed toroidal current density is satisfied. Since $\alpha \propto 1/{\varPsi _0}$, the eigenvalue determines the unknown value of the flux on axis.

Second, for a given geometry, the parameter $\nu$ is assumed to be a known quantity since ${p_0}$ and ${B_0}\delta B$ are, at this point, known inputs. It is critical to recognize that values of both ${p_0}$ and ${B_0}\delta B$ are required for post-processing in order to evaluate various plasma quantities of interest. However, only the value of a single parameter $\nu$ is needed to determine the solution for $\psi $ itself. Furthermore, it is shown in § 6 that $\nu \approx {\beta _P}$, the poloidal beta, thus making it easy to choose reasonable values. To be more specific, the interesting range of $\nu$ is defined by

(2.10)\begin{equation}\left. {\begin{array}{ll@{}} {0 \le \nu \le {\nu_{\max }},}&{}\\ {\nu = 0}&{\textrm{Force}\;\textrm{free}\;\textrm{(}{p_0} = 0),}\\ {\nu = 1}&{{\beta_p} = 1\;\textrm{(}\delta B = 0),}\\ {\nu = {\nu_{\max }} = 1/\hat{\varepsilon }}&{\textrm{Inboard}\;\textrm{current}\;\textrm{reversal}\;\textrm{at}\;R = {R_0} - a\textrm{ (high}\;\beta ).} \end{array}} \right\}\end{equation}

Note that solutions can be easily found for $\nu > {\nu _{\max }}$. However, these are considered to be physically uninteresting since they correspond to a reversal in direction of the toroidal current density near the inboard midplane.

Third, the normalization requirement $\psi (\textrm{axis}) = 1$ is simple to satisfy. Since (2.9) is linear in $\psi $, the condition just corresponds to a rescaling of the eigenfunction. Lastly, observe that in the large-aspect-ratio limit, $\varepsilon = \hat{\varepsilon } = 0$, the equation reduces to the basic Cartesian form of Poisson's equation. The separated solutions are then intuitively simple sines, cosines or hyperbolic sines, cosines. The toroidal modifications are proportional to $\varepsilon < 1$. For most practical configurations of interest $\varepsilon $ is not very large: $\varepsilon \sim 0.2 - 0.5$. This, we shall see shortly, suggests a rapidly converging power-series representation of the radial separated variable functions. It also provides a primary motivation for our choice of independent variable x.

2.4. Separation of variables and a mathematical insight

The next step in the analysis is to solve (2.9) using separation of variables:

(2.11)\begin{equation}\psi (x,y) = \sum {{X_n}(x){Y_n}(y)} .\end{equation}

For the moment we do not specify the number of terms to be included in the summation. The equations for ${X_n},\;{Y_n}$ are given by

(2.12)\begin{equation}\left. \begin{gathered} {\dfrac{{{\textrm{d}^2}{Y_n}}}{{\textrm{d}{y^2}}} + h_n^2{Y_n} = 0,}\\ {(1 + \hat{\varepsilon }x)\dfrac{{{\textrm{d}^2}{X_n}}}{{\textrm{d}{x^2}}} + (k_n^2 + {\lambda^2}x){X_n} = 0,}\\ {k_n^2 = {\alpha^2} - \dfrac{{h_n^2}}{{1 + {\varepsilon^2}}},}\\ {{\lambda^2} = {\alpha^2}\hat{\varepsilon }\nu .} \end{gathered} \right\}\end{equation}

In these equations, $h_n^2$ is the separation constant. Also, for physical systems recall that ${\alpha ^2} > 0$ and correspondingly ${\lambda ^2} > 0$.

At this point we do not know how to choose the separation constants $h_n^2$, which may be positive or negative with magnitudes small, medium or large. It is here that mathematical insight proves crucial. The key point is the recognition that physically interesting solutions must have closed flux surfaces everywhere within the plasma. Therefore, if we make each term in our expansion separately produce closed flux surfaces about the geometric axis $x = 0,\;y = 0$, then when summing the terms, the total flux surfaces will automatically be closed around $x = 0,\;y = 0$, and almost certainly around the magnetic axis as well.

Mathematically, the surfaces for each expansion term must be shifted ellipses with the centre located approximately on the magnetic axis. Since the solution for each expansion term behaves like ${\psi _n} = {c_0} + {c_1}x + {c_2}y + {c_3}{x^2} + {c_4}{y^2} + \cdots $ near the geometric axis, the implication is that ${c_3},\;{c_4}$ must have the same sign, either positive or negative, for closed surfaces. Now, if $h_n^2 < 0$, then $k_n^2 > 0$ and the signs for ${c_3},\;{c_4}$ are opposite and the surfaces are open. If, however, we choose $h_n^2 > 0$ but limit its magnitude so that $k_n^2 > 0$, then each and every term in the summation produces closed flux surfaces near the axis. This insight translates into the following constraints on the separation constants:

(2.13)\begin{equation}\left. \begin{gathered} {0 \le h_n^2 \le h_{\max }^2,\quad {h_{\max }} = {{(1 + {\varepsilon^2})}^{1/2}}\alpha ,}\\ {k_{\max }^2 \ge k_n^2 \ge 0,\quad {k_{\max }} = \alpha .} \end{gathered} \right\}\end{equation}

Choosing the separation constants substantially outside this range leads to separate terms in the expansion having open flux surfaces near the axis, which then have to be compensated by other terms with comparable and cancelling open flux surfaces. This is often a delicate balance, sometimes requiring many terms in the solution or even possibly not having a solution at all.

Implementation of these constraints is the reason why realistic solutions require only a few terms in the summation. Specifically, in the sections that follow we show that realistic solutions require only 4 separation constants and either 7 or 12 terms in the summation depending on whether the configuration of interest is up–down symmetric or asymmetric. These numerical values are not empirical but are based on firm mathematical constraints.

Having laid the groundwork, we close this section by deriving the actual solutions for the separation functions ${X_n}(x)$ and ${Y_n}(y)$. The ${Y_n}(y)$ solutions are straightforward:

(2.14)\begin{equation}{Y_n}(y) = \left\{\begin{array}{@{}l} \cos ({h_n}y)\\ \sin ({h_n}y). \end{array} \right.\end{equation}

The ${X_n}(x)$ solutions are more complicated. We could write them in terms of Whittaker functions, but that would negate the usefulness of many of the steps so far introduced. Instead, let us take a step back and re-examine the ${X_n}(x)$ equation, given by (2.12). This equation has a regular singular point at $x ={-} 1/\hat{\varepsilon }\;(R = 0)$ and an irregular singular point at $x = R = \infty $. However, the regime of interest $- 1 \le x \le 1$ does not include either singularity. From a mathematical point of view, the equation is ‘boring’. This behaviour coupled with the fact that ${X_n}(x)$ reduces to sines and cosines as $\varepsilon \to 0$ suggests that simple intuitive representations of its two solutions can be written as

(2.15)\begin{equation}{X_n}(x) = \left\{\begin{array}{@{}l} {C_n}(x) = \sum {[{{\hat{a}}_m}{x^m}\,\textrm{cos}({k_n}x) + {{\hat{b}}_m}{x^m}\,\textrm{sin}({k_n}x)]} \\ {S_n}(x) = \sum {[{{\tilde{a}}_m}{x^m}\,\textrm{cos}({k_n}x) + {{\tilde{b}}_m}{x^m}\,\textrm{sin}({k_n}x)]} \,. \end{array} \right.\end{equation}

We expect each power series to be rapidly converging since the mth coefficients are proportional to ${\hat{\varepsilon }^m}$. It is a simple matter to calculate the recursion relations for the coefficients ${\hat{a}_m},\;{\hat{b}_m},\;{\tilde{a}_m},\;{\tilde{b}_m}$ as well as analytic first and second derivatives. The details are described in Appendix A. The end results are indeed highly accurate and computationally fast representations of ${C_n}(x)$ and ${S_n}(x)$ plus first and second derivatives. Stated differently, ${C_n}(x)$ and ${S_n}(x)$ can be viewed as known mathematical functions. They are Whittaker functions that just happen to be expressed in a particularly convenient way that exploits the ‘boring’ behaviour of the solutions in the regime of interest.

The tokamak equilibrium problem has now been formulated. As such, we are now in a position to investigate various applications of interest.

3. Smooth, up–down symmetric limiter equilibria

3.1. Number of terms in the series

The number of terms and corresponding separation constants in the $\psi (x,y)$ series is determined by the number of constraints that must be applied. The number of constraints is based on the following mathematical intuition. We start by recognizing that the flux surfaces near the axis are, by construction, closed surfaces. It then makes sense to require the flux, its slope and its curvature to match a specified desired model surface at four critical points around the circumference: (i) the inner midplane point, (ii) the outer midplane point, (iii) the upper maximum point and (iv) the lower maximum point. This, in principle, corresponds to 12 constraints.

For an up–down symmetric configuration, the lower maximum point constraint is automatically satisfied by choosing only the $\cos ({h_n}y)$ terms in the ${Y_n}(y)$ series. This reduces the number of constraints from 12 to 9. Furthermore, since $\cos ({h_n}y)$ is up–down symmetric, the slope constraints at the inner and outer midplane points are also automatically satisfied. This reduces the number of constraints from nine to seven.

The conjecture is that seven carefully chosen terms in the series for $\psi (x,y)$ should produce a set of closed flux surfaces whose boundary closely approximates the specified model surface.

3.2. Choice of the separation constants

How many separation constants are needed and what values should they have in order to be able to satisfy the seven constraints? Now, each separation constant leads to two terms in the series, one proportional to ${C_n}(x)$, the other to ${S_n}(x)$. Allowing three separation constants corresponds to six terms in the series (not enough by one) while four separation constants produce eight terms (one too many). The minimum number of required separation constants is, therefore, equal to four although we need to decide the fate of the one extra term.

One possibility is to add another constraint. A simpler idea, and the one utilized in our studies, is to choose one of the separation constants, say $h_1^2$, to make $k_1^2 = 0$:

(3.1)\begin{equation}{h_1} = {h_{\textrm{max}}} = {(1 + {\varepsilon ^2})^{1/2}}\alpha \to {k_1} = 0.\end{equation}

This is useful, because as shown in Appendix A, the function ${S_1}(x)$ with ${k_1} = 0$ is identically zero, thereby eliminating one term in the series, reducing the total number from eight to the desired seven. The conclusion is that for a smooth up–down symmetric configuration, the solution for $\psi (x,y)$ and its first and second derivatives can be written as

(3.2)\begin{equation}\left. {\begin{array}{l@{}} {\psi (x,y) = {c_1}\,\textrm{cos}({h_1}y)\,{C_1}(x) + \sum\limits_2^4 {\textrm{cos}({h_n}y)[{c_n}{C_n}(x) + {s_n}{S_n}(x)]} ,}\\ {\dfrac{{\partial \psi (x,y)}}{{\partial x}} = {c_1}\,\textrm{cos}({h_1}y)\,{{C^{\prime}}_1}(x) + \sum\limits_2^4 {\textrm{cos}({h_n}y)[{c_n}{{C^{\prime}}_n}(x) + {s_n}{{S^{\prime}}_n}(x)]} ,}\\ {\dfrac{{\partial \psi (x,y)}}{{\partial y}} ={-} {c_1}{h_1}\,\textrm{sin}({h_1}y)\,{C_1}(x) - \sum\limits_2^4 {{h_n}\,\textrm{sin}({h_n}y)[{c_n}{C_n}(x) + {s_n}{S_n}(x)]} ,}\\ {\dfrac{{{\partial^2}\psi (x,y)}}{{\partial {x^2}}} = {c_1}\,\textrm{cos}({h_1}y)\,{{C^{\prime\prime}}_1}(x) + \sum\limits_2^4 {\textrm{cos}({h_n}y)[{c_n}{{C^{\prime\prime}}_n}(x) + {s_n}{{S^{\prime\prime}}_n}(x)]} ,}\\ {\dfrac{{{\partial^2}\psi (x,y)}}{{\partial {y^2}}} ={-} {c_1}h_1^2\,\textrm{cos}({h_1}y)\,{C_1}(x) - \sum\limits_2^4 {h_n^2\,\textrm{cos}({h_n}y)[{c_n}{C_n}(x) + {s_n}{S_n}(x)]} .} \end{array}} \right\}\end{equation}

The seven unknown expansion coefficients to be determined by the constraints are $\boldsymbol{u} = [{c_1},{c_2},{s_2},{c_3},{s_3},{c_4},{s_4}]$.

The next step is the choice of the other separation constants ${h_2},\;{h_3},\;{h_4}$. Intuitively, we might expect that equal separation of their values would be useful so as to prevent two terms from being too similar, thereby reducing the amount of independent information in the series. But, this turns out to not exactly be the case.

We have, in fact, tried a wide range of ${h_n}$, all satisfying the constraint given by (2.13). A large fraction, but not all, of these choices yield very similar flux surfaces. Keep in mind that they are all analytic solutions to the Grad–Shafranov equation, the only difference being how closely they approximate the model surface. However, since the model surface is only a model surface, there is nothing fundamental about it. Each viable choice of ${h_n}$ produces very reasonable looking surfaces, all with the same elongation and triangularity, and similar values for beta and safety factor. In short, there is not a strong motivation to devote computational effort choosing ${h_n}$ to get as close a match to the model surface as possible.

Instead, our strategy has been to develop a simple, empirical, essentially universal algorithm for choosing ${h_n}$ that works robustly over a very wide range of magnetic geometries and plasma parameters. Very many test cases have shown that the choices below accomplish this goal:

(3.3)\begin{equation} \left.{\begin{array}{ll@{}} {{h_1} = {{(1 + {\varepsilon^2})}^{1/2}}\alpha }&{{k_1} = 0,}\\ {{h_2} = {{\left( {\dfrac{{35}}{{36}}} \right)}^{1/2}}{{(1 + {\varepsilon^2})}^{1/2}}\alpha }&{{k_2} = \left( {\dfrac{1}{6}} \right)\alpha ,}\\ {{h_3} = {{\left( {\dfrac{{13}}{{49}}} \right)}^{1/2}}{{(1 + {\varepsilon^2})}^{1/2}}\alpha }&{{k_3} = \left( {\dfrac{6}{7}} \right)\alpha ,}\\ {{h_4} = 0}&{{k_4} = \alpha .} \end{array}} \right\}\end{equation}

Observe that the optimal choices in ${k_n}$ are bunched near the maximum and minimum values and not equally spaced as might have been expected.

The number of terms in the series plus the separation constants have now been defined.

3.3. Model surface and matching constraints

A good choice for a model surface that produces smooth, up–down symmetric configurations is the well-known Miller profile (Miller et al. 1998). Although known as the ‘Miller profile’, Miller et al. actually credit the profile to an earlier paper (Todd et al. 1979). For the actual surface coordinates ${R_S},\;{Z_S}$, plus the corresponding normalized surface coordinates ${x_S},\;{y_S}$, the Miller profile is defined in terms of an angle-like parameter $\theta $, $0 \le \theta \le 2{\rm \pi} $, as follows:

(3.4)\begin{equation}\left. \begin{gathered} {{R_S} = {R_0} + a\,\textrm{cos}(\theta + \hat{\delta }\,\textrm{sin}\,\theta ) \to {x_S}(\theta ) = \textrm{cos}(\theta + \hat{\delta }\,\textrm{sin}\,\theta ) - \dfrac{\varepsilon }{2}\,\textrm{si}{\textrm{n}^2}(\theta + \hat{\delta }\,\textrm{sin}\,\theta ),}\\ {{Z_S} = a\kappa \,\textrm{sin}\,\theta \to {y_S}(\theta ) = \kappa \,\textrm{sin}\,\theta .} \end{gathered} \right\}\end{equation}

Here, $\kappa $ is the elongation, $\delta $ is the triangularity and $\hat{\delta } = \textrm{si}{\textrm{n}^{ - 1}}\delta $.

The next task is to define the seven constraints at the matching points: inner midplane point $\theta = {\rm \pi},\;{x_S} ={-} 1,\;{y_S} = 0$; outer midplane point $\theta = 0,\;{x_S} = 1,\;{y_S} = 0$; and upper maximum point $\theta = {\rm \pi}/2$, ${x_S} ={-} \delta - (\varepsilon /2)(1 - {\delta ^2})$, ${y_S} = \kappa $. The flux and slope constraints are easy to specify:

(3.5)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{a})\;\psi ( - 1,0) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\theta = {\rm \pi} ,}\\ {(\textrm{b})\;\psi (1,0) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\theta = 0,}\\ {(\textrm{c})\;\psi ( - {x_\delta },\kappa ) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\theta = {\rm \pi} /2,}\\ {(\textrm{d})\;{\psi_x}( - {x_\delta },\kappa ) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{slope}\;\textrm{at}\;\theta = {\rm \pi}/2,} \end{array}} \right\}\end{equation}

where we have defined ${x_\delta } = \delta + (\varepsilon /2)(1 - {\delta ^2})$.

The curvature constraints require a short analysis. Their evaluation requires the local expansion of the surface about the three critical matching points. A simple calculation yields

(3.6)\begin{equation}\left. {\begin{array}{ll@{}} {{x_S} + 1 \approx \dfrac{{(1 - \varepsilon ){{(1 - \hat{\delta })}^2}}}{{2{\kappa^2}}}y_S^2}&{\theta = {\rm \pi},}\\ {{x_S} - 1 \approx{-} \dfrac{{(1 + \varepsilon ){{(1 + \hat{\delta })}^2}}}{{2{\kappa^2}}}y_S^2}&{\theta = 0,}\\ {{y_S} - \kappa \approx{-} \dfrac{\kappa }{{2{{(1 - \varepsilon \delta )}^2}(1 - {\delta^2})}}{{({x_S} + {x_\delta })}^2}}&{\theta = {\rm \pi}/2.} \end{array}} \right\}\end{equation}

These relations can be converted into flux function constraints by expanding $\psi (x,y)$ about the same three points:

(3.7)\begin{equation}\left. {\begin{array}{ll@{}} {\psi ({x_S},{y_S}) \approx \dfrac{{\partial \psi ( - 1,0)}}{{\partial x}}({x_S} + 1) + \dfrac{{{\partial^2}\psi ( - 1,0)}}{{\partial {y^2}}}\dfrac{{y_S^2}}{2} = 0}&{\theta = {\rm \pi},}\\ {\psi ({x_S},{y_S}) \approx \dfrac{{\partial \psi (1,0)}}{{\partial x}}({x_S} - 1) + \dfrac{{{\partial^2}\psi (1,0)}}{{\partial {y^2}}}\dfrac{{y_S^2}}{2} = 0}&{\theta = 0,}\\ {\psi ({x_S},{y_S}) \approx \dfrac{{{\partial^2}\psi ( - {x_\delta },\kappa )}}{{\partial {x^2}}}\dfrac{{{{({x_S} + {x_\delta })}^2}}}{2} + \dfrac{{\partial \psi ( - {x_\delta },\kappa )}}{{\partial y}}({y_S} - \kappa ) = 0}&{\theta = {\rm \pi}/2.} \end{array}} \right\}\end{equation}

Combining (3.6) and (3.7) we see that the curvature constraints reduce to

(3.8)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{e})\;\dfrac{{{\psi_{yy}}( - 1,0)}}{{{\psi_x}( - 1,0)}} ={-} \dfrac{{(1 - \varepsilon ){{(1 - \hat{\delta })}^2}}}{{{\kappa^2}}} \equiv{-} {\varLambda_1}}&{\theta = {\rm \pi} ,}\\ {(\textrm{f})\;\dfrac{{{\psi_{yy}}(1,0)}}{{{\psi_x}(1,0)}} = \dfrac{{(1 + \varepsilon ){{(1 + \hat{\delta })}^2}}}{{{\kappa^2}}} \equiv {\varLambda_2}}&{\theta = 0,}\\ {(\textrm{g})\;\dfrac{{{\psi_{xx}}( - {x_\delta },\kappa )}}{{{\psi_y}( - {x_\delta },\kappa )}} = \dfrac{\kappa }{{{{(1 - \varepsilon \delta )}^2}(1 - {\delta^2})}} \equiv {\varLambda_3}}&{\theta = {\rm \pi}/2.} \end{array}} \right\}\end{equation}

Equations (3.5) and (3.8) are the seven constraints on the flux function that determine ${\alpha ^2}$ and the unknown expansion coefficients. In terms of obtaining realistic analytic equilibria, it is worth emphasizing that the only purpose of the model surface is to provide reasonable values for the surface curvature at the three critical matching points.

3.4. Solution

Using the results obtained in Appendix A, the constraint relations lead to a set of linear homogeneous algebraic equations for the expansion coefficients:

(3.9)\begin{equation}\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over A} }(\alpha ) \boldsymbol{\cdot} \boldsymbol{u} = 0.\end{equation}

Here, $\alpha $ appears as an eigenvalue but in a complex way, both as a ‘standard’ eigenvalue and also in the separation constants.

The following procedure can be used to determine $\alpha $ and the expansion coefficients. Assume that values for $\varepsilon ,\;\kappa ,\;\delta ,\;\nu $ have been specified. Set the coefficient ${c_1} = 1$. Now, solve the first six constraint equations, (a)–(f), for the six unknowns ${c_2},\;{c_3},\;{c_4},\;{s_2},\;{s_3},\;{s_4}$. This is a simple linear solve. The key point is that only if ${\alpha ^2}$ just happens to be the correct value will the seventh equation be satisfied. In general, it will not be satisfied and some form of iteration on $\alpha $ is required. Schematically, constraint (g) can be written as

(3.10)\begin{equation}{Q_1} + {Q_2}{c_2} + {Q_3}{s_2} + {Q_4}{c_3} + {Q_5}{s_3} + {Q_6}{c_4} + {Q_7}{s_4} = 0.\end{equation}

Just looking for a zero crossing of this function should work but in practice there are often spurious results. The reason is that sign changes occur when the function changes from plus to minus infinity at certain values of $\alpha $. These appear like spurious zero crossings when plotted versus $\alpha $. A better approach is to define an error function

(3.11)\begin{equation}E(\alpha ) = {\left( {\frac{{{Q_1} + {Q_2}{c_2} + {Q_3}{s_2} + {Q_4}{c_3} + {Q_5}{s_3} + {Q_6}{c_4} + {Q_7}{s_4}}}{{|{{Q_1}} |+ |{{Q_2}{c_2}} |+ |{{Q_3}{s_2}} |+ |{{Q_4}{c_3}} |+ |{{Q_5}{s_3}} |+ |{{Q_6}{c_4}} |+ |{{Q_7}{s_4}} |}}} \right)^2} = 0.\end{equation}

Clearly this function is always positive and bounded by unity. We carry out a scan in $\alpha $ and look for that value for which $E(\alpha )$ has a minimum corresponding to $E(\alpha ) = 0$. Our scanning procedure is straightforward, although somewhat ‘brute force’ in nature. It is a two-step process. First we perform a fine-scale scan over a hundred points to bracket the lowest $\alpha $ corresponding to minimum in $E(\alpha )$ with $|{E(\alpha )} |\le \textrm{Tol} \sim {10^{ - 3}}$. Second, we then use MATLAB's standard minimizing option to find the actual $\alpha $ to whatever accuracy is desired, typically six significant figures. The choice of a ‘brute force’ procedure is made to keep the search algorithm robust.

Even so, the numerical determination of $\alpha$ and the expansion coefficients is very fast with $\alpha $ typically calculated to six-figure accuracy (although this can easily be increased). This represents the numerical solution to the Grad–Shafranov equation. The CPU times, with no attempt at numerical optimization, depend only weakly on how many terms are maintained in the series for ${C_n}(x)$ and ${S_n}(x)$. For typical inverse aspect ratios $\varepsilon \sim 0.35$, approximately 10 terms are required for an accuracy of 5–6 significant figures in $\psi $, sufficient for practical applications. The CPU time is typically 0.5 s. For an accuracy of 15 significant figures, the number of terms increases to about 25 while the CPU time is about the same (a nice feature of MATLAB linear algebra calculations). Even in an extreme spherical tokamak case with $\varepsilon = 0.75$, 220 terms produce 15-figure accuracy with the CPU time again about the same. Overall, the ultimate accuracy of the solution is determined by the accuracy of the eigenvalue, not the accuracy of the expansion functions. Specifically, the numerical solutions satisfy the Grad–Shafranov equation to a comparable number of significant figures as maintained in the expansion functions, while the boundary condition accuracy is comparable to the number of significant figures in $\alpha $.

Since the CPU times are so short, we admittedly use overkill, and for simplicity always maintain 150 terms in the series. Overall, the total CPU time required to calculate an equilibrium is largely determined by the amount of post-processing desired. Typically, the evaluations of $\alpha $, the expansion coefficients, a flux surface plot and the key parameters of physical interest require about 2 s of CPU time. This corresponds to a single evaluation of $q(\psi )$ at the 95% surface. The post-processing can increase to about 1 min if a large number of ultrahigh-resolution $q(\psi )$ surfaces are required. Note also that to carry out due diligence with respect to verification of our numerical code, we have compared results against an established Grad–Shafranov solver FLOW that includes plasma flow (Guazzotto et al. 2004). In the limit of zero flow the agreement is very good, thereby verifying our numerical implementation.

3.5. Example: circle, ellipse, elongated D, inverse D

To demonstrate the solutions we calculate and plot the flux surfaces for five model plasma surfaces: (a) a circle, (b) an ellipse, (c) an elongated D, (d) a high-triangularity D and (e) an inverse D. Also illustrated are the midplane plots for the pressure and current density, normalized to their maximum values. The vertical axis is $Z/a$ while the horizontal axis is $R/{R_0}$. The circle and the ellipse are just reference cases to demonstrate the credibility plus versatility of the procedure. The elongated, high-triangularity and inverse D shapes are more relevant to modern tokamaks. All correspond to special choices of the Miller profile parameters. The actual post-processing physical parameters of interest are presented in § 6, after they have been defined. The flux surfaces and midplane plots are illustrated in figures 1–5 with the corresponding specific profile parameters listed in table 1.

Figure 1. Circle: $\varepsilon = 0.33,\;\kappa = 1,\;\delta = 0,\;\nu = 1$.

Figure 2. Ellipse: $\varepsilon = 0.25,\;\kappa = 2,\;\delta = 0,\;\nu = 1$.

Figure 3. Elongated D: $\varepsilon = 0.33,\;\kappa = 1.8,\;\delta = 0.4,\;\nu = 0.3$.

Figure 4. High-triangularity D: $\varepsilon = 0.4,\;\kappa = 2,\;\delta = 0.75,\;\nu = 1$.

Figure 5. Inverse D: $\varepsilon = 0.33,\;\kappa = 1.9,\;\delta ={-} 0.6,\;\nu = 0.5$.

Table 1. Input parameters for smooth, up–down symmetric limiter tokamaks.

There are two main points to observe. First, flux surfaces (c)–(e) have reasonable, realistic shapes, quite comparable to existing tokamak experiments, and similarly for the midplane profiles. Second, plotted as a red curve is the model Miller surface for each configuration. We see that the surfaces for the analytic equilibria closely match the Miller surface. The analysis for smooth limiter surfaces appears to work quite well over a wide range of parameters.

4. Up–down symmetric double-null divertor equilibria

The generalization of the analysis to include the treatment of double-null divertor systems is straightforward. Only two modifications are necessary. First, a new X-point model surface must be introduced, enabling us to obtain modified curvature conditions at the inner and outer midplane points. A new surface is needed because the Miller shape does not allow X-points. Using its midplane curvatures resulted in flux surfaces that were overly distorted at the desired X-points. Second, the curvature constraint at the maximum point must be replaced by the X-point constraint that ${B_R}$ vanishes at this point.

4.1. The double-null model surface

A simple form for the shape of a double-null surface consists of two intersecting ellipses. Note that when the current density vanishes on the plasma surface then not only must ${\psi _x}$ and ${\psi _y}$ vanish simultaneously at the X-point, but there is a further constraint that must also be imposed. Specifically, the ellipses must intersect perpendicular to each other. The double-null model surface satisfying these constraints can be expressed in terms of an angle-like parameter $\theta $ and is defined as follows.

• Inboard ellipse:

  (4.1)\begin{equation}{Z_S} = a{\kappa _0}\,\textrm{sin}\,\theta \quad {R_S} = {R_0} + a[{x_1} + (1 + {x_1})\,\textrm{cos}\,\theta ]\quad {\rm \pi}- {\theta _0} \le \theta \le {\rm \pi}+ {\theta _0}.\end{equation}

• Outboard ellipse:

  (4.2)\begin{equation}{Z_S} = a{\kappa _0}\,\textrm{sin}\,\theta \quad {R_S} = {R_0} + a[ - {x_2} + (1 + {x_2})\,\cos \,\theta ]\quad - {\theta _0} \le \theta \le {\theta _0}.\end{equation}
  Here,
  (4.3)\begin{equation}\left. \begin{gathered} {\dfrac{{{\kappa_0}}}{{{\kappa_X}}} = \dfrac{1}{{{{(1 - {\xi^2})}^{1/2}}}},}\\ {{x_1} = \dfrac{{\xi - {\delta_X}}}{{1 - \xi }},}\\ {{x_2} = \dfrac{{\xi + {\delta_X}}}{{1 - \xi }},}\\ {{\theta_0} = {{\tan }^{ - 1}}\left[ {\dfrac{{{{(1 - {\xi^2})}^{1/2}}}}{\xi }} \right],}\\ {\xi = \dfrac{{{{(1 - \delta_X^2)}^{1/2}}}}{{{\kappa_X} - {{(1 - \delta_X^2)}^{1/2}}}},} \end{gathered} \right\}\end{equation}

and ${\kappa _X},\;{\delta _X}$ at the X-point are known inputs replacing $\kappa ,\;\delta $ for the smooth surface.

Observe that for real solutions to exist we require $0 \le \xi \le 1$. The left limit is always satisfied but the right limit requires that

(4.4)\begin{equation}{\kappa _X} > 2{(1 - \delta _X^2)^{1/2}}.\end{equation}

This constraint arises because in low-elongation, circular-like plasmas it is difficult to shrink the X-point angle to ${\rm \pi} /2$. In practice if we assume ${\kappa _X} > 2$ the condition is automatically satisfied, and this is not a very strong practical limitation. A cleverer choice of model surface could lower this limit if important.

4.2. The ${B_R} = 0$ constraint

In a double-null system there is no longer a curvature constraint at the X-point as there is for a smooth surface. Instead, we replace the curvature constraint given by term (g) in (3.8) with

(4.5)\begin{equation}{\psi _y}( - {x_X},{\kappa _X}) = 0\quad {x_X} = {\delta _X} + (\varepsilon /2)(1 - \delta _X^2).\end{equation}

4.3. The full set of double-null constraints

The full set of double-null constraints is summarized below:

(4.6)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{a})\;\psi ( - 1,0) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\textrm{inner}\;\textrm{midplane,}}\\ {(\textrm{b})\;\psi (1,0) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\textrm{outer}\;\textrm{midplane,}}\\ {(\textrm{c})\;\psi ( - {x_X},{\kappa_X}) = 0}&{\textrm{Flux}\;\textrm{surface}\;\textrm{at}\;\textrm{X - point,}}\\ {(\textrm{d})\;{\psi_x}( - {x_X},{\kappa_X}) = 0}&{{B_Z} = 0\;\textrm{at}\;\textrm{X - point,}}\\ {(\textrm{e})\;\dfrac{{{\psi_{yy}}( - 1,0)}}{{{\psi_x}( - 1,0)}} ={-} {\varLambda_1}}&{\textrm{Curvature}\;\textrm{at}\;\textrm{inner}\;\textrm{midplane,}}\\ {(\textrm{f})\;\dfrac{{{\psi_{yy}}(1,0)}}{{{\psi_x}(1,0)}} = {\varLambda_2}}&{\textrm{Curvature}\;\textrm{at}\;\textrm{outer}\;\textrm{midplane,}}\\ {(\textrm{g})\;{\psi_y}( - {x_X},{\kappa_X}) = 0}&{{B_R} = 0\;\textrm{at}\;\textrm{X - point}\textrm{.}} \end{array}} \right\}.\end{equation}

In this equation, a short calculation shows that the values of ${\varLambda _1},\;{\varLambda _2}$ for the intersecting ellipse model are given by

(4.7)\begin{equation}\left. \begin{gathered} {{\varLambda_1} = \dfrac{{(1 - \varepsilon )(1 - {\delta_X})(1 + \xi )}}{{\kappa_X^2}},}\\ {{\varLambda_2} = \dfrac{{(1 + \varepsilon )(1 + {\delta_X})(1 + \xi )}}{{\kappa_X^2}}.} \end{gathered} \right\}\end{equation}

Our expansion for $\psi (x,y)$ maintains the same form as (3.2). The constraint relations determining the expansion coefficients are nearly identical to the smooth surface constraints except for the ${B_R} = 0$ replacement and the new values for ${\varLambda _1},\;{\varLambda _2}$. Thus, the same solution procedure with similar CPU timing can be used to calculate equilibria.

4.4. Examples: standard tokamak, spherical tokamak

Double-null divertor analytic equilibria are illustrated for two examples: (a) a tokamak with standard aspect ratio, low ${\beta _p}$ and (b) a spherical tokamak with high aspect ratio, high ${\beta _p}$. The flux surfaces and midplane plots are illustrated in figures 6 and 7. The input parameters are listed in table 2.

Figure 6. Standard double null: $\varepsilon = 0.33,\;{\kappa _X} = 2,\;{\delta _X} = 0.5,\;\nu = 0.4.$

Figure 7. Spherical tokamak: $\varepsilon = 0.75,\;{\kappa _X} = 2.4,\;{\delta _X} = 0.8,\;\nu = 1.04$.

Table 2. Input parameters for up–down symmetric, double-null divertor tokamaks.

We again see that the analytic equilibrium model produces realistic looking results, and has no difficulty reproducing exact X-points. Also, the fit of the analytic solutions closely matches the two-ellipse model surface.

5. Up–down asymmetric single-null divertor equilibria

The analysis thus far presented can be further generalized to describe an up–down asymmetric single-null divertor configuration. This generalization requires combining many of the results already derived. The one main new step is the addition of several additional terms (but not additional separation constants) in the expansion to account for the asymmetry.

5.1. The single-null model surface

To begin, we introduce a convenient shape for the single-null divertor surface that can again be expressed parametrically in terms of an angle-like variable $\theta $. This is a combination of our previous surfaces. The inputs to the shape are $\kappa ,\;\delta $ for the smooth upper half and ${\kappa _X},\;{\delta _X}$ for the single-null lower half. We also need to specify the inverse aspect ratio $\varepsilon $. The shape for the upper half $(0 \le \theta \le {\rm \pi})$ corresponds to the Miller profile

(5.1)\begin{equation}{Z_S} = a\kappa \,\textrm{sin}\,\theta \quad {R_S} = {R_0} + a\,\textrm{cos}(\theta + \hat{\delta }\,\textrm{sin}\,\theta )\quad 0 \le \theta \le {\rm \pi} .\end{equation}

The lower half is subdivided into the two intersecting ellipses model:

(5.2)\begin{equation}\left. \begin{gathered} {{Z_S} = a{\kappa_0}\,\textrm{sin}\,\theta \quad {R_S} = {R_0} + a[{x_1} + (1 + {x_1})\,\textrm{cos}\,\theta ]\quad {\rm \pi}\le \theta \le {\rm \pi}+ {\theta_0},}\\ {{Z_S} = a{\kappa_0}\,\textrm{sin}\,\theta \quad {R_S} = {R_0} + a[ - {x_2} + (1 + {x_2})\,\textrm{cos}\,\theta ]\quad - {\theta_0} \le \theta \le 0.} \end{gathered} \right\}\end{equation}

5.2. The constraints

The constraints for an asymmetric single-null divertor configuration again involve the flux, the slope and the curvature at four critical points: outer midplane, inner midplane, upper maximum point and lower X-point. There are, however, some differences from symmetric systems that are described below. We begin with the flux constraints on $\psi (x,y)$, which are straightforward:

(5.3)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{a})\;\psi ( - 1,0) = 0}&{\textrm{Inner}\;\textrm{midplane,}}\\ {(\textrm{b})\;\psi (1,0) = 0}&{\textrm{Outer}\;\textrm{midplane,}}\\ {(\textrm{c})\;\psi ( - {x_\delta },\kappa ) = 0}&{\textrm{Upper}\;\textrm{maximum}\;\textrm{point,}}\\ {(\textrm{d})\;\psi ( - {x_X}, - {\kappa_X}) = 0}&{\textrm{Lower}\;\textrm{X - point}\textrm{.}} \end{array}} \right\}\end{equation}

The slope constraints are modified because $\textrm{d}x/\textrm{d}y$ is no longer automatically zero at the midplane points since the configuration no longer has up–down symmetry. Our model surfaces for the upper and lower sections both require that $\textrm{d}x/\textrm{d}y = 0$ at the outer and inner midplane points, as can be seen from (5.1) and (5.2). Nevertheless, this requirement must be enforced on our analytic solution since, by construction, $\psi (x,y)$ is not automatically up–down symmetric. A second modification occurs at the lower X-point. Here, it is not actually a slope constraint that applies, but instead the requirement that both components of magnetic field simultaneously vanish. Thus, there are two first-derivative constraints at this point. The total first-derivative constraints can be written as

(5.4)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{e})\;\dfrac{{\partial \psi ( - 1,0)}}{{\partial y}} = 0}&{\textrm{Inner}\;\textrm{midplane,}}\\ {(\textrm{f})\;\dfrac{{\partial \psi (1,0)}}{{\partial y}} = 0}&{\textrm{Outer}\;\textrm{midplane,}}\\ {(\textrm{g})\;\dfrac{{\partial \psi ( - {x_\delta },\kappa )}}{{\partial x}} = 0}&{\textrm{Upper}\;\textrm{maximum}\;\textrm{point,}}\\ {(\textrm{h})\;\dfrac{{\partial \psi ( - {\delta_X}, - {\kappa_X})}}{{\partial x}} = 0}&{\textrm{Lower}\;\textrm{X - point,}}\\ {(\textrm{i})\;\dfrac{{\partial \psi ( - {\delta_X}, - {\kappa_X})}}{{\partial y}} = 0}&{\textrm{Lower}\;\textrm{X - point}.} \end{array}} \right\}\end{equation}

The curvature constraints at the midplane points are also modified. The reason is that the second derivatives ${\partial ^2}x/\partial {y^2}$ are different for the upper and lower surfaces. If we imagine a smooth, rather than abrupt transition from one surface to the other, then the curvature constraint reduces to the average of the two separate curvature constraints. Formally, we write $x = ({x_{\textrm{Upper}}} + {x_{\textrm{Lower}}})/2$ to obtain the midplane curvature constraints. Also, at the X-point a curvature constraint is unnecessary. Consequently, there are three curvature constraints that must be applied, which are given by

(5.5)\begin{equation} \left. {\begin{array}{ll@{}} {(\textrm{j})\;\dfrac{{{\psi _{yy}}( - 1,0)}}{{{\psi _x}( - 1,0)}} = - {\varLambda _1}}&{\textrm{Curvature}\;\textrm{at}\;\textrm{inner}\;\textrm{midplane}}\\ {(\textrm{k})\;\dfrac{{{\psi _{yy}}(1,0)}}{{{\psi _x}(1,0)}} = {\varLambda _2}}&{\textrm{Curvature}\;\textrm{at}\;\textrm{outer}\;\textrm{midplane}}\\ {(\textrm{l})\;\dfrac{{{\psi _{xx}}( - {x_\delta },\kappa )}}{{{\psi _y}( - {x_\delta },\kappa )}} = {\varLambda _3}}&{\textrm{Curvature}\;\textrm{at}\;\textrm{top}\;\textrm{maximum}} \end{array}} \right\}, \end{equation}

where

(5.6)\begin{equation}\left. \begin{gathered} {{\varLambda_1} = \dfrac{1}{2}\left[ {\dfrac{{(1 - \varepsilon ){{(1 - \hat{\delta })}^2}}}{{{\kappa^2}}} + \dfrac{{(1 - \varepsilon )(1 - {\delta_X})(1 + \xi )}}{{\kappa_X^2}}} \right],}\\ {{\varLambda_2} = \dfrac{1}{2}\left[ {\dfrac{{(1 + \varepsilon ){{(1 + \hat{\delta })}^2}}}{{{\kappa^2}}} + \dfrac{{(1 + \varepsilon )(1 + {\delta_X})(1 + \xi )}}{{\kappa_X^2}}} \right],}\\ {{\varLambda_3} = \dfrac{\kappa }{{{{(1 - \varepsilon \delta )}^2}(1 - {\delta^2})}}.} \end{gathered} \right\}\end{equation}

All told, there are 12 constraints that must be satisfied by our analytic solution.

5.3. The single-null flux function expansion

As before, we write the analytic solution as a sum of terms containing four separation constants. The most general solution satisfying the four-separation-constant constraint can be written as

(5.7)\begin{equation}\begin{aligned} \psi (x,y) & =\sum\limits_1^4 {\textrm{cos}({h_n}y)[{c_n}{C_n}(x) + {s_n}{S_n}(x)]} \\ & \quad + \,\sum\limits_1^4 {\textrm{sin}({h_n}y)[{c_{n + 4}}{C_n}(x) + {s_{n + 4}}{S_n}(x)].} \end{aligned}\end{equation}

Observe that the expansion in principle contains 16 unknown coefficients whereas our formulation has only 12 constraints. This problem can be resolved by choosing the separation constants in the same way as for the up–down symmetric systems. Specifically, one separation constant is chosen so that ${k_1} = 0$ while a second is chosen so that ${h_4} = 0$. With these two choices plus the same intermediate separation constants, the total number of unknown coefficients is reduced to 12. The counting of constraints and unknown coefficients is now consistent. The modified expansion for the flux can now be written as

(5.8)\begin{equation}\begin{aligned} \psi (x,y) & ={c_1}\,\textrm{cos}({h_1}y)\,{C_1}(x) + \sum\limits_2^4 {\textrm{cos}({h_n}y)[{c_n}{C_n}(x) + {s_n}{S_n}(x)]} \\ & \quad+ \,{c_5}\,\textrm{sin}({h_1}y)\,{C_5}(x) + \sum\limits_2^3 {\textrm{sin}({h_n}y)[{c_{n + 4}}{C_n}(x) + {s_{n + 4}}{S_n}(x)].} \end{aligned}\end{equation}

5.4. Solution procedure

The solution procedure is the same as for the symmetric case except that we now have a 12 × 12 set of linear homogeneous algebraic equations with $\alpha $ as the eigenvalue:

(5.9)\begin{equation}\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over A} }(\alpha )\boldsymbol{\cdot }\boldsymbol{u} = 0.\end{equation}

Here, $\boldsymbol{u} = [{c_1},{c_2},{s_2},{c_3},{s_3},{c_4},{s_4},{c_5},{c_6},{s_6},{c_7},{s_7}]$. We again set ${c_1} = 1$, solve the first 11 constraint equations for the remaining ${c_j},{s_j}$ and then iterate the value of $\alpha $ until the 12th equation is satisfied. The equilibrium CPU time is only slightly changed, increasing from about 0.5 to about 1 s for an equilibrium calculation where 150 terms are maintained in the series for ${C_n},\;{S_n}$.

5.5. Examples: standard tokamak, high-${\beta _p}$ tokamak

Single-null divertor analytic equilibria are illustrated for two examples: (a) a tokamak with standard $\varepsilon $, standard ${\beta _p}$ and (b) a tokamak with high ${\beta _p}$, small $\varepsilon $. The flux surfaces and midplane plots are illustrated in figures 8 and 9. The input parameters are listed in table 3.

Figure 8. Standard single null: $\varepsilon = 0.33,\;\kappa = 1.6,\;\delta = 0.4,\;{\kappa _X} = 2,\;{\delta _X} = 0.5,\;\nu = 1$.

Figure 9. High-${\beta _P}$ single null: $\varepsilon = 0.25,\;\kappa = 1.8,\;\delta = 0.6$, ${\kappa _X} = 2.1,\;{\delta _X} = 0.8,\;\nu = 2.1$.

Table 3. Input parameters for up–down asymmetric, single-null divertor tokamaks.

Once again, we see that the analytic solutions have no difficulty generating realistic single-null divertor flux surfaces. The outer surfaces are quite close to the two-ellipse model surface and the midplane profiles smoothly approach zero on the surface.

6. Plasma parameters

The examples presented have shown that our analysis produces realistic flux surfaces over a wide range of interesting tokamak plasma configurations. In addition to flux surface plots, the results include analytic expressions for the flux $\psi (x,y)$ normalized so that $\psi (\textrm{axis}) = 1$, its first and second derivatives and a value for the eigenvalue $\alpha $. We can use these results to evaluate important plasma parameters of physical interest. To carry out this task several additional input quantities must be specified. A list of the original input parameters plus the new post-processing parameters is given below.

• Initial equilibrium parameters:

  (6.1)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{a})\;\varepsilon }&{\textrm{Inverse}\;\textrm{aspect}\;\textrm{ratio,}}\\ {(\textrm{b})\;\kappa ,\;{\kappa_X}}&{\textrm{Elongation,}}\\ {(\textrm{c})\;\delta ,\;{\delta_X}}&{\textrm{Triangularity,}}\\ {(\textrm{d})\;\nu }&{\textrm{Parameter}\;\textrm{approximately}\;\textrm{equal}\;\textrm{to}\;{\beta_P}.} \end{array}} \right\}\end{equation}

• Post-processing parameters:

  (6.2)\begin{equation}\left. {\begin{array}{ll@{}} {(\textrm{e})\;{R_0}}&{\textrm{Major}\;\textrm{radius,}}\\ {(\textrm{f})\;{B_0}}&{\textrm{Vacuum}\;\textrm{toroidal}\;\textrm{field}\;\textrm{at}\;R = {R_0},}\\ {(\textrm{g})\;{q_0}}&{\textrm{Safety}\;\textrm{factor}\;\textrm{on}\;\textrm{the}\;\textrm{magnetic}\;\textrm{axis}\textrm{.}} \end{array}} \right\}\end{equation}

With these as inputs, we can evaluate the plasma parameters of interest. It is convenient for mathematical compactness to express the results in terms of ${\beta _0}$, the toroidal beta on axis (i.e. ${\beta _0} = 2{\mu _0}{p_0}/B_0^2$) rather than ${q_0}$. However, there is a one-to-one relation between ${q_0}$ and ${\beta _0}$ that can easily be derived without the need for any iteration on the equilibrium parameters. The explicit relation between ${q_0}$ and ${\beta _0}$ is given at the end of this subsection. The desired plasma parameters are as follows.

• Plasma pressure on axis ${p_0}$:

  (6.3)\begin{equation}{p_0} = {\beta _0}\left( {\frac{{B_0^2}}{{2{\mu_0}}}} \right).\end{equation}

• Plasma diamagnetism $\delta B/{B_0}$:

  (6.4)\begin{equation}\frac{{\delta B}}{{{B_0}}} = \frac{{{\beta _0}(1 + {\varepsilon ^2})}}{2}\left( {\frac{{1 - \nu }}{\nu }} \right).\end{equation}

• Magnetic flux on axis ${\varPsi _0}$:

  (6.5)\begin{equation}{\varPsi _0} = \frac{{\varepsilon {B_0}R_0^2}}{\alpha }{\left( {\frac{{{\beta_0}}}{\nu }} \right)^{1/2}}.\end{equation}

• Volume -averaged toroidal beta ${\beta _T}$:

  (6.6)\begin{equation}{\beta _T} = \frac{{2{\mu _0}}}{{B_0^2}}{\langle p\rangle _V} = {\beta _0}\frac{{\int {{\psi ^2}\,\textrm{d}x\,\textrm{d}y} }}{{\int {\textrm{d}x\,\textrm{d}y} }}.\end{equation}

• Volume -averaged poloidal beta ${\beta _P}$:

  (6.7)\begin{equation}{\beta _P} = \frac{{2{\mu _0}{{\langle p\rangle }_V}}}{{{{\langle B_P^2\rangle }_V}}} = \frac{{2{\mu _0}{{\langle p\rangle }_V}}}{{{{\langle 2{\mu _0}p + ({F^2}/R_0^2 - B_0^2)/r\rangle }_V}}} = \frac{{\nu \int {{\psi ^2}\,\textrm{d}x\,\textrm{d}y} }}{{\int {\left( {\frac{{1 + \nu \hat{\varepsilon }x}}{{1 + \hat{\varepsilon }x}}} \right){\psi ^2}\,\textrm{d}x\,\textrm{d}y} }}.\end{equation}

• Toroidal plasma current $I$:

  (6.8)\begin{equation}{\mu _0}I = {\mu _0}\int_A {{J_\phi }\,\textrm{d}A} = \varepsilon {B_0}{R_0}\alpha {\left( {\frac{{{\beta_0}}}{\nu }} \right)^{1/2}}\int {\left( {\frac{{1 + \nu \hat{\varepsilon }x}}{{1 + \hat{\varepsilon }x}}} \right)\psi \,\textrm{d}x\,\textrm{d}y} .\end{equation}

• Normalized internal inductance per unit length ${l_i}$:

  (6.9)\begin{equation}{l_i} = \frac{{2{L_i}}}{{{\mu _0}{R_0}}} = \frac{2}{{\mu _0^2{I^2}{R_0}}}\int {B_P^2\,\textrm{d}{\kern 1pt} r} = \frac{{4{\rm \pi} }}{{{\alpha ^2}}}\frac{{\int {\left( {\frac{{1 + \nu \hat{\varepsilon }x}}{{1 + \hat{\varepsilon }x}}} \right){\psi ^2}\,\textrm{d}x\,\textrm{d}y} }}{{{{\left[ {\int {\left( {\frac{{1 + \nu \hat{\varepsilon }x}}{{1 + \hat{\varepsilon }x}}} \right)\psi \,\textrm{d}x\,\textrm{d}y} } \right]}^2}}}.\end{equation}

• Kink safety factor ${q_\ast }$ ($\kappa = {\kappa _{95}}$ for divertor surfaces):

  (6.10)\begin{equation}{q_\ast } = \frac{{2{\rm \pi} {a^2}{B_0}}}{{{\mu _0}{R_0}I}}\left( {\frac{{1 + {\kappa^2}}}{2}} \right) = \frac{{{\rm \pi} \varepsilon }}{\alpha }{\left( {\frac{\nu }{{{\beta_0}}}} \right)^{1/2}}\frac{{1 + {\kappa ^2}}}{{\int {\left( {\frac{{1 + \nu \hat{\varepsilon }x}}{{1 + \hat{\varepsilon }x}}} \right)\psi \,\textrm{d}x\,\textrm{d}y} }}.\end{equation}

• Local safety factor $q(\psi )$:

  (6.11)\begin{equation}q(\psi ) = \frac{{F(\psi )}}{{2{\rm \pi} }}\int {\frac{{\textrm{d}l}}{{{R^2}{B_P}}}} = \frac{{F(\psi )}}{{{R_0}{B_0}}}{\left( {\frac{\nu }{{{\beta_0}}}} \right)^{1/2}}\frac{{\varepsilon \alpha }}{{2{\rm \pi} }}{\oint {\left( {\frac{{x_\theta^2 + \rho y_\theta^2}}{{\psi_y^2 + \rho \psi_x^2}}} \right)} ^{1/2}}\frac{{\textrm{d}\theta }}{r}.\end{equation}

• Safety factor on axis ${q_0}$ and inverse relation ${\beta _0} = {\beta _0}({q_0})$:

  (6.12)\begin{equation}\left. \begin{gathered} {{q_0} = {{\left[ {\dfrac{{{B_\phi }}}{{{{({\varPsi_{RR}}{\varPsi_{ZZ}})}^{1/2}}}}} \right]}_{\varPsi = {\varPsi_0}}} = {{\left[ {\dfrac{{F(\psi )}}{{{R_0}{B_0}}}{{\left( {\dfrac{\nu }{{{\beta_0}}}} \right)}^{1/2}}\dfrac{{\varepsilon \alpha }}{{r{{({\psi_{xx}}{\psi_{yy}})}^{1/2}}}}} \right]}_{\psi = 1}}}\\ {\textrm{or}\quad {\beta_0} = \dfrac{{\nu {\varepsilon^2}{\alpha^2}}}{{q_0^2{{({r^2}{\psi_{xx}}{\psi_{yy}})}_{\psi = 1}} - (1 + {\varepsilon^2})(1 - \nu ){\varepsilon^2}{\alpha^2}}}.} \end{gathered} \right\}\textrm{ }\end{equation}

In these expressions recall that $r(x) = 1 + {\varepsilon ^2} + 2\varepsilon x$. Also, note from (6.12) that the ratio ${\beta _0}/\nu $ is finite as $\nu \to 0$ so there is no difficulty taking this limit. The plasma parameters of interest have now been defined.

As examples we have evaluated the dimensionless parameters ${\beta _T}$, ${\beta _P}$, ${l_i}$, ${q_{95}}$ and ${q_\ast }$ for the test cases previously presented, assuming that ${q_0} = 1$. Also given are the corresponding eigenvalues $\alpha $ in truncated form. The actual $\alpha $ values are accurate to at least six significant figures. Lastly, keep in mind that ${q_{95}} = q(\psi = 0.05)$. The plasma parameters are listed in table 4.

Table 4. Dimensionless plasma parameters for the test examples.

We observe that the analytic equilibria have no difficulty modelling configurations covering a wide range of operating space for smooth, double-null and single-null plasmas. One interesting value worth mentioning is ${q_{95}} = 18.2$ for the high-${\beta _P}$ spherical tokamak. This high value occurs because the inboard poloidal magnetic field is very small because of the high ${\beta _P}$. The implication is that a separatrix is closely approaching the inboard midplane, corresponding to the equilibrium beta limit.

7. Conclusions

We have derived simple, general, realistic, robust, analytic tokamak equilibria which are solutions to the Grad–Shafranov equation. They apply to a wide range of configurations: smooth surface, double-null divertor, single-null divertor, arbitrary $\varepsilon ,\;\kappa ,\;\delta ,\;{\beta _T},\;{\beta _P}$. Simple expressions for the plasma parameters of interest have also been derived. The model provides analytic expressions for the flux and its first and second derivatives. A simple code has been developed that typically takes of the order of two seconds to plot the flux surfaces and evaluate the plasma parameters of interest. The code is available to all interested readers. The code can also switch from a lower-null to an upper-null divertor, or to a smooth surface with different upper and lower shapes. These are trivial changes and as such there has been no need to explicitly discuss them in the text.

This first paper assumes simple but realistic boundary constraints in which the pressure, pressure gradient and current density all vanish on the plasma surface. The follow-on paper generalizes the analysis to include pedestals in these same quantities. It also includes the highly localized fraction of the bootstrap current appearing at the plasma edge in the presence of a density pedestal. A further generalization allows toroidal plasma flow. Again, all of these new effects are included in simple, analytic solutions to the Grad–Shafranov equation (including flow modifications when appropriate).

Editor William Dorland thanks the referees for their advice in evaluating this article.