The Impact of Collisionality, FLR and Parallel Closure Effects on Instabilities in the Tomakak Pedestal: Numerical Studies with the NIMROD code

TL;DR

This study uses the NIMROD extended-MHD code to analyze how collisionality, FLR, and parallel closure effects influence plasma instabilities in the tokamak pedestal, comparing results with ideal-MHD models and analytic calculations.

Contribution

It provides a detailed numerical verification of NIMROD against the ELITE code and explores the effects of various physical factors on plasma stability in the tokamak pedestal.

Findings

Resistive effects are destabilizing.

Finite-Larmor-radius effects are stabilizing.

Analytic models overpredict destabilization and stabilization.

Abstract

The extended-MHD NIMROD code [C.R. Sovinec and J.R. King, J. Comput. Phys. 229, 5803 (2010)] is verified against the ideal-MHD ELITE code [H.R. Wilson et al. Phys. Plasmas 9, 1277 (2002)] on a diverted tokamak discharge. When the NIMROD model complexity is increased incrementally, resistive and first-order finite-Larmour radius effects are destabilizing and stabilizing, respectively. The full result is compared to local analytic calculations which are found to overpredict both the resistive destabilization and drift stabilization in comparison to the NIMROD computations.

Equations

$$\begin{gathered} {\mathbf{\Pi }}_{\parallel i}=m_i{n}_i{\nu }_{\parallel i}\left(\hat{\mathbf{b}}\hat{\mathbf{b}}-\frac{1}{3}\mathbf{I}\right) \\ \times \left(3\hat{\mathbf{b}}\cdot \nabla {\mathbf{v}}_i\cdot \hat{\mathbf{b}}-\nabla \cdot {\mathbf{v}}_i\right), \end{gathered}$$

$${\mathbf{q}}_{\parallel \alpha }=-n_i{\chi }_{\parallel \alpha }\hat{\mathbf{b}}\hat{\mathbf{b}}\cdot \nabla T_{\alpha },$$

$$\begin{gathered} {\mathbf{\Pi }}_{\times i}=\frac{m_i{p}_i}{4ZeB}[\hat{\mathbf{b}}\times \mathbf{W}\cdot \left(\mathbf{I}+3\hat{\mathbf{b}}\hat{\mathbf{b}}\right) \\ -\left(\mathbf{I}+3\hat{\mathbf{b}}\hat{\mathbf{b}}\right)\cdot \mathbf{W}\times \hat{\mathbf{b}}], \end{gathered}$$

$$\mathbf{W}=\nabla \mathbf{v}+\nabla {\mathbf{v}}^T-(2/3)\mathbf{I}\nabla \cdot \mathbf{v},$$

$$\mathbf{E}=-\mathbf{v}\times \mathbf{B}+\frac{\mathbf{J}\times \mathbf{B}}{n_{e}e}-\frac{\nabla p_e}{n_{e}e}+\eta \mathbf{J}+\frac{m_e}{n_e{e}^2}\frac{\partial \mathbf{J}}{\partial t},$$

$${\mathbf{q}}_{\times \alpha }=\frac{5p_{\alpha }}{2q_{\alpha }B}\hat{\mathbf{b}}\times \nabla T_{\alpha },$$

$$\omega (\omega -{\omega }_{\ast i})=-{\gamma }_I^2$$

Full text

edge localized modes, peeling-ballooning modes, extended-MHD, code verification

I Introduction

Understanding the complex physics of instabilities in the tokamak pedestal region is essential for ensuring that tokamaks achieve their highest performance. Computation plays a critical role in furthering the knowledge of edge dynamics, including advances in the understanding of type-I edge-localized modes (ELMs) Leonard et al. (2006). Nonlinear computation of ELM dynamics remains challenging, but is becoming tractable with modern computing resources Pankin et al. (2007); Sugiyama and Strauss (2010); Xu et al. (2010, 2013); Xi et al. (2013); Huijsmans and Loarte (2013); Pamela et al. (2013). ELM modeling can be roughly broken into three regimes: the early linear regime, the nonlinear regime, and an energy loss regime. Two important aspects of the linear stage are the stability boundary within an operational space; e.g., pedestal current and pressure gradient; and the linear growth rate for each toroidal-mode number; i.e., the growth-rate or mode spectrum. The dominant instability drive that sets the linear-stability boundary is from the peeling-ballooning mode (PBM) Connor et al. (1998); Snyder et al. (2002, 2004) based on the ideal-MHD model. The linear spectrum is critical to the next phase, the nonlinear regime as it seeds the dynamics. Non-ideal effects can greatly impact the growth-rate spectrum and are thus important in determining the nonlinear ELM dynamics. In the final regime, high-fidelity simulation of the energy loss requires modeling of the fast conductive and convective transport, and accurate coupling to neutrals and the wall. State-of-the-art computations are beginning to push the boundary of current capabilities to modeling that captures accurate dynamics with these additional effects that are relevant during energy loss Huijsmans and Loarte (2013); Pamela et al. (2013); Gui et al. (2014).

As the subsequent nonlinear evolution and energy loss in simulations is greatly impacted by the amplitudes of initial linear perturbations; a large source of uncertainty arises from the linear phase of ELM evolution. Given the difficulty (or impossibility) of measuring the small perturbations that seed an ELM, a reasonable model is to seed it from unstable PBMs that grow from small-wave perturbations. Although the stability boundary is well described by ideal MHD, details of the mode spectrum depend upon specific non-ideal effects. In particular, resistive effects may destabilize and first-order finite-Larmour-radius (FLR) effects may stabilize the modes. These latter effects, part of which is drift stabilization, require two-fluid modeling with advanced FLR closures. It is important to verify that codes accurately capture these effects within regimes of interest. In this paper, we study these effects on the linear modes with the extended-MHD NIMROD code Sovinec et al. (2004); Sovinec and King (2010).

Type-I ELM stability boundaries are dominated by drives described by ideal MHD as the ideal-MHD operator is stiff; i.e., near the ideal-MHD stability boundary small changes in pressure and current-density gradients lead to large changes in the growth rate due to the ideal-MHD terms becoming large. Thus only small profile changes are required to cross the boundary, and non-ideal modifications to the growth rate become secondary when describing the boundary location within operational phase space. This ideal-MHD stiffness is a general property that impacts the behavior of a wide class of plasma-instability phenomena. For example, in addition to ELM stability, it also describes the dynamics of disruptions caused by core-mode activity Callen et al. (1999); Kruger et al. (2005). However, type-I ELMs are more complicated than core-mode disruptions given the broad toroidal-mode spectrum and associated nonlinear coupling. Similar to studies of core-mode disruptions, comprehensive understanding of ELM onset requires evolving the equilibrium on the transport time scale through the instability boundary with a verified code that includes non-ideal terms to calculate an accurate initial linear-mode spectrum.

MHD codes that study ELM instabilities can be placed into two broad categories: 1) extended-MHD codes (e.g. NIMROD and M3D-C1 Jardin et al. (2007)) and 2) reduced-MHD-based codes (e.g. BOUT++ Dudson et al. (2009) and JOREK Huysmans et al. (2009)). The definition of extended MHD used in this work includes the full ideal-MHD force operator (see Ref. Schnack et al. (2006) and references therein). Using the full operator retains the fast compressional-Alfvén waves that enforce force balance and thus are critical to extended-MHD transport calculations with sources and sinks (e.g. Jenkins et al. (2010)). Relative to reduced modeling, the extended-MHD approach also has advantages in avoiding low-$\beta$ approximations and a straight-forward implementation of the FLR closures as discussed in Section III. However, the inclusion of the full force operator comes at a cost of increased numerical difficulties when solving the equations associated with the fast waves Gruber and Rappaz (2012); Degtyarev and Medvedev (1986); Bondeson and Fu (1991); Lütjens and Luciani (1996); Sovinec (2015). Verification for all regimes is important, and previous work Burke et al. (2010) verified the ideal-MHD terms for a non-diverted, peeling-ballooning unstable test case. As one expects that transport calculations across the ELM stability boundary coupled with an instability model that captures FLR drift stabilization significantly impacts dynamics, emphasis on verifying linear physics of extended-MHD codes thus naturally follows.

In this work, a verification of the full extended-MHD, initial-value NIMROD code with the ideal-MHD ELITE code Wilson et al. (2002); Snyder et al. (2002, 2004) is performed using the Meudas-1-case comparison Snyder et al. (2009). This is a diverted tokamak case that is peeling-ballooning unstable. As we later discuss, the inclusion of the x-point significantly modifies the mode structure; making the comparison in full magnetic geometry an important additional verification. To assess the impact of the two-fluid, FLR effects, we compare our results with drift stabilization to a calculation of drift stabilization on the resistive PBM from Ref. Hastie et al. (2003). Although the approximations made during the analysis of first-order FLR PBMs in Ref. Hastie et al. (2003) preclude quantitative comparison, we make a qualitative comparison of our first-order FLR computations to analytics.

The paper proceeds as follows: In Sec. II we summarize the benchmark and give an overview of our methodology. The exhaustive details of the methods used in this comparison, in particular the mechanism whereby a vacuum region is effectively included through extended-MHD modeling, are previously available in Refs. Burke et al. (2010) and Ferraro et al. (2010) and thus our description of this aspect is brief. In Sec. III the model complexity is gradually increased from the ideal-like model through a resistive-MHD model with experimentally relevant profiles to a full extended-MHD model which includes both parallel-closure and FLR-drift effects. The resulting impact on the growth-rate spectrum in analyzed. Comparisons of the computed drift-stabilized growth rates are made to analytic treatments of the PBM in Sec. IV. In Sec. V the density and temperature are varied at fixed plasma $\beta$ (thus the ion gyroradius, $\rho_{i}$, is modified) in order to ascertain how the FLR effects are impacted. In this study, the bootstrap current is fixed (the zeroth-order effect from variations of density modify the bootstrap current) and the impact of FLR drift effects are isolated. Sec. VI summarizes this work and places it within the broader context of modeling of ELM relaxation events in the tokamak edge with transport effects.

II Verification benchmark

We begin with a study of a high resolution, lower-single-null, JT-60U-like equilibrium (‘Meudas-1’), which was originally employed in a benchmark of the MARG2D and ELITE codes Aiba et al. (2007), including a close approach to the X-point Snyder et al. (2009). This extends previous benchmarks Burke et al. (2010) of ELITE and NIMROD as it includes diverted magnetic topology and a higher edge safety factor ($q_{95}=6.74$, the safety factor at 95% of the normalized poloidal flux) that leads to increased resolution requirements. An ideal-MHD limit is achieved in NIMROD by using flat density and resistivity profiles inside the last closed flux surface (LCFS) with small resistivity, $S=10^{8}$ where $S$ is the Lundquist number ($S=\tau_{R}/\tau_{A}$), $\tau_{A}$ is the Alfvén time ($\tau_{A}=R_{o}/v_{A}$), $v_{A}$ is the Alfvén velocity ($B/\sqrt{m_{i}n_{i}\mu_{0}}$), $\tau_{R}$ is the resistive diffusion time ($\tau_{R}=R_{o}^{2}\mu_{0}/\eta$), $R_{o}=2.936m$ is the radius of the magnetic axis, $\eta$ is the electrical resistivity, $\mu_{0}$ is the permeability of free space, $m_{\alpha}$ is a species mass (the $\alpha$ subscript denotes ions or electrons in this work), and $n_{\alpha}$ is a species density. The deuteron mass ($m_{i}=3.34\times 10^{-27}kg$) is used. In order to reproduce the vacuum response model outside the LCFS that is used by ELITE, a low density ($0.01$ of the core density) and high resistivity ($10^{7}$ times the core resistivity) is prescribed beyond the LCFS (more details on these approximations are in Ref. Burke et al. (2010)).

The normalized growth rates ($\gamma\tau_{A}$ where the linearized mode grows as $\text{exp}[\gamma t]$) vs. toroidal mode number ($n_{\phi}$) from NIMROD and ELITE are plotted in Fig. 1. There is good agreement between the codes except for $n_{\phi}$=4 where there is a 27% relative difference. All other cases have a relative difference of less than 8% with typical differences of 5%. The NIMROD convergence in terms of the maximum polynomial order of the spectral elements is shown in Fig. 2. Convergence is most challenging at high wavenumbers where the resolution requirements are most stringent (the poloidal mesh is composed of $72\times 512$ spectral elements). These cases converge from the unstable side where the growth rate decreases with enhanced resolution. Thus the excellent agreement between NIMROD and ELITE at high $n_{\phi}$ in Fig. 1 may be spurious and indicate that slightly more resolution is required for $n_{\phi}$>25, however, the shown growth rates are likely within 5% of their converged values. Studying nearly ideal cases with extended MHD codes such as NIMROD is challenging given the vanishingly small dissipation operators, and convergence is achieved more quickly with the additional non-ideal terms in the extended-MHD equations, as in the cases in Sec. III.

Figure 3 shows a poloidal cross section of the magnetic ($B_{R}$) eigenmode. The mode develops an ‘interference-pattern’ structure near the X-point when inboard and outboard finger-like structures overlap. The finite-element-mesh nodes are superimposed atop the smallest-scale sub-figure. As established by Fig. 2, this simulation is spatially and temporally converged. The high resolution required to resolve these high-$q_{95}$, diverted cases motivated development of memory-scaling improvements in the NIMROD code.

III Increasing the model complexity

Initial-value extended-MHD computations, while more computationally expensive than ideal-MHD calculations with ELITE, are able to include additional effects that are more representative of experimental conditions:

• •

Instead of a magneto-static vacuum outside the LCFS, NIMROD includes a cold-plasma region.

• •

Finite-dissipation effects such as resistivity, viscosity and thermal conduction.

• •

Density and temperature profiles that prescribe the associated drifts and resistivity profile.

• •

Parallel-closure effects that represent large diffusivities oriented along the magnetic field.

• •

First-order FLR-closure and two-fluid effects.

Reference Ferraro et al. (2010) investigates the first three of these effects with the Meudas-1 case. We further this work by including the last two effects. The first three effects are also re-investigated and different modifications to the growth rates are found. These differences are expected as slightly different profiles for density and temperature are used in order to avoid spurious electrostatic modes with the two-fluid, first-order FLR model. Ref. Ferraro et al. (2010) avoided these modes by employing only the single-fluid resistive-MHD model.

In calculations beyond the ideal model, the density profile is $n_{e}(\psi_{n})=n_{e0}(p(\psi_{n})/p_{0})^{0.3}$ where $n_{e0}=6\times 10^{19}\;m^{-3}$ and $\psi_{n}$ is the normalized poloidal flux. Here $p$ is the pressure as prescribed by the ideal gas law, $p=n_{e}T_{e}+n_{i}T_{i}$, where $T_{\alpha}$ is a species’ temperature. This choice does not match Ref. Ferraro et al. (2010), but effectively avoids spurious electrostatic modes with the two-fluid, first-order FLR model. These spurious modes are discussed in more detail in Sec. V. The pressure-profile gradient is specified by the Meudas-1 case and setting the edge temperature to $200eV$ leads to a core temperature of $5806eV$. The edge temperature is modified independent of the MHD-stability (determined by $p^{\prime}$) by the transformation $p(\psi_{n})\rightarrow p(\psi_{n})+p_{c}$ where $p_{c}$ is a constant. Without scrape-off-layer (SOL) profiles of density and temperature, which are not included in standard reconstructions, this choice of the edge temperature also specifies the temperature outside the LCFS. SOL profiles are required to get both the correct profiles outside the LCFS, which determines the vacuum response, and simultaneously set the correct resistivity at the mode resonant location when Spitzer resistivity is used. The non-ideal viscous, conductive and particle diffusivities are set to a small value, one tenth of the resistivity at the LCFS.

Figure 4 plots the growth rates computed by NIMROD on the ‘Meudas-1’ case with four different models: ideal MHD (ideal), resistive MHD with reconstructed density, temperature, and Spitzer-resistivity profiles (1f), the single-fluid model with parallel Braginskii closures (1f-Par) where the parallel viscosity is

[TABLE]

and the parallel heat flux vector is

[TABLE]

and a two-fluid model with parallel closures that includes ion gyroviscosity

[TABLE]

where $\mathbf{W}$ is the rate-of-strain tensor

[TABLE]

a full extended-MHD Ohm’s law,

[TABLE]

and separate temperature evolutions with cross-heat fluxes,

[TABLE]

(2f-Par-FLR). Here $v_{i}$ is the bulk-ion flow, $\nu_{\parallel i}$ is the ion parallel diffusivity, $\chi_{\parallel\alpha}$ is the species’ parallel diffusivity, $q_{\alpha}$ ($e$) is the species’ (electron) electric charge, $\mathbf{J}$ is the current density, $\hat{\mathbf{b}}$ is the magnetic unit vector ($\hat{\mathbf{b}}=\mathbf{B}/B$), and $\mathbf{I}$ is the identity tensor. An effective ion charge ($Z$) of three is assumed. Each model builds on the previously listed and includes all prior terms. Electron viscosity is not included.

When a density profile is included, the normalized growth rate of the mode remains constant with the single-fluid model when the Alfvén time is computed at the mode resonant surface. Our calculations are normalized by the Alfvén time as computed with the values at the magnetic axis. These values remain constant when a density profile is included such that the normalized mode growth rate is effectively increased.

The effect of the resistivity profile is more nuanced. Finite resistivity allows for reconnection within the model by relaxing the frozen-flux constraint. This permits resistive-ballooning modes that grow faster than their ideal counterparts. Alternatively, the resistive dissipation can stabilize the modes by acting as a dissipative term and modifying the response outside the LCFS from that of a vacuum to that of a plasma. Figure 5 of Ref. Ferraro et al. (2010) shows that the plasma response is stabilizing relative to a vacuum region. Overall, these effects stabilize the ballooning modes at low-$n_{\phi}$ and destabilize the modes at high-$n_{\phi}$ as seen in Fig. 4 when comparing the single fluid and ideal normalized growth rates. The effect of including small particle, viscous, and thermal diffusivities on the resistive-MHD mode is small (not shown).

Including the Braginskii parallel closures with large coefficients ($\chi_{\parallel e}=10^{9}\times\eta_{0}/\mu_{0}$ and $\chi_{\parallel i}=\nu_{\parallel i}=10^{8}\times\eta_{0}/\mu_{0}$) has little effect on the growth rates. When thermal conduction is large, the temperature quickly equilibrates along field lines to produce an isothermal response. Thus NIMROD computations that show little effect from the parallel closures are consistent with the ELITE results that find the shape of the mode growth-rate spectrum is largely not sensitive to the value of the ratio of specific heats, $\Gamma$, as seen in Fig. 1. When the response to the compressible motion associated with the sound wave is eliminated ($\Gamma=0$ in the ELITE computations), the growth rate is slightly enhanced relative to the adiabatic limit. If the growth rates are modified beyond this small effect, it would be an indication that changes in the energy equation can impact the MHD response, and thus that closure effects are significant. Although the short-mean-free path Braginskii-like closure Braginskii (1965); Catto and Simakov (2004) is used in our NIMROD compuatations, we note that other closures, such as long-mean-free path Ramos (2011) or the general approach of solving drift-kinetic equation Held et al. (2015) for a closure can be applied. However, for this case where the parallel-closure effects are negligible, different parallel closures are unlikely to significantly modify the results.

With the full two-fluid, FLR model, there is a stabilizing effect on the intermediate and high-$n_{\phi}$ modes as expected from analytic treatments Hastie et al. (2003). A common reduced-MHD approximation for first-order FLR closures is to assume that the ion-gyroviscous force and the advection by the diamagnetic drift exactly cancel (this is colloquially known as the gyroviscous cancellation, see, e.g. Coppi (1964)). However, as discussed by Ramos in Ref. Ramos (2007), “these cancellations are only partial and not very useful in practice for general magnetic geometries…”. This cancellation is only valid with a large, uniform guide field, without curvature terms (slab approximation), and in the low-$\beta$ limit (see, e.g. King and Kruger (2014)). Use of the full gyroviscous operator is critical as in addition to effects that appear as diamagnetic drifts, it also contains terms proportional to first-order FLR magnetic-curvature and grad-B drifts King et al. (2011). Our modeling contains not just the first-order FLR drifts from ion gyroviscosity, but also the drifts from the cross-heat fluxes that enter the equations on the same order.

Just as the large balance of the ideal-MHD forces can cause numerical difficulties with the ideal-MHD force operator, the partial gyroviscous cancellation leads to similar numerical pitfalls and verification is important. The implementation of these terms in NIMROD are verified in cylindrical and slab magnetic geometries against analytic calculations for tearing Sovinec and King (2010); King et al. (2011); King and Kruger (2014) and ion-temperature-gradient Schnack et al. (2013) (ITG) modes. To test the effect of these terms against analytic theory in full tokamak magnetic geometry, we perform a qualitative comparison with analytic theory in Sec. IV. Full verification is precluded by the approximations made in the analytic theories.

The first-order FLR model is valid for the toroidal modes presented although this is not necessarily the case for general edge mode modeling. Assuming $k\simeq(q/r+1/R)n_{\phi}$ and evaluating local quantities at the peak of the mode eigenfunction ($\psi_{n}=0.969$) on the outboard midplane, the normalized ion-gyroradius ($\rho_{i}=\sqrt{\Gamma m_{i}T_{i}}/ZeB$) as a function of toroidal mode is $k\rho_{i}\simeq 0.0055n_{\phi}$. If the first-order FLR assumption is violated, more-complex full-ion-orbit kinetic modeling is required to simulate large-fluctuation ELM dynamics. For this case, computations produce approximately the same results with and without the cross-heat fluxes. This is consistent with drift analytics King and Kruger (2014) that shows that the cross-heat flux is only significant at very low values of plasma $\beta$ (here the $\beta=2\mu_{0}p/B^{2}$ is 0.1%). Perhaps coincidentally for this case, the effects of the resistive destabilization and drift stabilization largely counteract one another. The drift stabilized growth rates are not less than the ideal calculation until $n_{\phi}\gtrsim 24$.

IV Comparison to drift analytics

To gain insight into our two-fluid computations, we compare with analytic descriptions of drift stabilization of the ballooning mode. The dispersion relation for drift stabilizations of the ideal, incompressible ballooning mode is

[TABLE]

where $\gamma_{I}$ is the ideal growth rate in the absence of drift effects and $\omega$ is the complex frequency ($\omega=i\gamma$) Tang et al. (1982). The diamagnetic drift velocity is typically defined as the perpendicular-to-the-magnetic-field flow contribution from the pressure gradient within the context of a two-fluid, or generalized Ohm’s law. For a tokamak, neoclassical effects damp the poloidal diamagnetic flow contribution. Thus in our case the remaining toroidal diamagnetic flow determines the diamagnetic frequency, $\omega_{*i}=\mathbf{k}\cdot\mathbf{v}_{*i}=n_{\phi}/(n_{i}q_{i})\partial p_{i}/\partial\psi$ where $\mathbf{k}$ is the mode wavenumber. Figure 5 shows the comparison to this dispersion relation with full and reduced diamagnetic drift values. The value of $\omega_{*i}$ is computed at the radial location where the mode structure peaks. The normalized frequency from the diamagnetic drift velocity, $\omega_{*i}\tau_{A}$, varies linearly with toroidal mode number as $0.00483n_{\phi}$. This ‘local’ approximation, which neglects the radial variation of the $\omega_{*i}$ profile, is known to over-estimate the effect of drift stabilization Hastie et al. (2000); Snyder et al. (2011). We find that a reduction of the value of $\omega_{*i}$ (approximately by a factor of 2-4) gives the best agreement between Eqn. (7) and the two-fluid NIMROD calculation. A similar comparison of ELITE to the two-fluid BOUT++ model on a different case finds the drift stabilization is over-predicted Snyder et al. (2011) by a comparable factor.

One of the limitations of Eqn. (7) is that it does not include the effect of resistive destabilization on the high-$n_{\phi}$ modes, as is included in our two-fluid computations. The analytic treatment of Ref. Hastie et al. (2003) provides a dispersion relation (Eqn. (43) of the reference; referred to as HRP Eqn. (43) in this discussion) that includes the effect of resistive destabilization in addition to drift effects. In order to compare with this theoretical treatment, parameters are computed on the outboard midplane at the radial location where the mode structure peaks, consistent with the previous $\omega_{*i}$ calculation. This yields local values of the safety factor, $q=7.34$, normalized magnetic shear, $s=8.8$ and normalized pressure gradient, $\alpha=28.8$, where the definitions of Ref. Hastie et al. (2003) are used. The $\Delta^{\prime}_{B}$ stability parameter is determined by using the ideal growth rate ($\Gamma=5/3$) and solving Eqn. 35 of Ref. Hastie et al. (2003). $\Delta^{\prime}_{B}$ varies with toroidal mode number but falls in the range of $-17.4$ to $-12.1$.

Comparison of HRP Eqn. (43) with $\omega_{*}=0$ and the NIMROD resistive-MHD results differ by a factor approximately 2.5 at large toroidal mode numbers where the analytic calculation results in greater growth rates and there is no stabilization at low $n_{\phi}$. As the analytic calculation is a local calculation, it does not capture the stabilizing influence on the global mode of finite resistivity, in particular it misses the effect of modifying the response outside the LCFS from that of a vacuum to that of a resistive plasma. Our interest is largely to compare and contrast the predicted drift stabilization and thus we reduce the resistivity to 35% of the value used in the NIMROD computations when calculating the analytic values of HRP Eqn. (43). This produces growth rates with $\omega_{*}=0$ that roughly match the NIMROD resistive-MHD computations at high $n_{\phi}$ as shown in Fig. 6.

More significantly, Fig. 6 compares the results of our linear NIMROD computations with the drift-stabilized growth rates computed from Eqn. (43) of Ref. Hastie et al. (2003). Again, the analytics of HRP Eqn. (43) over-estimate the effect of drift stabilization where the reduction of $\omega_{*}$ that leads HRP Eqn. (43) to qualitatively agree with the two-fluid NIMROD computations is between 2 to 4 (dependent on $n_{\phi}$). There is little effect on the low-$n_{\phi}$ modes and the ion drift-wave resonance predicted in Ref. Hastie et al. (2003) is not observed.

V FLR parameter scan

In order to ascertain the role of two-fluid effects further, we vary the density and temperature while scaling parameters to maintain constant $\beta$ and $S$ (other dissipation parameters are scaled relative to resistivity). In experimental discharges, the zeroth-order effect on the edge of modifying the density is modify the plasma collisionality and thus the bootstrap current. This causes a transition from ballooning-like modes at high density to peeling-like modes at low density. Our computations use a fixed equilibrium current and thus are not sensitive to this effect. Instead, we isolate the effect of the modification of $\rho_{i}$ (and the associated ion skin depth, $d_{i}$) within the context of our two-fluid, first-order FLR model. This contrasts with Ref. Xu et al. (2014) that examines these effects in concert and Ref. Zhu et al. (2012) that considers only the effect of modifications to the current profile.

Figure 7 shows the linear growth rate from NIMROD computations for resistive-MHD and three two-fluid cases with varying densities (the core density is varied from $3\times 10^{19}\;m^{-3}$ to $1.2\times 10^{20}\;m^{-3}$ where $6\times 10^{19}\;m^{-3}$ is the value used in the prior computations discussed in the paper). The effect of drift-stabilization is stronger at low densities ($\rho_{i}$ scales as $1/\sqrt{n_{i}}$ at constant $\beta$) as expected. The cases with $n_{e}=3\times 10^{19}\;m^{-3}$ are not plotted for $n_{\phi}>18$ as the most unstable mode is no longer the peeling-ballooning mode but rather a dominatly electro-static mode related to the ITG as studied in the context NIMROD two-fluid advance in Ref. Schnack et al. (2013). This region is excluded as examination of ITG mode dynamics is outside the scope of this work.

VI Discussion and Conclusions

Our comparisons to the ideal PBM growth rates computed by the ELITE code and analytic descriptions of the resistive, two-fluid PBM provide another verification test of the NIMROD code as a means to study the tokamak edge. Quantitative agreement with ELITE is achieved for these cases, extending prior work Burke et al. (2010) that did not include diverted magnetic topology.

Two-fluid and FLR-drift effects can both produce stabilizing (drift stabilization through the sheared motion of the ion and electron fluid responses) and destabilizing (species decoupling and enhanced growth through mediation by the more mobile electron fluid) effects (see, e.g., King et al. (2011); King and Kruger (2014)). At intermediate and large mode numbers, the inclusion of finite-resistivity is destabilizing while the FLR-drift effects are stabilizing. These effects are largely offsetting, resulting in a toroidal-mode growth-rate spectrum with the full extended-MHD model that resembles the ideal calculation. This result is not general but rather case dependent. Consistent with ELITE calculations that vary the ratio of specific heats from the adiabatic through the isothermal limit, we find the addition of anisotropic thermal conduction produces only a small modification of the growth rate.

Relative to the NIMROD calculations, analytics Hastie et al. (2003) over-predict both the effects of resistive destabilization and FLR-drift stabilization for the case studied. Qualitative agreement with the analytic descriptions of the peeling-ballooning mode is found where high-$n_{\phi}$ modes are susceptible to resistive destabilization and drift stabilization in both the analytics and NIMROD computations. One of the limits of the analytic derivation is a ‘local’ approximation which leads to the over prediction of FLR-drift stabilization Hastie et al. (2000). Given the highly complex nature of both the configuration of the tokamak edge and the two-fluid model equations, quantitative agreement is not expected as the local approximation neglects the global configuration and profile effects.

Regarding our investigation into the impact of the size of the ion gyroradius within the context of NIMROD’s full first-order FLR, two-fluid model, we find significant changes on the growth rates of high-$n_{\phi}$ modes at varying densities but constant $\beta$. However, the most unstable mode number is largely unchanged. Thus we conclude the zeroth-order effect of modifications to the bootstrap current largely dominates the determination of the most unstable mode during collisionality scans while lower densities (large $\rho_{i}$) increase the magnitude of the drift stabilization on the high-$n_{\phi}$ modes consistent with the discussion of Ref. Snyder et al. (2007).

Acknowledgements.

We thank Carl Sovinec, Chris Hegna and Nate Ferraro for discussions involving this paper and N. Aiba for providing the Meudas-1 equilibrium. This material is based on work supported by US Department of Energy, Office of Science, Office of Fusion Energy Sciences under award numbers DE-FC02-06ER54875 and DE-FC02-08ER54972. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.