3.1. Sources of surface current

As stated in § 1, there are two sources of surface currents. First, even if the pressure gradient and toroidal current density vanish smoothly at the plasma edge, there can still be a finite edge pressure, creating a pressure pedestal. From MHD pressure balance across the surface, it follows that this pressure jump must be balanced by a perpendicular surface current; that is, the current direction is perpendicular to both the surface normal and the average magnetic field across the current sheet. Its magnitude is proportional to the difference in fields across the current sheet. Mathematically, the pressure pedestal can be modelled by leaving the free function ${F^2}(\varPsi /{\varPsi _0})$ unchanged but modifying $p(\varPsi /{\varPsi _0})$ in the plasma core as follows:

(3.1)\begin{equation}p(\varPsi /{\varPsi _0}) = {p_0}{\left( {\frac{\varPsi }{{{\varPsi _0}}}} \right)^2} \to p(\varPsi /{\varPsi _0}) = {p_0}\left[ {{f_P} + (1 - {f_P}){{\left( {\frac{\varPsi }{{{\varPsi _0}}}} \right)}^2}} \right].\end{equation}

Note that $p(\textrm{axis}) = p(1) = {p_0}$ and $p(\textrm{surf}) = p(0) = {f_P}{p_0}$. We see that ${f_P} = p(\textrm{surf})/p(\textrm{axis)}$ is the fractional height of the pressure pedestal, and is an additional input quantity. A pressure jump leaves the normalized Grad–Shafranov equation unchanged, but with modified definitions of $\alpha $ and $\nu $:

(3.2)\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 ,}\\ {{\alpha^2} = \dfrac{{2R_0^2{a^2}}}{{\varPsi_0^2}}\left[ {{\mu_0}{p_0}(1 - {f_P}) + \dfrac{{{B_0}\delta B}}{{1 + {\varepsilon^2}}}} \right],}\\ {\nu = \dfrac{{{\mu_0}{p_0}(1 - {f_P})}}{{{\mu_0}{p_0}(1 - {f_P}) + {B_0}\delta B/(1 + {\varepsilon^2})}}.} \end{gathered} \right\}\end{equation}

The second component of surface current flows perpendicular to the surface but parallel to the average field across the surface. It is assumed to be generated by the portion of the total bootstrap current localized near the edge. This can be seen from the large-aspect-ratio approximation for the flux surface averaged bootstrap current given by (Helander & Sigmar Reference Helander and Sigmar2002)

(3.3)\begin{equation}{J_B}(\psi ) \equiv \left\langle {\frac{{{\boldsymbol{J}_B}\boldsymbol{\cdot }\boldsymbol{B}}}{B}} \right\rangle \approx{-} 2.42{\varepsilon ^{1/2}}{R_0}p\left[ {\frac{1}{n}\frac{{\textrm{d}n}}{{\textrm{d}\psi }} + 0.054\frac{1}{T}\frac{{\textrm{d}T}}{{\textrm{d}\psi }}} \right].\end{equation}

Near the edge, the density and temperature pedestals produce a large localized contribution to the bootstrap current, which corresponds to the parallel component of the surface current.

A convenient definition of the bootstrap fraction is described below. This definition focuses on the ${\boldsymbol{e}_\phi }$ component of bootstrap current, which is simpler to measure experimentally as compared to the actual bootstrap current which flows parallel to the magnetic field. We denote the parallel surface current by ${\boldsymbol{K}_\parallel } = {K_{{\parallel} \phi }}{\boldsymbol{e}_\phi } + {K_{{\parallel} P}}{\boldsymbol{e}_P}$, where ${\boldsymbol{e}_P} = {\boldsymbol{B}_P}/{B_P}$ and ${K_{{\parallel} \phi }}$ is the component of interest. The corresponding toroidal current generated by ${K_{{\parallel} \phi }}$ is given by

(3.4)\begin{equation}{I_{{\parallel} \phi }} = \oint {{K_{{\parallel} \phi }}\,\textrm{d}{l_P},}\end{equation}

with $\textrm{d}{l_P}$ the differential poloidal arc length. Next, note that the total toroidal current flowing in the plasma is $\hat{I} = I + {I_{{\parallel} \phi }}$, where I is the contribution from the plasma core calculated in Part 1. The definition of the edge-localized portion of bootstrap current ${f_B}$ is taken as

(3.5)\begin{equation}{f_B} \equiv \frac{{{I_{{\parallel} \phi }}}}{{\hat{I}}} = \frac{{\hat{I} - I}}{{\hat{I}}} = 1 - \frac{I}{{\hat{I}}}.\end{equation}

The quantity ${f_B}$ is an additional input parameter.

3.2. Surface current analysis

The surface current analysis can now be carried out as follows. Assume we have obtained a smooth edge solution based on the analysis in Part 1. Specifically, we assume that just inside the surface current we know ${B_\phi }(\theta )$ and ${B_P}(\theta )$. Next, we give values for the two new surface current fractions ${f_P}$ and ${f_B}$. Using this information we calculate the fields just outside the surface current layer ${\hat{B}_\phi }(\theta )$ and ${\hat{B}_P}(\theta )$, leading to a self-consistent surface current constraint.

To begin, we make use of Ampere's law and pressure balance across the surface to derive expressions for ${\hat{B}_\phi }(\theta )$ and ${\hat{B}_P}(\theta )$. Consider first the toroidal magnetic field. The profiles just inside and outside the surface current are given by

(3.6)\begin{align}\left. \begin{aligned} {B_\phi }(\theta ) &= {B_0}\dfrac{{{R_0}}}{{R(\theta )}}{\textrm{Inside}\;\textrm{field,}}\\ {{\hat{B}}_\phi }(\theta ) &= {{\hat{B}}_0}\dfrac{{{R_0}}}{{R(\theta )}}{\textrm{Outside}\;\textrm{field,}} \end{aligned} \right\}\end{align}

where the only new quantity is the as yet undetermined exterior central field amplitude ${\hat{B}_0}$. Observe that ${\hat{B}_0}$ is the actual vacuum toroidal field on axis, although for mathematical convenience we formulate the analysis in terms of ${B_0}$, the vacuum field on axis in the absence of surface currents. The relationship between the outside and inside toroidal fields thus defines the poloidal surface current:

(3.7)\begin{equation}{\mu _0}\frac{{{K_P}}}{{{B_0}}} = \frac{{{{\hat{B}}_\phi }}}{{{B_0}}} - \frac{{{B_\phi }}}{{{B_0}}} = ({\varDelta _B} - 1)\frac{{{R_0}}}{R}\quad {\varDelta _B} = \frac{{{{\hat{B}}_0}}}{{{B_0}}}.\end{equation}

Hereafter, the unknown ${\varDelta _B}$ replaces ${\hat{B}_0}$ as the constant relating the two toroidal magnetic field profiles.

The relationship between the inside and outside poloidal magnetic fields is determined from the pressure balance jump condition:

(3.8)\begin{equation}{\left[\kern-0.15em\left[ {p + \frac{{{B^2}}}{{2{\mu_0}}}} \right]\kern-0.15em\right]_{\textrm{Surf}}} = 0 \to 2{\mu _0}{f_P}{p_0} + B_p^2 + B_\phi ^2 = \hat{B}_p^2 + \hat{B}_\phi ^2,\end{equation}

which can be rewritten as

(3.9)\begin{equation}{\mu _0}\frac{{{K_\phi }}}{{{B_0}}} = \frac{{{{\hat{B}}_P}}}{{{B_0}}} - \frac{{{B_P}}}{{{B_0}}} = {\left[ {\frac{{B_P^2}}{{B_0^2}} + {\beta_0}{f_P} - (\Delta_B^2 - 1)\frac{{R_0^2}}{{{R^2}}}} \right]^{1/2}} - \frac{{{B_P}}}{{{B_0}}},\end{equation}

where again ${\beta _0} = 2{\mu _0}{p_0}/B_0^2$. If for the moment we assume that ${\varDelta _B}$ is known, then (3.9) expresses the outside poloidal field in terms of known quantities.

3.3. Evaluation of the jump conditions

We are now in a position to complete the evaluation of the surface current constraint. The practical procedure is as follows.

1. a. Calculate a smooth edge equilibrium as described in Part 1 assuming the geometry and $\nu $ have been specified. Evaluate ${B_P}(\theta )/{B_0},\;{\beta _0},\;I/{B_0}$.

2. b. Specify the additional pressure pedestal and bootstrap fractions ${f_P}$ and ${f_B}$.

3. c. ‘Guess’ a value for ${\varDelta _B}$.

4. d. Evaluate ${\hat{B}_P}(\theta )/{B_0}$ from (3.9).

5. e. From (3.5) calculate

   (3.10)\begin{equation}{f_B} = 1 - \frac{I}{{\hat{I}}} = 1 - \frac{{\oint {\dfrac{{{B_P}}}{{{B_0}}}\,\textrm{d}{l_P}} }}{{\oint {\dfrac{{{{\hat{B}}_P}}}{{{B_0}}}\,\textrm{d}{l_P}} }}.\end{equation}

6. f. Iterate ${\varDelta _B}$ until ${f_B}$ is equal to the desired input value.

7. g. Evaluate ${\hat{B}_\phi }(\theta )/{B_0}$ from (3.7).

Assuming this procedure is successfully carried out, we then have the required expressions for the field jumps across the surface currents. We emphasize that no modifications have been made to the solution procedure in Part 1 that calculates the core profiles. The same expansions and constraints still apply. Once the solution for $\psi (x,y)$ is obtained, we can then make use of the jump conditions to define and evaluate modified plasma parameters consistent with standard usage in the fusion community. Several plasma parameters are calculated and compared with the corresponding no-pedestal cases obtained in Part 1. However, we delay presenting these results until the end of the analysis, after a pressure gradient pedestal, current density pedestal and toroidal flow are included.

4.1. Mathematical modifications to account for pedestals

Mathematically, the procedure to include pressure gradient and current density pedestals requires modifications to the two free functions $p(\psi )$ and ${F^2}(\psi )$. Solov'ev-like terms, linear in $\psi $, must be added. The modified free functions can be written as

(4.1)\begin{equation}\left. \begin{gathered} {p(\varPsi /{\varPsi _0}) = {p_0}\left\{ {{f_P} + (1 - {f_P}){{\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right)}^2} + {A_1}\left[ {\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right) - {{\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right)}^2}} \right]} \right\},}\\ {{F^2}(\varPsi /{\varPsi _0}) = R_0^2B_0^2\left\{ {1 + \dfrac{{2\delta B}}{{{B_0}}}{{\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right)}^2} + 2{A_2}\left[ {\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right) - {{\left( {\dfrac{\varPsi }{{{\varPsi _0}}}} \right)}^2}} \right]} \right\}.} \end{gathered} \right\}\end{equation}

As before $\varPsi = 0$ is the flux on the plasma surface and $\varPsi = {\varPsi _0}$ is the flux on the magnetic axis. Thus, $1 \ge \varPsi /{\varPsi _0} \ge 0$. As in § 3, ${f_P} = p(\textrm{surf})/p(\textrm{axis)}$ is the fractional height of the pressure profiles. The ratio ${F^2}(\textrm{axis})/{F^2}(\textrm{surf}) = 1 + 2\delta B/{B_0}$ again represents the toroidal field paramagnetism/diamagnetism. The remaining two parameters, ${A_1},{A_2}$, can for the moment be considered free. They will ultimately be related to the as yet undefined current pedestal parameter ${f_J}$.

4.2. Modified Grad–Shafranov equation

The pedestal analysis begins by introducing the same normalized quantities as previously defined: $\psi = \varPsi /{\varPsi _0}$, $R = {R_0}{r^{1/2}}$ and $Z = ay$. The Grad–Shafranov equation can now be written as

(4.2)\begin{gather}\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}}} ={-} {j_\phi }(x,\psi ),}\\ {{j_\phi }(x,\psi ) = \dfrac{{2R_0^2{a^2}}}{{\varPsi _0^2}}\left\{ {\left[ {{\mu_0}{p_0}(1 - {f_P} - {A_1})r + B_0^2\left( {\dfrac{{\delta B}}{{{B_0}}} - {A_2}} \right)} \right]\psi + \left[ {\dfrac{{{\mu_0}{p_0}}}{2}{A_1}r + \dfrac{{B_0^2}}{2}{A_2}} \right]} \right\},}\\ {r = 1 + {\varepsilon^2} + 2\varepsilon x.} \end{gathered} \right\}\end{gather}

A key step in the analysis is to recognize that by choosing ${A_2}/{A_1}$ to make the ratio of the r term to the constant term in the $\psi $ square bracket identical to the same ratio as in the constant square bracket, the result will lead to a simple particular solution for $\psi $. This pedestal constraint can be written as

(4.3)\begin{equation}\frac{{(1 - {f_P} - {A_1})}}{{(\delta B/{B_0} - {A_2})}} = \frac{{{A_1}}}{{{A_2}}} \to {A_2} = \frac{{\delta B}}{{{B_0}}}\frac{{{A_1}}}{{1 - {f_P}}}.\end{equation}

The quantity ${j_\phi }$ then reduces to

(4.4)\begin{equation}{j_\phi }(x,\psi ) = \frac{{R_0^2}}{{\varPsi _0^2}}\left( {{\mu_0}{p_0}r + \frac{{{B_0}\delta B}}{{1 - {f_P}}}} \right)[2(1 - {f_P} - {A_1})\psi + {A_1}].\end{equation}

We are now in a position to introduce the current pedestal parameter ${f_J}$. Since the current density is not a pure flux function, its pedestal height varies around the plasma surface. A simple but convenient definition that represents an average of this variation is given by

(4.5)\begin{equation}{f_J} = \frac{{{J_\phi }(R,\textrm{surf})}}{{{J_\phi }(R,\textrm{axis})}} = \frac{{{j_\phi }(x,0)/{r^{1/2}}}}{{{j_\phi }(x,1)/{r^{1/2}}}} = \frac{{{A_1}}}{{2(1 - {f_P} - {A_1}) + {A_1}}} \to {A_1} = \frac{{2(1 - {f_P}){f_J}}}{{1 + {f_J}}},\end{equation}

which leads to

(4.6)\begin{equation}{j_\phi }(x,\psi ) = \frac{{2R_0^2{a^2}}}{{\varPsi _0^2(1 + {f_J})}}[{\mu _0}{p_0}(1 - {f_P})r + {B_0}\delta B][(1 - {f_J})\psi + {f_J}].\end{equation}

Observe that the pedestal constraint introduced in (4.3) has led to a form of ${j_\phi }$ that implies that $\psi ={-} {f_J}/(1 - {f_J}) = \textrm{const}\textrm{.}$ is a particular solution to the Grad–Shafranov equation.

The final form of the Grad–Shafranov equation is obtained by substituting (4.6) into (4.2), defining modified constants $\alpha ,\nu $ plus a shifted flux function ${\psi _J}(x,y)$:

(4.7)\begin{equation}\left. \begin{gathered} {(1 + \hat{\varepsilon }x)\dfrac{{{\partial^2}{\psi_J}}}{{\partial {x^2}}} + \dfrac{1}{{1 + {\varepsilon^2}}}\dfrac{{{\partial^2}{\psi_J}}}{{\partial {y^2}}} ={-} {\alpha^2}(1 + \hat{\varepsilon }\nu x){\psi_J},}\\ {{\alpha^2} = \dfrac{{2R_0^2{a^2}}}{{\varPsi_0^2}}\left( {\dfrac{{1 - {f_J}}}{{1 + {f_J}}}} \right)\left[ {{\mu_0}{p_0}(1 - {f_P}) + \dfrac{{{B_0}\delta B}}{{1 + {\varepsilon^2}}}} \right],}\\ {\nu = \dfrac{{{\mu_0}{p_0}(1 - {f_P})}}{{{\mu_0}{p_0}(1 - {f_P}) + {B_0}\delta B/(1 + {\varepsilon^2})}},}\\ {{\psi_J} = \psi + \dfrac{{{f_J}}}{{1 - {f_J}}}.} \end{gathered} \right\}\end{equation}

This equation is formally identical to the no-pedestal equation given by (2.1). The only difference is in the flux function boundary constraint. We now require

(4.8)\begin{equation}\psi (\textrm{surf}) = 0 \to {\psi _J}(\textrm{surf}) = \frac{{{f_J}}}{{1 - {f_J}}}.\end{equation}

4.3. Boundary constraints

Consider next the boundary constraints at the matching points. Since $\psi $ and ${\psi _J}$ only differ by a constant, the slope constraints and curvature constraints are the same with or without an edge current pedestal, assuming that the model surfaces (i.e. the Miller model (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998; Todd et al. Reference Todd, Manickam, Okabayashi, Chance, Grimm, Greene and Johnson1979) or the intersecting ellipse model) also remain unchanged.

However, there is one subtlety with the curvature constraints at the inner and outer midplane points with a separatrix. Recall that without a pedestal the midplane surface shapes make use of the fact that at the X-point the two intersecting ellipses must be perpendicular to each other. With a finite edge current the angle is no longer $\mathrm{\pi }/2$ and in fact cannot be determined until the whole problem is solved. If the current pedestal is not too large, it is not important to modify the midplane curvature constraints from their no-edge-current values. The only result is that the analytic solution will not match the intersecting ellipse model very accurately at the X-point because of the differing angles. Nevertheless, the solution still satisfies the exact Grad–Shafranov equation, and there is nothing fundamental about the intersecting ellipse shape. Most importantly, the analytic solution will automatically produce the correct angle at the X-point, while the intersecting ellipse model will always have the wrong angle if the edge current is non-zero.

4.4. Solution procedure

Since the boundary condition on ${\psi _J}(\textrm{surf)}$ is now inhomogeneous, we might at first think that solutions can be found for arbitrary values of ${\alpha ^2}$. This is not the case, as there is a second normalization condition on the flux that is now non-trivial and must be satisfied, namely that

(4.9)\begin{equation}{\psi _J}(\textrm{axis}) = \frac{1}{{1 - {f_J}}}.\end{equation}

A useful way to understand the system is as follows. Assume that the geometry and a value of $\nu $ are specified. ‘Guess’ a value for ${\alpha ^2}$ and temporarily ignore the normalization condition on the axis. Because the boundary condition on the surface is inhomogeneous, it is indeed true that we can always find a non-trivial solution, subject to (4.8), that satisfies the 7 or 12 constraint relations at the required surface points. For example, for symmetric plasmas, the constraint relations are now given by

(4.10)\begin{equation}\left. \begin{gathered} {\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\leftrightarrow$}} \over A} }(\alpha )\boldsymbol{\cdot }\boldsymbol{u} = \boldsymbol{v},}\\ {\boldsymbol{u} = [{c_1},{c_2},{s_2},{c_3},{s_3},{c_4},{s_4}],}\\ {\boldsymbol{v} = [{\psi_J}(\textrm{surf}),{\psi_J}(\textrm{surf}),{\psi_J}(\textrm{surf}),0,0,0,0]\quad {\psi_J}(\textrm{surf}) = \dfrac{{{f_J}}}{{1 - {f_J}}}.} \end{gathered} \right\}\end{equation}

A similar relation holds for the asymmetric case.

However, the normalization condition given by (4.9) will not in general be satisfied. Furthermore, even though the equation is linear we cannot simply scale the solution to satisfy (4.9). If we do, then the inhomogeneous boundary condition at the surface will no longer be satisfied. In other words, one must carefully choose the value of $\alpha $ so that (4.9) just happens to be satisfied when the constraint relations are solved subject to (4.8). The parameter $\alpha $ is an output of the analysis, not an input.

The value of $\alpha $ can be determined in a manner similar to that for the no pedestal case. The 7 × 7 or 12 × 12 set of inhomogeneous constraint relations are solved for a given ‘guess’ for $\alpha $ which is then iterated until the error

(4.11)\begin{equation}E(\alpha ) = {\left[ {\frac{{{\psi _J}(\textrm{axis}) - 1/(1 - {f_J})}}{{\left| {{\psi _J}(\textrm{axis})} \right| + \left| {1/(1 - {f_J})} \right|}}} \right]^2}\end{equation}

has a minimum value at $E(\alpha ) = 0$.

4.5. Examples

To demonstrate the effects of pedestals we illustrate two examples with geometric parameters identical to the pedestal-free spherical tokamak equilibrium discussed in Part 1: $\varepsilon = 0.75,\;{\kappa _X} = 2.4,\;{\delta _X} = 0.8,\;\nu = 1.04$. The results are shown in figure 1. The first example (top profile curves) has a pedestal only in current density $({f_P} = {f_B} = 0,{f_J} = 0.25)$. The second example (bottom profile curves) has pedestals in the pressure $({f_P} = 0.2)$ and current density $({f_J} = 0.25)$. There is also a bootstrap surface current $({f_B} = 0.35)$. As stated, the flux surface plots for both examples are the same since they are not affected by the pressure pedestal or bootstrap current.

Figure 1.

Flux surface and midplane profile plots for a spherical tokamak including pedestals. Black curves correspond to the pedestal-free solution while the white curves contain pedestals. In the midplane plots each curve is normalized to unity. The upper plot has only a ${J_\phi }$ pedestal while the lower plot also has a pressure pedestal and bootstrap current. The labels I and II refer to Part 1 (no pedestals) and Part 2 (with pedestals).

In the flux surface plot, the red boundary curve shows the intersecting ellipse model surface. The thinner black curves are the flux surfaces for the pedestal-free case discussed in Part 1. They are illustrated just for purposes of comparison. The white flux surfaces represent the new solutions including the effects of a current density pedestal. The same colours (black for no-pedestal, no-flow solution, white for the new solution in this paper) are used in the remainder of the paper.

Observe that this pedestal produces a slight inward shift of the flux surfaces. A current density with the same ${q_0} = 1$ but containing a ${J_\phi }$ pedestal, has a flatter radial profile with approximately the same ${J_\phi }$ on axis. The total current has increased implying that ${J_\phi }$ is larger near the outer boundary, which pushes the pressure inwards so the shift is a little smaller. Also, recall that the slightly larger mismatch at the X-point is due to setting the X-point angle in the intersecting ellipse model surface to $\mathrm{\pi }/2$ for convenience, although this is not the correct, self-consistent value. The analytic shaded boundary edge has the correct X-point angle which therefore represents an output of the analysis.

Although the flux surfaces are the same, the two examples have different midplane pressure and current density profiles. This is clearly seen in figure 1. In both plots the edge pedestal in the current density on the outboard side has a similar height. The pedestal on both inboard sides is zero as the parameter $\nu $ has been chosen to be at its maximum value, corresponding to zero current density at $R = {R_0} - a$. The upper curves show the edge pressure smoothly going to zero while the lower curves show the finite pressure pedestal. These pedestals will have an impact on the plasma parameters of interest as shown shortly.

5.1. The starting equation

The final case of interest adds toroidal flow to the analysis. A general form of the Grad–Shafranov equation including both toroidal and poloidal flow has been known for many years (Zehrfeld & Green Reference Zehrfeld and Green1972; Maschke & Perrin Reference Maschke and Perrin1980; Morozov & Solov'ev Reference Morozov and Solov'ev1980; Hameiri Reference Hameiri1983; Semenzato et al. Reference Semenzato, Gruber and Zehrfeld1984; Iacono et al. Reference Iacono, Bondeson, Troyon and Gruber1990) and nonlinear numerical codes have been developed to solve this equation (Semenzato et al. Reference Semenzato, Gruber and Zehrfeld1984; Belien et al. Reference Belien, Botchev, Goedbloed, Van der Holst and Keppens2002; Guazzotto et al. Reference Guazzotto, Betti, Manickam and Kaye2004). For readers unfamiliar with the generalized Grad–Shafranov equation in the limit of purely toroidal flow, we present a derivation in Appendix A that leads to a form allowing analytic solutions. This somewhat complicated form is given by

(5.1)\begin{gather}\left. \begin{gathered} {{\Delta ^*}{\varPsi } = - FF' - {\mu _0}{R^2}{{\left[ {1 + {M^2}\left( {\frac{{{R^2}}}{{R_0^2}} - 1} \right)} \right]}^{1/(\gamma - 1)}}\left\{ {\bar{p}' + \left[ {{M^2}\bar{p}' + \frac{\gamma }{{\gamma - 1}}\bar{p}({M^2})'} \right]\left( {\frac{{{R^2}}}{{R_0^2}} - 1} \right)} \right\},}\\ {{B_\phi } = \frac{{F({\varPsi })}}{R},}\\ {\frac{p}{{\bar{p}({\varPsi })}} = {{\left[ {1 + {M^2}\left( {\frac{{{R^2}}}{{R_0^2}} - 1} \right)} \right]}^{\gamma /(\gamma - 1)}},}\\ {\frac{\rho }{{\bar{\rho }({\varPsi })}} = {{\left[ {1 + {M^2}\left( {\frac{{{R^2}}}{{R_0^2}} - 1} \right)} \right]}^{1/(\gamma - 1)}},}\\ {{\Omega ^2}({\varPsi })R_0^2 = \frac{{2\gamma }}{{\gamma - 1}}\frac{{\bar{p}({\varPsi })}}{{\bar{\rho }({\varPsi })}}{M^2}({\varPsi }).} \end{gathered} \right\}\end{gather}

Observe that three free functions must be specified: the familiar toroidal field and pressure functions $F(\varPsi ),\bar{p}(\varPsi )$ as well as a new function representing the square of the effective thermal Mach number ${M^2}(\varPsi ) = [2\gamma /(\gamma - 1)]M_T^2$ with ${M_T}(\varPsi ) = {(\bar{\rho }{\varOmega ^2}R_0^2/\bar{p})^{1/2}}$ the standard thermal Mach number. The density function $\bar{\rho }(\varPsi )$ and toroidal angular velocity function $\varOmega (\varPsi ){R_0}$ are separate free functions but only the combination $\bar{\rho }{\varOmega ^2}R_0^2$ appears in the Grad–Shafranov equation. Also, note that $\gamma $ is the ratio of specific heats.

We obtain analytic solutions by special choices of $F(\varPsi ),\bar{p}(\varPsi ),M(\varPsi )$ that again lead to a linear partial differential equation in $\varPsi $. The coefficient of the $\varPsi $ term can be reduced to a simple power series in ${R^2}$ which is convenient for our method of solution. To make this reduction we must restrict $\gamma $ to two values: $\gamma = 2$ corresponding to adiabatic behaviour and $\gamma = \infty $ corresponding to incompressible flow. Recall that the adiabatic index $\gamma = (N + 2)/N$, where N is the number of degrees of freedom for fluid motion. Thus, $\gamma = 2$ corresponds to two degrees of freedom, parallel and perpendicular to the magnet field. Two other interesting cases,$\gamma = 5/3,\gamma = 1$, do not lead to simple forms amenable to an analytic solution, similar to the ones presented in this paper. However, different approaches may in the future also provide analytic solutions.

The analysis proceeds in two steps. First, we treat the case of flow with smooth edge conditions in which the pressure, pressure gradient and current density all vanish smoothly at the plasma edge. Second, we add the effects of pedestals and surface currents. The solution proceeds as follows.

5.2. General reduction leading to a linear equation in $\varPsi $

We can reduce (5.1) to a linear equation in $\varPsi $ for smooth plasma edge conditions and arbitrary $\gamma $, by choosing the three free functions as follows:

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

Besides the familiar ${R_0},{B_0}$, two new parameters are introduced: ${p_0},{M_0}$. These parameters have the following qualitative interpretations. If we approximate the magnetic axis by setting $R = {R_0}$, then a short calculation shows that

(5.3)\begin{equation}\left. \begin{gathered} {{p_0} \approx p{|_{\textrm{axis}}}\quad \textrm{Pressure}\;\textrm{on}\;\textrm{axis,}}\\ {M_0^2 \approx \dfrac{{\gamma - 1}}{{2\gamma }}{{\left. {\dfrac{{\bar{\rho }{\varOmega^2}R_0^2}}{{\bar{p}}}} \right|}_{\textrm{axis}}}\quad \textrm{Mach}\;\textrm{number}\;\textrm{squared}\;\textrm{on}\;\textrm{axis}\textrm{.}} \end{gathered} \right\}\end{equation}

Observe that ${p_0}$ is defined the same way in terms of the flux function as in Part 1. Nevertheless, since the actual pressure is not an exact flux function with non-zero flow, the quantity ${p_0}$ is only approximately equal to the pressure on axis. The approximations in (5.3) are accurate for large aspect ratio. They are less accurate but still qualitatively correct for tight-aspect-ratio configurations such as the spherical tokamak where the shift between the magnetic and geometric axes can be finite.

We now substitute the free functions into (5.1). A straightforward calculation leads to the desired form of the Grad–Shafranov equation for arbitrary $\gamma $:

(5.4)\begin{equation}{\varDelta ^\ast }\varPsi ={-} \frac{{2R_0^2B_0^2}}{{\varPsi _0^2}}\left[ {\frac{{\delta B}}{{{B_0}}} + \frac{{{\mu_0}{p_0}}}{{B_0^2}}\frac{{{R^2}}}{{R_0^2}}{{\left( {1 - M_0^2 + M_0^2\frac{{{R^2}}}{{R_0^2}}} \right)}^{\gamma /(\gamma - 1)}}} \right]\varPsi = 0.\end{equation}

Observe that the term in parentheses reduces to a simple polynomial in ${R^2}$ for $\gamma = 2$ (adiabatic) or $\gamma = \infty $ (incompressible) or more generally $\gamma = n/(n - 1)$ with n an integer.

5.3. Normalized equations

As for the case of zero flow, we again introduce normalized variables:

(5.5)\begin{equation}\left. \begin{gathered} {\varPsi = {\varPsi_0}\psi \quad 1 \ge \psi \ge 0,}\\ {Z = ay\quad - \kappa \le y \le \kappa ,}\\ {R = {R_0}{r^{1/2}} = {R_0}{{(1 + {\varepsilon^2} + 2\varepsilon x)}^{1/2}}\quad - 1 \le x \le 1.} \end{gathered} \right\}\end{equation}

The modified Grad–Shafranov equation reduces to

(5.6)\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 + \nu G(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[ {\dfrac{{{B_0}\delta B}}{{1 + {\varepsilon^2}}} + {\mu_0}{p_0}{{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}}} \right] \ge 0,}\\ {\nu = \dfrac{{{\mu_0}{p_0}{{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}}}}{{{\mu_0}{p_0}{{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}} + {B_0}\delta B/(1 + {\varepsilon^2})}} \ge 0,}\\ {G(x) = (1 + \hat{\varepsilon }x){{\left( {1 + \dfrac{{2M_0^2\varepsilon }}{{1 + M_0^2{\varepsilon^2}}}x} \right)}^{\gamma /(\gamma - 1)}} - 1.} \end{gathered} \right\}\end{equation}

Observe that for small $\varepsilon $, $G(x) \propto \varepsilon x$. The magnetic axis is defined by the condition $\psi (\textrm{axis}) = 1$. Also, the range of $\nu $ to prevent current reversal on the inboard midplane is given by

(5.7)\begin{equation}\left. \begin{gathered} {0 \le \nu \le {\nu _{\max }},}\\ {{\nu _{\max }}(\varepsilon ,M_0^2) = - \frac{1}{{G( - 1)}} = \frac{1}{{1 - (1 - \hat{\varepsilon }){{(1 - \hat{\varepsilon }\zeta )}^{\gamma /(\gamma - 1)}}}},}\\ {\zeta (\varepsilon ,M_0^2) = \frac{{M_0^2(1 + {\varepsilon ^2})}}{{1 + M_0^2{\varepsilon ^2}}}.} \end{gathered} \right\}\end{equation}

5.4. Two special cases of interest

The two special cases of interest that lead to simple analytic solutions correspond to $\gamma = \infty $ (incompressible) and $\gamma = 2$ (adiabatic). In both cases (5.6) reduces to an equation that has the form

(5.8)\begin{equation}\left. \begin{gathered} {(1 + \hat{\varepsilon }x)\frac{{{\partial ^2}\psi }}{{\partial {x^2}}} + \frac{1}{{1 + {\varepsilon ^2}}}\frac{{{\partial ^2}\psi }}{{\partial {y^2}}} = - {\alpha ^2}\left[ {1 + \hat{\varepsilon }{\nu _1}x + {{\hat{\varepsilon }}^2}{\nu _2}{x^2} + {{\hat{\varepsilon }}^3}{\nu _3}{x^3}} \right]\psi \quad \psi (\textrm{surf)} = 0,}\\ {\gamma = \infty :\quad {\nu _1}\textrm{ = (1 + }\zeta )\nu \quad {\nu _2}\textrm{ = }\zeta \nu \quad {\nu _3}\textrm{ = 0, }}\\ {\gamma = 2:\quad {\nu _1}\textrm{ = (1 + 2}\zeta )\nu \quad {\nu _2}\textrm{ = }\zeta (2 + \zeta )\nu \quad {\nu _3}\textrm{ = }{\zeta ^2}\nu .} \end{gathered} \right\}\end{equation}

5.5. The solution

Following the procedure described in Part 1 we again find a solution by separation of variables:

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

where ${Y_n},{X_n}$ satisfy

(5.10)\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_1^2x + \lambda_2^2{x^2} + \lambda_3^2{x^3}){X_n} = 0,}\\ {k_n^2 = {\alpha^2} - \dfrac{{h_n^2}}{{1 + {\varepsilon^2}}},}\\ {\lambda_1^2 = \hat{\varepsilon }{\alpha^2}{\nu_1},}\\ {\lambda_2^2 = {{\hat{\varepsilon }}^2}{\alpha^2}{\nu_2},}\\ {\lambda_3^2 = {{\hat{\varepsilon }}^3}{\alpha^2}{\nu_3},} \end{gathered} \right\}\end{equation}

and $h_n^2$ are the separation constants.

The solutions for ${Y_n}$ are given by

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

The equation for ${X_n}$ is slightly more complicated than in Part 1. Recall that the solutions in Part 1 could be written in terms of Whittaker functions although this was of no great help. Instead, a straightforward analytic power series was used to obtain solutions for ${X_n}$ which was much faster and equally accurate to numerically evaluate, as compared to standard Whittaker function packages. Equally important, analytic forms for the first and second derivatives were derived.

Following the same reasoning as in Part 1, we can again obtain solutions for ${X_n}$ in terms of power series. These series have a similar but slightly more complicated recursion relation for the expansion coefficients. The main difference is that the solutions no longer have a well-known name (e.g. Whittaker functions) but this is of little consequence. The solutions for ${X_n}$ are thus written as

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

We 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 in Appendix B. Hereafter, we can assume that ${C_n},{S_n}$ are known solutions.

The conclusion from this analysis is that the flux function solution procedures described in Part 1 for smooth, double-null and single-null systems without pedestals apply directly to systems with toroidal flow. The only difference is the generalized form of the expansion functions ${C_n},{S_n}$.

5.6. Flow with pedestals

Adding edge pedestals to the calculation of flux surfaces with flow is straightforward. The results are nearly identical and only slightly more complicated than the case without flow. The starting point is (5.4), the generalized Grad–Shafranov equation including toroidal flow. To include pedestals we require modified forms of the free functions that contain Solov'ev-type additions. Specifically, for systems with toroidal flow and pedestals we choose the same free functions as without flow (see (4.1)). After introducing normalized coordinates, we see that the Grad–Shafranov equation reduces to

(5.13)\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}}} ={-} {j_\phi }(\rho ,\psi ),}\\ {j_\phi }(x,\psi ) = \dfrac{{2R_0^2{a^2}}}{{(1 + {\varepsilon^2})\varPsi_0^2}}\left\{ {\left[ {{\mu_0}{p_0}(1 - {f_P} - {A_1})g(x) + B_0^2\left( {\dfrac{{\delta B}}{{{B_0}}} - {A_2}} \right)} \right]\psi } \right.\\ \hspace{-7.5pc}\left. + \left[ {\dfrac{{{\mu_0}{p_0}}}{2}{A_1}g(x) + \dfrac{{B_0^2}}{2}{A_2}} \right] \right\},\\ {g(x) = r\frac{\rho }{{\bar{\rho }}}(1 - M_0^2 + M_0^2r) = r{{(1 - M_0^2 + M_0^2r)}^{\gamma /(\gamma - 1)}},}\\ {r = 1 + {\varepsilon^2} + 2\varepsilon x.} \end{gathered} \right\}\end{equation}

We again choose the constants ${A_1},{A_2}$ to (i) produce a simple particular solution for $\psi $ and (ii) introduce the height of the current density pedestal. The required relations are given by (4.3) and (4.5). These values are substituted into (5.13). The result is the desired form of the Grad–Shafranov equation with flow and pedestals:

(5.14)\begin{equation}\left. \begin{gathered} {(1 + \hat{\varepsilon }x)\dfrac{{{\partial^2}{\psi_J}}}{{\partial {x^2}}} + \dfrac{1}{{1 + {\varepsilon^2}}}\dfrac{{{\partial^2}{\psi_J}}}{{\partial {y^2}}} ={-} {\alpha^2}(1 + \nu G){\psi_J},}\\ {{\alpha^2} = \dfrac{{2R_0^2{a^2}}}{{\varPsi_0^2}}\left( {\dfrac{{1 - {f_J}}}{{1 + {f_J}}}} \right)\left[ {{\mu_0}{p_0}(1 - {f_P}){{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}} + \dfrac{{{B_0}\delta B}}{{1 + {\varepsilon^2}}}} \right],}\\ {\nu = \dfrac{{{\mu_0}{p_0}(1 - {f_P}){{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}}}}{{{\mu_0}{p_0}(1 - {f_P}){{(1 + M_0^2{\varepsilon^2})}^{\gamma /(\gamma - 1)}} + {B_0}\delta B/(1 + {\varepsilon^2})}},}\\ {{\psi_J}(x,y) = \psi + \dfrac{{{f_J}}}{{1 - {f_J}}}\quad {\psi_J}(\textrm{surf}) = \dfrac{{{f_J}}}{{1 - {f_J}}}\quad {\psi_J}(\textrm{axis}) = \dfrac{1}{{1 - {f_J}}},}\\ {G(x) = (1 + \hat{\varepsilon }x){{\left( {1 + \dfrac{{2M_0^2\varepsilon }}{{1 + M_0^2{\varepsilon^2}}}x} \right)}^{\gamma /(\gamma - 1)}} - 1.} \end{gathered} \right\}\end{equation}

Note that this equation is identical in form to the one with pedestals but no flow. The only differences in the equation are slightly modified definitions of $\alpha ,\nu $ and the more general form for $G(x)$ which corresponds to the generalized ${C_n}(x),{S_n}(x)$. The conclusion is that with these minor modifications we can use the same solution procedure as for the system without flow.

5.7. Examples

We have run a large number of cases with flow to convince ourselves that the solution procedure works and is robust. For due diligence wherever possible, we also compared our results with the established FLOW code to verify the numerical implementation of our code. A specific example is presented in Appendix C where it is shown that all tests proved highly satisfactory. To demonstrate this statement, we present two examples which simultaneously include all the effects of our generalized analysis in the most complicated geometry. The configurations include pedestals in pressure and current density, an edge bootstrap current, substantial flow, in a single-null diverted tokamak. The two cases have nearly identical parameters: $\varepsilon = 0.33,\;\kappa = 1.6,\;\delta = 0.4,\;{\kappa _X} = 2,\;{\delta _X}$ $= 0.5,\;\nu = 1,\,{f_B} = 0.35,\;{f_P} = 0.2,\;{f_J} = 0.25,\;{M_0} = 0.6$. The difference is that the case in figure 2 is adiabatic $(\gamma = 2)$ while the one in figure 3 is incompressible $(\gamma = \infty )$.

Figure 2.

Flux surface plot and midplane profiles for the adiabatic case $(\gamma = 2)$. Black curves correspond to the reference pedestal-free, flow-free solution while the white curves contain both effects. In the midplane plots each curve is normalized to unity. The blue and red curves represent the solutions with pedestals and flow while the black curves are the reference case with neither effect. The labels I and II refer to Part 1 (no pedestals, no flow) and Part 2 (with pedestals and flow).

Figure 3.

Flux surface plot and midplane profiles for the incompressible case $(\gamma = \infty )$. Black curves correspond to the reference pedestal-free, flow-free solution while the white curves contain both effects. In the midplane plots each curve is normalized to unity. The blue and red curves represent the solutions with pedestals and flow while the black curves are the reference case with neither effect. The labels I and II refer to Part 1 (no pedestals, no flow) and Part 2 (with pedestals and flow).

We see that the solution procedure works well when flow is included. The most pronounced effect is the not unexpected large outward shift of the flux surfaces due to the flow-generated centrifugal force. A comparison of the two examples indicates that the adiabatic case has a slightly larger shift than the incompressible case.