2.1. Starting equations and boundary conditions

The tokamak equilibria of interest are described by the familiar Grad–Shafranov equation (Lüst & Schlüter Reference Lüst and Schlüter1957; Grad & Rubin Reference Grad and Rubin1958; Shafranov Reference Shafranov1966):

(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 Reference Goedbloed1984; Takeda & Tokuda Reference Takeda and Tokuda1991; LoDestro & Pearlstein Reference LoDestro and Pearlstein1994; Pataki et al. Reference Pataki, Cerfon, Freidberg, Greengard and O'Neil2013).

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.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. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998). Although known as the ‘Miller profile’, Miller et al. actually credit the profile to an earlier paper (Todd et al. Reference Todd, Manickam, Okabayashi, Chance, Grimm, Greene and Johnson1979). 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. Reference Guazzotto, Betti, Manickam and Kaye2004). 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.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.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.