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

1.7 MB·

by  J. R. King

The Impact of Collisionality, FLR and Parallel Closure Effects on Instabilities
in the Tomakak Pedestal: Numerical Studies with the NIMROD code
J. R. King,1 A. Y. Pankin,1 S. E. Kruger,1 and P. B. Snyder2
1Tech-X Corporation, 5621 Arapahoe Ave. Boulder, CO 80303, USA
2General Atomics, PO Box 85608, San Diego, CA 92186–5608, USA
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. Plas-
mas 9, 1277 (2002)] on a diverted tokamak discharge. When the NIMROD model complex-
7
ity is increased incrementally, resistive and first-order finite-Larmour radius effects are destabi-
1
lizing and stabilizing, respectively. The full result is compared to local analytic calculations
0
which are found to overpredict both the resistive destabilization and drift stabilization in com-
2
parison to the NIMROD computations. Published version: Phys. Plasmas 23, 062123 (2016)
n [http://dx.doi.org/10.1063/1.4954302]
a
J
PACSnumbers: 52.30.Ex52.35.Py,52.55.Fa,52.55.Tn,52.65.Kj
1 Keywords: edgelocalizedmodes,peeling-ballooningmodes,extended-MHD,codeverification
3

] I. INTRODUCTION uncertaintyarisesfromthelinearphaseofELMevo-
h
lution. Giventhedifficulty(orimpossibility)ofmea-
p
- suring the small perturbations that seed an ELM, a
m Understanding the complex physics of instabili- reasonable model is to seed it from unstable PBMs
ties in the tokamak pedestal region is essential for
s thatgrowfromsmall-waveperturbations. Although
a ensuring that tokamaks achieve their highest per- the stability boundary is well described by ideal
l formance. Computation plays a critical role in fur-
p MHD, details of the mode spectrum depend upon
. thering the knowledge of edge dynamics, includ- specific non-ideal effects. In particular, resistive ef-
s
ing advances in the understanding of type-I edge-
c fects may destabilize and first-order finite-Larmour-
i localizedmodes(ELMs)[1]. Nonlinearcomputation radius(FLR)effectsmaystabilizethemodes. These
s
ofELMdynamicsremainschallenging,butisbecom-
y latter effects, part of which is drift stabilization, re-
h ing tractable with modern computing resources [2– quire two-fluid modeling with advanced FLR clo-
p 8]. ELM modeling can be roughly broken into sures. Itisimportanttoverifythatcodesaccurately
[ threeregimes: theearlylinearregime, thenonlinear
capture these effects within regimes of interest. In
regime, and an energy loss regime. Two important
1 thispaper,westudytheseeffectsonthelinearmodes
v aspectsofthelinearstagearethestabilityboundary with the extended-MHD NIMROD code [13, 14].
2 within an operational space; e.g., pedestal current
Type-I ELM stability boundaries are dominated
4 and pressure gradient; and the linear growth rate
by drives described by ideal MHD as the ideal-
0 foreachtoroidal-modenumber; i.e.,thegrowth-rate
MHD operator is stiff; i.e., near the ideal-MHD
0
or mode spectrum. The dominant instability drive
0 stability boundary small changes in pressure and
that sets the linear-stability boundary is from the
. current-density gradients lead to large changes in
2 peeling-ballooningmode(PBM)[9–11]basedonthe
the growth rate due to the ideal-MHD terms be-
0 ideal-MHDmodel. Thelinearspectrumiscriticalto
7 coming large. Thus only small profile changes are
the next phase, the nonlinear regime as it seeds the
1 requiredtocrosstheboundary, andnon-idealmodi-
dynamics. Non-ideal effects can greatly impact the
: ficationstothegrowthratebecomesecondarywhen
v growth-ratespectrumandarethusimportantinde-
describingtheboundarylocationwithinoperational
Xi terminingthenonlinearELMdynamics. Inthefinal phase space. This ideal-MHD stiffness is a general
regime, high-fidelity simulation of the energy loss
r property that impacts the behavior of a wide class
a requires modeling of the fast conductive and con-
of plasma-instability phenomena. For example, in
vective transport, and accurate coupling to neutrals
addition to ELM stability, it also describes the dy-
and the wall. State-of-the-art computations are be-
namics of disruptions caused by core-mode activity
ginningtopushtheboundaryofcurrentcapabilities
[15, 16]. However, type-I ELMs are more compli-
to modeling that captures accurate dynamics with
cated than core-mode disruptions given the broad
these additional effects that are relevant during en-
toroidal-mode spectrum and associated nonlinear
ergy loss [7, 8, 12].
coupling. Similar to studies of core-mode disrup-
Asthesubsequentnonlinearevolutionandenergy tions, comprehensive understanding of ELM onset
lossinsimulationsisgreatlyimpactedbytheampli- requires evolving the equilibrium on the transport
tudesofinitiallinearperturbations;alargesourceof time scale through the instability boundary with a
2
verified code that includes non-ideal terms to calcu- modelwhichincludesbothparallel-closureandFLR-
late an accurate initial linear-mode spectrum. drift effects. The resulting impact on the growth-
MHD codes that study ELM instabilities can be ratespectruminanalyzed. Comparisonsofthecom-
placedintotwobroadcategories: 1)extended-MHD puted drift-stabilized growth rates are made to ana-
codes (e.g. NIMROD and M3D-C1 [17]) and 2) lytictreatmentsofthePBMinSec.IV.InSec.Vthe
reduced-MHD-based codes (e.g. BOUT++ [18] and density and temperature are varied at fixed plasma
JOREK[19]). ThedefinitionofextendedMHDused β (thus the ion gyroradius, ρi, is modified) in order
inthisworkincludesthefullideal-MHDforceopera- to ascertain how the FLR effects are impacted. In
tor (see Ref. [20] and references therein). Using the thisstudy,thebootstrapcurrentisfixed(thezeroth-
full operator retains the fast compressional-Alfvén order effect from variations of density modify the
waves that enforce force balance and thus are crit- bootstrap current) and the impact of FLR drift ef-
ical to extended-MHD transport calculations with fectsareisolated. Sec.VIsummarizesthisworkand
sources and sinks (e.g. [21]). Relative to reduced places it within the broader context of modeling of
modeling, theextended-MHDapproachalsohasad- ELM relaxation events in the tokamak edge with
vantages in avoiding low-β approximations and a transport effects.
straight-forwardimplementationoftheFLRclosures
as discussed in Section III. However, the inclusion
of the full force operator comes at a cost of in- II. VERIFICATION BENCHMARK
creasednumericaldifficultieswhensolvingtheequa-
tions associated with the fast waves [22–26]. Veri-
Webeginwithastudyofahighresolution,lower-
fication for all regimes is important, and previous
single-null, JT-60U-like equilibrium (‘Meudas-1’),
work [27] verified the ideal-MHD terms for a non-
whichwasoriginallyemployedinabenchmarkofthe
diverted, peeling-ballooning unstable test case. As
MARG2D and ELITE codes [32], including a close
one expects that transport calculations across the
approach to the X-point [29]. This extends previ-
ELM stability boundary coupled with an instability
ous benchmarks [27] of ELITE and NIMROD as it
model that captures FLR drift stabilization signifi-
includes diverted magnetic topology and a higher
cantly impacts dynamics, emphasis on verifying lin-
edge safety factor (q = 6.74, the safety factor at
95
ear physics of extended-MHD codes thus naturally
95% of the normalized poloidal flux) that leads to
follows.
increased resolution requirements. An ideal-MHD
In this work, a verification of the full extended- limit is achieved in NIMROD by using flat density
MHD, initial-value NIMROD code with the ideal- andresistivityprofilesinsidethelastclosedfluxsur-
MHDELITEcode[10,11,28]isperformedusingthe face (LCFS) with small resistivity, S =108 where S
Meudas-1-case comparison [29]. This is a diverted is the Lundquist number (S = τ /τ ), τ is the
R A A
tokamak case that is peeling-ballooning unstable. Alfvén time (τ = R /v ), v is the Alfvén veloc-
√ A o A A
As we later discuss, the inclusion of the x-point sig- ity (B/ m n µ ), τ is the resistive diffusion time
i i 0 R
nificantly modifies the mode structure; making the (τ = R2µ /η), R = 2.936m is the radius of the
R o 0 o
comparison in full magnetic geometry an important magneticaxis,η istheelectricalresistivity,µ isthe
0
additional verification. To assess the impact of the permeabilityoffreespace,m isaspeciesmass(the
α
two-fluid, FLR effects, we compare our results with α subscript denotes ions or electrons in this work),
drift stabilization to a calculation of drift stabiliza- and n is a species density. The deuteron mass
α
tion on the resistive PBM from Ref. [30]. Although (m =3.34×10−27kg)isused. Inordertoreproduce
i
theapproximationsmadeduringtheanalysisoffirst- thevacuumresponsemodeloutsidetheLCFSthatis
order FLR PBMs in Ref. [30] preclude quantitative used by ELITE, a low density (0.01 of the core den-
comparison,wemakeaqualitativecomparisonofour sity) and high resistivity (107 times the core resis-
first-order FLR computations to analytics. tivity) is prescribed beyond the LCFS (more details
The paper proceeds as follows: In Sec. II we sum- on these approximations are in Ref. [27]).
marize the benchmark and give an overview of our The normalized growth rates (γτ where the lin-
A
methodology. The exhaustive details of the meth- earized mode grows as exp[γt]) vs. toroidal mode
odsusedinthiscomparison, inparticularthemech- number (n ) from NIMROD and ELITE are plot-
φ
anism whereby a vacuum region is effectively in- ted in Fig. 1. There is good agreement between the
cludedthroughextended-MHDmodeling, areprevi- codes except for n =4 where there is a 27% rela-
φ
ouslyavailableinRefs.[27]and[31]andthusourde- tivedifference. Allothercaseshavearelativediffer-
scriptionofthisaspectisbrief. InSec.IIIthemodel ence of less than 8% with typical differences of 5%.
complexityisgraduallyincreasedfromtheideal-like The NIMROD convergence in terms of the maxi-
model through a resistive-MHD model with exper- mum polynomial order of the spectral elements is
imentally relevant profiles to a full extended-MHD shown in Fig. 2. Convergence is most challenging
3
Figure 1: [Color online] Growth rates for the ‘Meudas-
1’ benchmark. ELITE with Γ = 5/3 and Γ = 0 are
Figure 3: [Color online] Poloidal cross section of the ra-
comparedagainstresultsfromNIMRODwithΓ=5/3).
dial magnetic field component of the n = 11 peeling-
Associated NIMROD data available in Ref. [33]. φ
ballooning mode from the ‘Meudas-1’ benchmark case.
Figure 3 shows a poloidal cross section of the
magnetic (B ) eigenmode. The mode develops
R
an ‘interference-pattern’ structure near the X-point
when inboard and outboard finger-like structures
overlap. The finite-element-mesh nodes are super-
imposed atop the smallest-scale sub-figure. As es-
tablished by Fig. 2, this simulation is spatially and
temporally converged. The high resolution required
to resolve these high-q , diverted cases motivated
95
developmentofmemory-scalingimprovementsinthe
NIMROD code.
III. INCREASING THE MODEL
Figure 2: [Color online] Spectral convergence of the COMPLEXITY
NIMRODcodefortheideal-likeparameters. Themaxi-
mum polynomial degree (P) of the basis functions com-
Initial-value extended-MHD computations, while
posing the spectral elements in increased in each subse-
more computationally expensive than ideal-MHD
quent line plotted. Associated NIMROD data available
calculations with ELITE, are able to include addi-
in Ref. [33].
tional effects that are more representative of exper-
imental conditions:
at high wavenumbers where the resolution require- • Instead of a magneto-static vacuum outside
ments are most stringent (the poloidal mesh is com- the LCFS, NIMROD includes a cold-plasma
posed of 72×512 spectral elements). These cases region.
converge from the unstable side where the growth • Finite-dissipation effects such as resistivity,
rate decreases with enhanced resolution. Thus the viscosity and thermal conduction.
excellent agreement between NIMROD and ELITE • Density and temperature profiles that pre-
at high nφ in Fig. 1 may be spurious and indicate scribetheassociateddriftsandresistivitypro-
that slightly more resolution is required for nφ>25, file.
however, the shown growth rates are likely within • Parallel-closureeffectsthatrepresentlargedif-
5% of their converged values. Studying nearly ideal fusivities oriented along the magnetic field.
cases with extended MHD codes such as NIMROD • First-order FLR-closure and two-fluid effects.
ischallenginggiventhevanishinglysmalldissipation
operators,andconvergenceisachievedmorequickly Reference[31]investigatesthefirstthreeoftheseef-
withtheadditionalnon-idealtermsintheextended- fects with the Meudas-1 case. We further this work
MHD equations, as in the cases in Sec. III. by including the last two effects. The first three
4
effects are also re-investigated and different modifi-
cations to the growth rates are found. These differ-
ences are expected as slightly different profiles for
density and temperature are used in order to avoid
spuriouselectrostaticmodeswiththetwo-fluid,first-
order FLR model. Ref. [31] avoided these modes
by employing only the single-fluid resistive-MHD
model.
In calculations beyond the ideal model, the den-
sity profile is n (ψ ) = n (p(ψ )/p )0.3 where
e n e0 n 0
n =6×1019 m−3andψ isthenormalizedpoloidal
e0 n
flux. Herepisthepressureasprescribedbytheideal
gaslaw,p=n T +n T ,whereT isaspecies’tem-
e e i i α
perature. This choice does not match Ref. [31], but
effectively avoids spurious electrostatic modes with
Figure4: [Coloronline]GrowthratescomputedbyNIM-
thetwo-fluid,first-orderFLRmodel. Thesespurious RODonthe‘Meudas-1’casewithfourdifferentmodels:
modes are discussed in more detail in Sec. V. The ideal MHD (ideal), resistive MHD with reconstructed
pressure-profile gradient is specified by the Meudas- density,temperature,andSpitzer-resistivityprofiles(1f),
1 case and setting the edge temperature to 200eV the 1f model with parallel Braginskii closures (1f-Par),
leads to a core temperature of 5806eV. The edge the 1f-Par model with ion gyroviscosity, a full extended
temperature is modified independent of the MHD- MHD Ohm’s law (including the Hall, ∇pe, and elec-
stability (determined by p0) by the transformation tron inertia terms) and separate temperature evolution
equations with cross-heat fluxes (2f-Par-FLR). Associ-
p(ψ )→p(ψ )+p wherep isaconstant. Without
n n c c ated NIMROD data available in Ref. [33].
scrape-off-layer (SOL) profiles of density and tem-
perature, which are not included in standard recon-
structions, this choice of the edge temperature also
where W is the rate-of-strain tensor
specifies the temperature outside the LCFS. SOL
profiles are required to get both the correct profiles W=∇v+∇vT −(2/3)I∇·v, (4)
outsidetheLCFS,whichdeterminesthevacuumre-
sponse, and simultaneously set the correct resistiv- a full extended-MHD Ohm’s law,
ity at the mode resonant location when Spitzer re-
sistivity is used. The non-ideal viscous, conductive J×B ∇p m ∂J
E=−v×B+ − e +ηJ+ e , (5)
andparticlediffusivitiesaresettoasmallvalue,one n e n e n e2 ∂t
e e e
tenth of the resistivity at the LCFS.
Figure4plotsthegrowthratescomputedbyNIM- andseparatetemperatureevolutionswithcross-heat
ROD on the ‘Meudas-1’ case with four different fluxes,
models: ideal MHD (ideal), resistive MHD with
5p
reconstructed density, temperature, and Spitzer- q = α bˆ×∇T , (6)
resistivity profiles (1f), the single-fluid model with ×α 2qαB α
parallel Braginskii closures (1f-Par) where the par-
(2f-Par-FLR).Herev isthebulk-ionflow, ν isthe
allel viscosity is i ki
ion parallel diffusivity, χ is the species’ parallel
kα
(cid:18) (cid:19) diffusivity, q (e) is the species’ (electron) electric
1 α
Πki =miniνki bˆbˆ− 3I charge, J is the current density, bˆ is the magnetic
unit vector (bˆ =B/B), and I is the identity tensor.
(cid:16) (cid:17)
× 3bˆ·∇v ·bˆ−∇·v , (1) Aneffectiveioncharge(Z)ofthreeisassumed. Each
i i
modelbuildsonthepreviouslylistedandincludesall
and the parallel heat flux vector is prior terms. Electron viscosity is not included.
Whenadensityprofileisincluded,thenormalized
q =−n χ bˆbˆ·∇T , (2)
kα i kα α growth rate of the mode remains constant with the
single-fluidmodelwhentheAlfvéntimeiscomputed
and a two-fluid model with parallel closures that in-
at the mode resonant surface. Our calculations are
cludes ion gyroviscosity
normalized by the Alfvén time as computed with
Π = mipi [bˆ×W·(cid:16)I+3bˆbˆ(cid:17) the values at the magnetic axis. These values re-
×i 4ZeB mainconstantwhenadensityprofileisincludedsuch
−(cid:16)I+3bˆbˆ(cid:17)·W×bˆ], (3) that the normalized mode growth rate is effectively
increased.
5
The effect of the resistivity profile is more nu- field,withoutcurvatureterms(slabapproximation),
anced. Finite resistivity allows for reconnection and in the low-β limit (see, e.g. [40]). Use of
within the model by relaxing the frozen-flux con- the full gyroviscous operator is critical as in addi-
straint. This permits resistive-ballooning modes tion to effects that appear as diamagnetic drifts,
that grow faster than their ideal counterparts. Al- it also contains terms proportional to first-order
ternatively,theresistivedissipationcanstabilizethe FLRmagnetic-curvatureandgrad-Bdrifts[41]. Our
modes by acting as a dissipative term and modify- modelingcontainsnotjustthefirst-orderFLRdrifts
ing the response outside the LCFS from that of a from ion gyroviscosity, but also the drifts from the
vacuum to that of a plasma. Figure 5 of Ref. [31] cross-heat fluxes that enter the equations on the
showsthattheplasmaresponseisstabilizingrelative same order.
to a vacuum region. Overall, these effects stabilize Just as the large balance of the ideal-MHD forces
the ballooning modes at low-nφ and destabilize the can cause numerical difficulties with the ideal-MHD
modes at high-nφ as seen in Fig. 4 when comparing force operator, the partial gyroviscous cancella-
the single fluid and ideal normalized growth rates. tion leads to similar numerical pitfalls and verifi-
The effect of including small particle, viscous, and cation is important. The implementation of these
thermal diffusivities on the resistive-MHD mode is terms in NIMROD are verified in cylindrical and
small (not shown). slab magnetic geometries against analytic calcula-
Including the Braginskii parallel closures with tions for tearing [14, 40, 41] and ion-temperature-
large coefficients (χ = 109 × η /µ and χ = gradient[42](ITG)modes. Totesttheeffectofthese
ke 0 0 ki
ν = 108 ×η /µ ) has little effect on the growth terms against analytic theory in full tokamak mag-
ki 0 0
rates. When thermal conduction is large, the tem- neticgeometry,weperformaqualitativecomparison
peraturequicklyequilibratesalongfieldlinestopro- with analytic theory in Sec. IV. Full verification is
duce an isothermal response. Thus NIMROD com- precluded by the approximations made in the ana-
putations that show little effect from the parallel lytic theories.
closures are consistent with the ELITE results that Thefirst-orderFLRmodelisvalidforthetoroidal
find the shape of the mode growth-rate spectrum modespresentedalthoughthisisnotnecessarilythe
is largely not sensitive to the value of the ratio of caseforgeneraledgemodemodeling. Assumingk ’
specific heats, Γ, as seen in Fig. 1. When the re- (q/r+1/R)n andevaluatinglocalquantitiesatthe
φ
sponse to the compressible motion associated with peak of the mode eigenfunction (ψ =0.969) on the
n
the sound wave is eliminated (Γ = 0 in the ELITE outboa√rd midplane, the normalized ion-gyroradius
computations), the growth rate is slightly enhanced (ρ = Γm T /ZeB) as a function of toroidal mode
i i i
relative to the adiabatic limit. If the growth rates iskρ ’0.0055n . Ifthefirst-orderFLRassumption
i φ
aremodifiedbeyondthissmalleffect,itwouldbean isviolated,more-complexfull-ion-orbitkineticmod-
indication that changes in the energy equation can eling is required to simulate large-fluctuation ELM
impacttheMHDresponse, andthusthatclosureef- dynamics. For this case, computations produce ap-
fects are significant. Although the short-mean-free proximately the same results with and without the
path Braginskii-like closure [34, 35] is used in our cross-heat fluxes. This is consistent with drift an-
NIMROD compuatations, we note that other clo- alytics [40] that shows that the cross-heat flux is
sures, such as long-mean-free path [36] or the gen- only significant at very low values of plasma β (here
eral approach of solving drift-kinetic equation [37] the β = 2µ p/B2 is 0.1%). Perhaps coincidentally
0
for a closure can be applied. However, for this case for this case, the effects of the resistive destabiliza-
where the parallel-closure effects are negligible, dif- tion and drift stabilization largely counteract one
ferent parallel closures are unlikely to significantly another. The drift stabilized growth rates are not
modify the results. less than the ideal calculation until n (cid:38)24.
φ
With the full two-fluid, FLR model, there is a
stabilizing effect on the intermediate and high-n
φ
modes as expected from analytic treatments [30].
IV. COMPARISON TO DRIFT ANALYTICS
A common reduced-MHD approximation for first-
order FLR closures is to assume that the ion-
gyroviscous force and the advection by the diamag- To gain insight into our two-fluid computations,
netic drift exactly cancel (this is colloquially known we compare with analytic descriptions of drift sta-
asthegyroviscouscancellation,see,e.g. [38]). How- bilization of the ballooning mode. The dispersion
ever,asdiscussedbyRamosinRef.[39],“thesecan- relation for drift stabilizations of the ideal, incom-
cellations are only partial and not very useful in pressible ballooning mode is
practice for general magnetic geometries...”. This
cancellationisonlyvalidwithalarge,uniformguide ω(ω−ω )=−γ2 (7)
∗i I
6
Figure5: [Coloronline]GrowthratescomputedbyNIM- Figure6: [Coloronline]GrowthratescomputedbyNIM-
ROD with the ideal, single-fluid (1f) and full FLR (2f- ROD with the ideal, single-fluid (1f) and full FLR (2f-
Par-FLR) models in comparison to the analytic dis- Par-FLR) models in comparison to the analytic disper-
persion relation of the FLR-stabilized, ideal ballooning sionrelationofRef.[30],Eqn.(43)withthesameparam-
modewithfullandreducedω values. AssociatedNIM- etersastheNIMRODcomputationsexceptreducedω
∗α ∗α
ROD data available in Ref. [33]. values and a reduced resistivity (0.35 of the NIMROD
value). Associated NIMRODdata availablein Ref.[33].
where γ is the ideal growth rate in the absence of
I
drifteffectsandω isthecomplexfrequency(ω =iγ) parewiththistheoreticaltreatment,parametersare
[43]. The diamagnetic drift velocity is typically computed on the outboard midplane at the radial
defined as the perpendicular-to-the-magnetic-field locationwherethemodestructurepeaks, consistent
flow contribution from the pressure gradient within with the previous ω∗i calculation. This yields lo-
the context of a two-fluid, or generalized Ohm’s cal values of the safety factor, q = 7.34, normalized
law. For a tokamak, neoclassical effects damp the magnetic shear, s = 8.8 and normalized pressure
poloidal diamagnetic flow contribution. Thus in our gradient, α=28.8, wherethedefinitionsofRef.[30]
case the remaining toroidal diamagnetic flow deter- areused. The∆0 stabilityparameterisdetermined
B
mines the diamagnetic frequency, ω = k ·v = byusingtheidealgrowthrate(Γ=5/3)andsolving
∗i ∗i
n /(n q )∂p /∂ψ where k is the mode wavenumber. Eqn. 35 of Ref. [30]. ∆0 varies with toroidal mode
φ i i i B
Figure5showsthecomparisontothisdispersionre- number but falls in the range of −17.4 to −12.1.
lation with full and reduced diamagnetic drift val- Comparison of HRP Eqn. (43) with ω = 0 and
∗
ues. The value of ω is computed at the radial lo- the NIMROD resistive-MHD results differ by a fac-
∗i
cation where the mode structure peaks. The nor- tor approximately 2.5 at large toroidal mode num-
malized frequency from the diamagnetic drift veloc- berswheretheanalyticcalculationresultsingreater
ity, ω τ , varies linearly with toroidal mode num- growth rates and there is no stabilization at low
∗i A
beras0.00483n . This‘local’approximation, which n . As the analytic calculation is a local calcula-
φ φ
neglects the radial variation of the ω profile, is tion, it does not capture the stabilizing influence on
∗i
known to over-estimate the effect of drift stabiliza- the global mode of finite resistivity, in particular it
tion[44,45]. Wefindthatareductionofthevalueof misses the effect of modifying the response outside
ω (approximatelybyafactorof2-4)givesthebest the LCFS from that of a vacuum to that of a resis-
∗i
agreementbetweenEqn.(7)andthetwo-fluidNIM- tive plasma. Our interest is largely to compare and
ROD calculation. A similar comparison of ELITE contrast the predicted drift stabilization and thus
tothetwo-fluidBOUT++modelonadifferentcase we reduce the resistivity to 35% of the value used
finds the drift stabilization is over-predicted [45] by intheNIMRODcomputationswhencalculatingthe
a comparable factor. analytic values of HRP Eqn. (43). This produces
One of the limitations of Eqn. (7) is that it does growth rates with ω∗ = 0 that roughly match the
not include the effect of resistive destabilization on NIMROD resistive-MHD computations at high nφ
the high-n modes, as is included in our two-fluid as shown in Fig. 6.
φ
computations. The analytic treatment of Ref. [30] More significantly, Fig. 6 compares the results of
provides a dispersion relation (Eqn. (43) of the ref- our linear NIMROD computations with the drift-
erence; referred to as HRP Eqn. (43) in this discus- stabilized growth rates computed from Eqn. (43) of
sion) that includes the effect of resistive destabiliza- Ref. [30]. Again, the analytics of HRP Eqn. (43)
tion in addition to drift effects. In order to com- over-estimate the effect of drift stabilization where
7
of drift-stabilization is stronger at low densities (ρ
√ i
scales as 1/ n at constant β) as expected. The
i
cases with n = 3×1019 m−3 are not plotted for
e
n > 18 as the most unstable mode is no longer
φ
thepeeling-ballooningmodebutratheradominatly
electro-static mode related to the ITG as studied in
thecontextNIMRODtwo-fluidadvanceinRef.[42].
ThisregionisexcludedasexaminationofITGmode
dynamics is outside the scope of this work.
VI. DISCUSSION AND CONCLUSIONS
Our comparisons to the ideal PBM growth rates
Figure7: [Coloronline]GrowthratescomputedbyNIM- computed by the ELITE code and analytic descrip-
RODonthe‘Meudas-1’casewiththe1fmodelwithpar- tionsoftheresistive,two-fluidPBMprovideanother
allel Braginskii closures (1f-Par; single-fluid limit), and
verificationtestoftheNIMRODcodeasameansto
three different densities with the two-fluid, FLR model
study the tokamak edge. Quantitative agreement
with parallel Braginskii closures (2f-Par-FLR): the core
with ELITE is achieved for these cases, extending
densitiesaren =3×1019,6×1019 and1.2×1020 m−3.
e prior work [27] that did not include diverted mag-
Associated NIMROD data available in Ref. [33].
netic topology.
Two-fluid and FLR-drift effects can both produce
stabilizing (drift stabilization through the sheared
the reduction of ω that leads HRP Eqn. (43)
∗
motion of the ion and electron fluid responses)
to qualitatively agree with the two-fluid NIMROD
and destabilizing (species decoupling and enhanced
computations is between 2 to 4 (dependent on n ).
φ
growth through mediation by the more mobile elec-
There is little effect on the low-n modes and the
φ
tron fluid) effects (see, e.g., [40, 41]). At intermedi-
iondrift-waveresonancepredictedinRef.[30]isnot
ate and large mode numbers, the inclusion of finite-
observed.
resistivityisdestabilizingwhiletheFLR-drifteffects
are stabilizing. These effects are largely offsetting,
resulting in a toroidal-mode growth-rate spectrum
V. FLR PARAMETER SCAN with the full extended-MHD model that resembles
the ideal calculation. This result is not general but
In order to ascertain the role of two-fluid effects rather case dependent. Consistent with ELITE cal-
further, we vary the density and temperature while culations that vary the ratio of specific heats from
scaling parameters to maintain constant β and S the adiabatic through the isothermal limit, we find
(other dissipation parameters are scaled relative to the addition of anisotropic thermal conduction pro-
resistivity). In experimental discharges, the zeroth- duces only a small modification of the growth rate.
order effect on the edge of modifying the density Relative to the NIMROD calculations, analytics
is modify the plasma collisionality and thus the [30] over-predict both the effects of resistive desta-
bootstrap current. This causes a transition from bilization and FLR-drift stabilization for the case
ballooning-likemodesathighdensitytopeeling-like studied. Qualitativeagreementwiththeanalyticde-
modes at low density. Our computations use a fixed scriptions of the peeling-ballooning mode is found
equilibriumcurrentandthusarenotsensitivetothis where high-n modes are susceptible to resistive
φ
effect. Instead, we isolate the effect of the modifi- destabilizationanddriftstabilizationinboththean-
cation of ρ (and the associated ion skin depth, d ) alyticsandNIMRODcomputations. Oneofthelim-
i i
within the context of our two-fluid, first-order FLR its of the analytic derivation is a ‘local’ approxima-
model. This contrasts with Ref. [46] that examines tion which leads to the over prediction of FLR-drift
these effects in concert and Ref. [47] that considers stabilization [44]. Given the highly complex nature
onlytheeffectofmodificationstothecurrentprofile. of both the configuration of the tokamak edge and
Figure 7 shows the linear growth rate from NIM- the two-fluid model equations, quantitative agree-
ROD computations for resistive-MHD and three ment is not expected as the local approximation ne-
two-fluid cases with varying densities (the core den- glects the global configuration and profile effects.
sity is varied from 3×1019 m−3 to 1.2×1020 m−3 Regarding our investigation into the impact of
where 6×1019 m−3 is the value used in the prior the size of the ion gyroradius within the context
computations discussed in the paper). The effect of NIMROD’s full first-order FLR, two-fluid model,
8
we find significant changes on the growth rates of N. Aiba for providing the Meudas-1 equilibrium.
high-n modes at varying densities but constant β. This material is based on work supported by US
φ
However, the most unstable mode number is largely Department of Energy, Office of Science, Office of
unchanged. Thus we conclude the zeroth-order ef- Fusion Energy Sciences under award numbers DE-
fectofmodificationstothebootstrapcurrentlargely FC02-06ER54875 and DE-FC02-08ER54972. This
dominates the determination of the most unstable research used resources of the National Energy Re-
mode during collisionality scans while lower densi- searchScientificComputingCenter,aDOEOfficeof
ties (large ρ ) increase the magnitude of the drift Science User Facility supported by the Office of Sci-
i
stabilization on the high-n modes consistent with ence of the U.S. Department of Energy under Con-
φ
the discussion of Ref. [48]. tract No. DE-AC02-05CH11231.
Acknowledgments
We thank Carl Sovinec, Chris Hegna and Nate
Ferraro for discussions involving this paper and
[1] A. W. Leonard, N. Asakura, J. A. Boedo, M. Be- Physics of Plasmas 21, 112302 (2014).
coulet, G. F. Counsell, T. Eich, W. Fundamenski, [13] C. Sovinec, A. Glasser, T. Gianakon, D. Barnes,
A. Herrmann, L. D. Horton, Y. Kamada, A. Kirk, R. Nebel, S. Kruger, D. Schnack, S. Plimpton,
B. Kurzan, A. Loarte, J. Neuhauser, I. Nunes, A.Tarditi, andM.Chu,JournalofComputational
N. Oyama, R. A. Pitts, G. Saibene, C. Silva, P. B. Physics 195, 355 (2004).
Snyder, H.Urano, M.R.Wade, H.R.Wilson, and [14] C. Sovinec and J. King, Journal of Computational
f. t. P. a. E. P. I. Group, Plasma Physics and Con- Physics 229, 5803 (2010).
trolled Fusion 48, A149 (2006). [15] J.Callen,C.Hegna,B.Rice,E.Strait, andA.Turn-
[2] A. Pankin, G. Bateman, D. Brennan, A. Kritz, bull, Phys. Plasmas 6, 2963 (1999).
S.Kruger,P.Snyder, andC.Sovinec,Plasma.Phys. [16] S.Kruger,D.Schnack, andC.Sovinec,Phys.Plas-
Control. Fusion 49, S63 (2007). mas 12, 056113 (2005).
[3] L. E. Sugiyama and H. R. Strauss, Physics of Plas- [17] S.C.Jardin,J.Breslau, andN.Ferraro,Journalof
mas 17, 062505 (2010). Computational Physics 226, 2146 (2007).
[4] X.Q.Xu,B.Dudson,P.B.Snyder,M.V.Umansky, [18] B. Dudson, M. Umansky, X. Xu, P. Snyder, and
andH.Wilson,PhysicalReviewLetters105(2010), H.Wilson,ComputerPhysicsCommunications180,
10.1103/PhysRevLett.105.175005. 1467 (2009).
[5] X. Q. Xu, P. W. Xi, A. Dimits, I. Joseph, M. V. [19] G. Huysmans, S. Pamela, E. Van der Plas, and
Umansky,T.Y.Xia,B.Gui,S.S.Kim,G.Y.Park, P. Ramet, Plasma Physics and Controlled Fusion
T.Rhee,H.Jhang,P.H.Diamond,B.Dudson, and 51, 124012 (2009).
P. B. Snyder, Phys. Plasmas 20, 056113 (2013). [20] D. D. Schnack, D. C. Barnes, D. P. Brennan, C. C.
[6] P. Xi, X. Xu, T. Xia, W. Nevins, and S. Kim, Nu- Hegna, E. D. Held, C. C. Kim, S. E. Kruger, A. Y.
clear Fusion 53, 113020 (2013). Pankin, andC.R.Sovinec,PhysicsofPlasmas13,
[7] G. Huijsmans and A. Loarte, Nuclear Fusion 53, 058103 (2006).
123023 (2013). [21] T. G. Jenkins, S. E. Kruger, C. C. Hegna, D. D.
[8] S. J. P. Pamela, G. T. A. Huijsmans, A. Kirk, I. T. Schnack, and C. R. Sovinec, Physics of Plasmas
Chapman,J.R.Harrison,R.Scannell,A.J.Thorn- 17, 012502 (2010).
ton, M. Becoulet, F. Orain, and the MAST Team, [22] R. Gruber and M. Rappaz, Finite element methods
Plasma Physics and Controlled Fusion 55, 095001 inlinearidealmagnetohydrodynamics(SpringerSci-
(2013). ence & Business Media, 2012).
[9] J.W.Connor,R.J.Hastie,H.R.Wilson, andR.L. [23] L. Degtyarev and S. Y. Medvedev, Computer
Miller, Physics of Plasmas 5, 2687 (1998). Physics Communications 43, 29 (1986).
[10] P. B. Snyder, H. R. Wilson, J. R. Ferron, L. L. [24] A.BondesonandG.Fu,Computerphysicscommu-
Lao, A. W. Leonard, T. H. Osborne, A. D. Turn- nications 66, 167 (1991).
bull, D. Mossessian, M. Murakami, and X. Q. Xu, [25] H. Lütjens and J. Luciani, Computer physics com-
Physics of Plasmas (1994-present) 9, 2037 (2002). munications 95, 47 (1996).
[11] P.B.Snyder,H.R.Wilson,J.R.Ferron,L.L.Lao, [26] C. Sovinec, submitted to Journal of Computational
A.W.Leonard,D.Mossessian,M.Murakami,T.H. Physics (2015).
Osborne, A. D. Turnbull, and X. Q. Xu, Nuclear [27] B.J.Burke,S.E.Kruger,C.C.Hegna,P.Zhu,P.B.
fusion 44, 320 (2004). Snyder, C. R. Sovinec, and E. C. Howell, Physics
[12] B.Gui,X.Q.Xu,J.R.Myra, andD.A.D’Ippolito, of Plasmas 17, 032103 (2010).
9
[28] H. R. Wilson, P. B. Snyder, G. T. A. Huysmans, [36] J. J. Ramos, Physics of Plasmas (1994-present) 18,
andR.L.Miller,PhysicsofPlasmas(1994-present) 102506 (2011).
9, 1277 (2002). [37] E. D. Held, S.E. Kruger, J.-Y. Ji, E.A. Belli, and
[29] P. Snyder, N. Aiba, M. Beurskens, R. Groebner, B.C.Lyons,PhysicsofPlasmas22,032511(2015).
L. Horton, A. Hubbard, J. Hughes, G. Huysmans, [38] B. Coppi, Physics of Fluids 7, 1501 (1964).
Y. Kamada, A. Kirk, C. Konz, A. Leonard, J. Lön- [39] J. J. Ramos, Physics of plasmas 14, 2506 (2007).
nroth,C.Maggi,R.Maingi,T.Osborne,N.Oyama, [40] J.R.KingandS.E.Kruger,PhysicsofPlasmas21,
A. Pankin, S. Saarelma, G. Saibene, J. Terry, 102113 (2014).
H. Urano, and H. Wilson, Nuclear Fusion 49, [41] J. R. King, C. R. Sovinec, and V. V. Mirnov,
085035 (2009). Physics of Plasmas 18, 042303 (2011).
[30] R. J. Hastie, J. J. Ramos, and F. Porcelli, Physics [42] D. D. Schnack, J. Cheng, D. C. Barnes, and S. E.
of Plasmas 10, 4405 (2003). Parker, Physics of Plasmas 20, 062106 (2013).
[31] N. M. Ferraro, S. C. Jardin, and P. B. Snyder, [43] W. Tang, R. Dewar, and J. Manickam, Nuclear
Physics of Plasmas 17, 102508 (2010). Fusion 22, 1079 (1982).
[32] N.Aiba,S.Tokuda,T.Fujita,T.Ozeki,M.S.Chu, [44] R. J. Hastie, P. J. Catto, and J. J. Ramos, Phys.
P.B.Snyder, andH.R.Wilson,PlasmaFusionRes. Plasmas 7, 4561 (2000).
2, 010 (2007). [45] P. Snyder, R. Groebner, J. Hughes, T. Osborne,
[33] J.R.King,“Nimrodgrowthrates(1/s)byfiniteele- M. Beurskens, A. Leonard, H. Wilson, and X. Xu,
mentpoly_degree(columns)andcasefrom“theim- Nuclear Fusion 51, 103016 (2011).
pact of collisionality, flr and parallel closure effects [46] X.Q.Xu,J.F.Ma, andG.Q.Li,PhysicsofPlas-
oninstabilitiesinthetomakakpedestal: Numerical mas 21, 120704 (2014).
studies with the nimrod code”,” (2016). [47] P.Zhu,C.C.Hegna, andC.R.Sovinec,Physicsof
[34] S. I. Braginskii, Transport Properties in a Plasma Plasmas 19, 032503 (2012).
in Review of Plasma Physics, edited by M. A. [48] P. Snyder, K. Burrell, H. Wilson, M. Chu, M. Fen-
Leontovich,Vol.1(ConsultantsBureau,NewYork, stermacher, A. Leonard, R. Moyer, T. Osborne,
1965). M. Umansky, W. West, and X. Xu, Nucl. Fusion
[35] P. J. Catto and A. N. Simakov, Physics of Plasmas 47, 961 (2007).
11, 90 (2004).