Magnetohydrodynamic calculations with a nonmonotonic q profile and equilibrium , sheared toroidal flow

The linear and nonlinear stability of a nonmonotonic q profile is examined using a reduced set of magnetohydrodynamic ~MHD! equations with an equilibrium, sheared toroidal flow. The reversed shear profile is shown to be unstable to a rich variety of resistive MHD modes including pressure-driven instabilities and tearing instabilities possessing a tearing/interchange character at low Lundquist number,S, and taking on a double/triple tearing structure at high S. Linear calculations show that the destabilizing effect of toroidal velocity shear on tearing modes is enhanced at finite pressure. In addition, this velocity shear decreases the stabilizing effect of finite pressure seen previously for tearing modes at high S. Nonlinear calculations show the generation of a large, m51, n50, Reynolds-stress-driven poloidal flow in the absence of significant flow damping. Calculations in which the poloidal flow was heavily damped show that sub-Alfve ́nic, sheared toroidal flows have a minimal effect on weakly coupled, localized instabilities. © 1999 American Institute of Physics. @S1070-664X~99!00603-5#


I. INTRODUCTION
The enhanced core confinement and high neutron rate of reversed shear plasmas present a promising operating regime for future tokamaks. 1Evidence suggests, however, that negative shear plasmas are susceptible to a variety of magnetohydrodynamic ͑MHD͒ activity including resistive interchange and tearing modes.These instabilities degrade confinement and may trigger violent plasma disruptions. 2For this reason, it is important to understand the MHD stability properties of nonmonotonic q profiles in order to better take advantage of this regime of enhanced confinement.
Reversed shear plasmas are characterized by large, sheared toroidal flows resulting from early neutral beam injection. 1This paper examines the influence of toroidal flow shear on resistive tearing mode growth rates in a finite-␤, reversed-shear plasma using a reduced set of MHD equations with a sub-Alfve ´nic, equilibrium toroidal flow.Previous studies regarding the linear stability of tearing modes have addressed separately the issues of finite-␤ 3,4 and shear flow effects 5,6 in toroidal plasmas.The linear results presented here show the importance of considering both effects simultaneously.
Reversed shear plasmas are also characterized by large, sheared poloidal flows that may provide a barrier to crossfield transport near the reversed shear region.In this work, the extremes of zero and infinite poloidal flow damping are assumed in order to study the nonlinear stability of negative shear profiles.In the absence of damping, it is shown that small-scale MHD activity can generate large, Reynoldsstress-driven, sheared poloidal flows.Because sheared flow may also affect the nonlinear coupling of localized instabilities, nonlinear calculations are performed in the presence of infinite poloidal flow damping to determine the extent to which sheared toroidal flow provides a stabilizing influence by decoupling modes of a common helicity.Because the toroidal flow in this work is neither evolved nor included self-consistently in the equilibrium, it is ordered at onehundredth the Alfve ´n speed so as to minimize its presence as a free energy source and its effect on the equilibrium.A reduced version of the resistive MHD stability code FAR, 7,8 with the effects of toroidal flow, is used to study the linear and nonlinear stability of the q profile shown in Fig. 1.
This paper is organized as follows.In Sec.II, we discuss how an equilibrium, sheared toroidal flow modifies the standard reduced MHD model.In Sec.III we examine the combined effects of sheared flow and finite-␤ on tearing mode growth rates.In Sec.IV we discuss the nonlinear stability of nonmonotonic q profiles in the presence of sheared toroidal flow and in Sec.V we summarize.

II. REDUCED MHD EQUATIONS WITH EQUILIBRIUM, SHEARED TOROIDAL FLOW
In this section we discuss reduced MHD equations 9 with the effects of an equilibrium, sheared toroidal flow.Reduced equations assume a flute mode ordering whereby the main magnetic and fluid perturbations lie in a plane perpendicular to the dominant equilibrium magnetic field.Although more refined models allow for arbitrary aspect ratio, 10 the equations used in this work depend heavily upon the standard, large aspect ratio assumption as well, ⑀ϭa/R 0 Ӷ1, where a and R 0 is the minor and major radius, respectively.The reduced MHD model incorporates the pertinent physics of full MHD in the following set of three scalar equations describ-ing the evolution of the poloidal flux, , the toroidal component of the vorticity, U, and the pressure, P: The equation for the poloidal flux, , is the toroidal component of the resistive Ohm's law.The vorticity equation is the result of operating on the momentum equation with •ٌϫ1/R and the pressure equation is simply the incompressible equation of state.To lowest order in , the parallel direction is in the toroidal direction denoted by the unit vector, ˆ, with the associated toroidal magnetic field curvature to lowest order, ٌ(1/R).
Additional definitions include

͑7͒
Here is the velocity stream function, J /R is the toroidal current density, and an inductive electric field, E /R, is added to Ohm's law to provide a means for calculating the equilibrium resistivity profile, .Perpendicular viscosity and pressure diffusivity appear as dissipative terms in the vorticity and pressure equations with the associated coefficients, and D P , respectively.Although these terms provide negligible physical dissipation, they represent a small but important amount of numerical dissipation in nonlinear calculations.The ratio of the resistive and poloidal Alfve ´n times, Sϭ R / Hp , is referred to as the Lundquist number, where R ϭa 2 0 / and Hp ϭR 0 ( 0 ) 1/2 /B 0 .The aspect ratio is again ϭa/R 0 , and the constant, ␤, is the ratio of the kinetic and magnetic pressures at the magnetic axis.Time and length scales have been normalized to R and a, respectively.
Before continuing, a few comments are in order concerning the nature of these reduced equations.First, in addition to the dominant EϫB flow velocity given by v Ќ , an equilibrium toroidal flow, v , has been retained in the advective terms of the vorticity and pressure equations.Although v is normally ordered out of reduced equations, it is shown later that, depending on the velocity shear, small levels of toroidal flow can have an important effect on resistive tearing modes.Note that the velocity appearing in the equation for is strictly perpendicular.
Second, we address the issue of the following energy integral for Eqs.͑1͒-͑3͒, where the equilibrium toroidal velocity component appears as a possible source of free energy on the right-hand side,

͑8͒
Here the integral extends over the plasma volume and and are assumed to vanish at the plasma boundary.The terms on the left-hand side represent the kinetic, magnetic, and stored energy of the plasma, respectively.The dissipative resistive Joule heating, diffusivity and viscosity terms on the right-hand side provide a small sink of free energy that is beneficial for nonlinear calculations.The advective term from the vorticity equation, however, may be positive or negative definite.This demands that the magnitude of the equilibrium flow be small enough to prevent it from entering as an energy source in nonlinear calculations, but large enough to have an effect on mode coupling.For this reason, the toroidal flow is ordered as v /v Hp р0.01, where v Hp ϭB 0 /( 0 ) 1/2 is the poloidal Alfve ´n velocity.Nonlinear calculations confirm that this ordering is sufficient to prevent the toroidal flow from entering as a source of free energy.In addition, this ordering is sufficient to prevent toroidal flow, which was not included self-consistently in the equilibrium, from competing with the pressure gradient in the equilibrium momentum equation.This effect appeared in the distortion of linear eigenmodes whenever the magnitude of the toroidal flow was too large.
One final comment concerning this reduced MHD model is that the toroidal momentum equation is not evolved.This ignores the coupling between the toroidal and poloidal flows.FIG. 1.This q profile is linearly unstable to a resistive, pressure-driven mode in the inner region of positive magnetic shear and a resistive tearing mode residing mainly in the outer region of positive shear.Note the pressure gradient, dp/dr, is steepest in the inner region of positive shear.
For linear calculations of resistive tearing modes with submillisecond growth times, it is safe to assume that the toroidal flow evolution, which occurs on a millisecond to second timescale, decouples from Eqs. ͑1͒-͑3͒.In addition, we show later that at the small levels of toroidal flow and flow shear used in nonlinear calculations, this coupling of the flows is unimportant.
The linear and nonlinear reduced versions of the FAR 7,8 stability code solve Eqs.͑1͒-͑3͒ using finite differences in the radial coordinate and Fourier series expansion in the poloidal and toroidal angles.The linear terms are handled implicitly by inverting the eigenvalue matrix, whereas nonlinear terms are treated explicitly.One advantage of the FAR code is its ability to run as an initial value code and as an eigenvalue solver.This feature is especially useful when studying the linear stability of nonmonotonic q profiles.Judicious choice of the convergence parameter ͑timestep͒ permits rapid convergence to different linearly unstable eigenmodes with different growth rates.

III. LINEAR RESULTS
In this section we examine the linear stability of the nonmonotonic q profile shown in Fig. 1.In order to produce a resistive, pressure-driven instability in the inner region of positive shear (rр0.4) the pressure profile was tailored to peak around rϭ0.30.Here r is the minor radius coordinate normalized to a. Figure 2 shows the typical eigenmode structure of the velocity stream function, , centered at the qϭ5/2 surface for Sϭ10 5 and ␤ϭ1.5%.This resistive, pressure-driven mode is unstable at the 2/1,7/3,12/5,..., helicities as well, with growth rates, ␥ Hp ϭ␥• Hp , increasing roughly linearly with toroidal mode number, n.It is important to note that this mode is not representative of the stability properties of reversed shear profiles, but rather was created to study nonlinear coupling of localized modes in the presence of sheared toroidal flow.For this reason we do not investigate its linear stability properties in depth.
By decreasing the convergence parameter, we locate a slower-growing, resistive tearing mode in the outer region of positive shear rу0.75 ͑see Fig. 3͒.For the 5/2 helicity, coupling to the middle and inner 5/2 surfaces is evident.At Sϭ10 5 and ␤ϭ1.5%, this mode is resistive tearing/ interchange in character in that it does not display a definite parity at the resonant surfaces, i.e., both and remain finite.Figure 4 shows, however, that at Sϭ10 8 and ␤ϭ0.0%, the mode takes on a ''triple tearing'' structure.Here we have chosen the more virulent 7/3 helicity to emphasize the strong coupling between surfaces of commensurate helicity in nonmonotonic q profiles.The 7/3 surface is more unstable than the 5/2 surface because of its proximity to the region of zero shear.Whereas Fig. 5 shows that at ␤ϭ1.5% this coupling is reduced by the stabilizing effect finite pressure has on tearing modes at high S, the triple tearing structure reappears in Fig. 6 when a toroidal flow profile like the one shown in Fig. 7 is imposed on the equilibrium with the shear parameter, s FIG. 2. Resistive, pressure-driven -eigenmode at the inner 5/2 surface for Sϭ10 5 and ␤ϭ1.5%.ϭ10.It is shown later that this result follows from the destabilizing effect that local flow shear has on the tearing mode at the outer surface.Linear calculations were carried out to study the combined effects of toroidal flow shear and finite pressure on the 5/2 harmonic of this resistive tearing mode.
Figure 8 shows the dominant tearing mode character at high S, with ␥ Hp falling off according to the familiar asymptotic scaling, S Ϫ3/5 , for ␤ϭ0.0% and zero toroidal flow. 12Grid-packing around the outer 5/2 surface has preserved the scaling even up to Sϭ10 10 .As seen previously, at ␤ϭ1.5%, ␥ Hp falls off faster than S Ϫ3/5 due to the stronger stabilizing effect of finite pressure at high S. 3,4 The falloff of ␥ Hp at lower S is due to the fact that the assumption, S →ϱ, used to derive the S 3/5 scaling is no longer valid.
Next, we examine the effect of toroidal flow shear on linear growth rates using the following equilibrium, toroidal velocity profile, v ϭv 0 tanh"s͑rϪr s ͒…, ͑9͒ where v 0 is the maximum flow and the parameter, s, controls the toroidal velocity shear at the rational surface, r s .Figure 7 shows the velocity profile for different s centered at the outer 5/2 surface, r s ϭ0.815.The effect of flow shear on ␥ Hp at different ␤ and S is shown in Figs. 9 and 10.For S ϭ10 8 , a small amount of flow shear destabilizes the tearing mode, whereas a larger amount returns ␥ Hp to the values without the flow in the case of ␤ϭ1.5% and cuts growth rates approximately in half in the case of ␤ϭ0% ͑see Fig. 10͒.Although the resistive layer is wider at Sϭ10 5 , the stabilizing and destabilizing effects of the toroidal flow are weaker for comparable amounts of flow shear ͑see Fig. 9͒.Because this mode takes on an interchange character at lower S, this is consistent with previous results showing little dependence between flow shear and growth rates for resistive interchange modes. 2 Although it has been demonstrated that poloidal flows on the order of the Alfve ´n speed can completely stabilize tearing modes, 6 we have chosen the ordering, v /v Hp р0.01, because it is more representative of the toroidal flows generated in experimental plasmas. 1 In addition to showing that modest levels of toroidal flow and flow shear (sр10) have a largely deleterious effect on ␥ Hp , Figs. 9 and 10 also display the relative effect of sheared flow at finite ␤.Finite-pressure stabilization is spoiled by flow shear as ␥ Hp quickly becomes comparable around sϭ5 for Sϭ10 8 and sϭ15 for Sϭ10 5 .This fact is reiterated in Figs.11 and 12, where ␥ Hp is plotted versus ␤ for sϭ0 and sϭ10.At Sϭ10 5 ͑see Fig. 11͒, toroidal flow shear reduces the stabilizing effect of increasing pressure, and at Sϭ10 8 ͑see Fig. 12͒, it nearly eliminates the stabilization entirely.In addition, flow shear destabilizes the tearing mode over a larger range in ␤.Note that these effects do not depend on the sign of the flow shear.These results give insight into the importance of the advective terms in the linear versions of Eqs.͑1͒-͑3͒.At low shear, coupling to the pressure equation via the curvature term in the vorticity equation dominates.This coupling is absent in the case of zero ␤.Thus the curvature term in the vorticity equation that is responsible for the finite-␤ stabilization effect is also responsible for the stronger destabilizing effect due to the advective term in the pressure equation.In other words, the destabilizing effect due to the advective term in the vorticity equation, included in zero-pressure models, takes a backseat to the pressure term at low flow shear.At larger values of shear the finite-and zero-␤ growth rates coincide as the advective term in the vorticity equation takes over.
The linear results presented here show qualitative agreement with previous models that treated flow shear and finite-␤ effects separately for single tearing modes.More importantly, however, they predict that for modest values of flow and flow shear, the coupling to the advective term in the pressure equation reduces the stabilizing effect of finite pressure and enhances the destabilizing effect of flow shear.Because of the double and triple tearing structure evident at Sϭ10 8 , it is important to distinguish the present results with previous models of double tearing modes that considered differential rotation between resonant surfaces but no local flow shear. 13In accordance with these results that showed differential flow to be strictly stabilizing, imposing a differential rotation between the middle and outer resonant surfaces revealed a strong stabilizing effect in the absence of local flow shear at the outer surface.Nonetheless, for the levels of flow, local flow shear, and differential flow considered in this paper, it is the local flow shear that dominates, i.e., the layer physics is more important than coupling between modes of commensurate helicity.In this regime, our results agree with previous studies showing that modest levels of flow and flow shear are strictly destabilizing.

IV. NONLINEAR RESULTS
In Sec.III we showed that the q-profile shown in Fig. 1 is linearly unstable to a resistive, pressure-driven instability in the inner region of positive shear and a resistive tearing mode in the outer region of positive shear with some coupling to the remaining resonant surfaces at Sϭ10 5 .The presence of linearly unstable modes at different rational surfaces permits the study of nonlinear coupling of localized instabilities for nonmonotonic q profiles.A typical set of 120 poloidal ͑m͒ and toroidal ͑n͒ Fourier components evolved in the nonlinear calculations is shown in Fig. 13, along with the q min ͑minimum q͒ and q max ͑maximum q͒ boundary lines.It is clear from Fig. 13 that the modes cover most of resonant space.For nonlinear calculations, a toroidal flow profile that varied linearly with radius from a maximum on axis to zero at the plasma boundary was imposed on the equilibrium.With this simple profile we investigate how sheared toroidal flow decouples resonant surfaces of the same helicity.
The first nonlinear calculation was performed at S ϭ10 5 , ␤ϭ2.0%, and an on-axis equilibrium flow, v /v Hp ϭ0.005.Here the flow shear at the outer resonant surfaces corresponds to sϭ0.5 ͑see Fig. 7͒. Figure 14 shows the evolution of the root mean square of the velocity stream function, rms , versus t R , where for nonlinear calculations, t R ϭt/ R .The Fourier coefficients have been squared and summed over an annulus of width 0.05 centered at rϭ0.25 and rϭ0.75.The faster-growing, pressure-driven instability grows and saturates nonlinearly around t R ϭ0.005.Although the tearing mode at rϭ0.75 grows linearly until t R ϭ0.004, it is soon overcome by the activity at the inner region of positive shear where the sharp increase in rms is caused by the pressure-driven mode suddenly being included in the volume average.Keep in mind that although the tearing mode is FIG.13.Plot of poloidal, m, versus toroidal, n, harmonics used in nonlinear calculations.The number of evolved modes ͑120͒ covers most of the resonant space between q min ϭ1.229 and q max ϭ3.755.FIG.14. Plot of log rms versus t R at ␤ϭ2.0% and Sϭ10 5 .The pressuredriven mode in the inner region of positive shear grows linearly and saturates (t R ϭ0.005).Reynolds stress generates a large, mϭ1, poloidal flow around t R ϭ0.0085.Note that the oscillations are due to the advecting of in and out of the volume over which rms is calculated.
present throughout this calculation, its effects are hidden first by the pressure-driven mode and later by the development of a large poloidal flow.
The exponential growth that begins around t R ϭ0.0085 is the result of a large, Reynolds-stress-driven poloidal flow.Recall that the classical viscosity present in the vorticity equation provides little dissipation to damp this flow.In addition, because we have ignored the poloidal flow damping that enters in the momentum equation via the viscous stress tensor, the flow continues to build and eventually terminates the calculation.The coherent oscillations seen during the growth of the poloidal flow are due to the region of large moving in and out of the volume over which the average is taken.
Figures 15͑a͒-15͑c͒ show contour plots of during the linear, saturated, and late growth stages.Note the fine scale structure that develops between t R ϭ0.005 and 0.0075.As nonlinear coupling continues in the saturated phase, small eddies generate a large mϭ1, nϭ0 poloidal flow which continues to grow past t R ϭ0.01 where the large mϭ1 component is evident ͓see Fig. 15͑c͔͒.It is important to note that the generation of mϭ1 poloidal flow is unique to the Reynolds stress in toroidal geometry; cylindrical Reynolds stress can only generate an mϭ0 component.
That even purely growing modes can lead to a poloidally asymmetric Reynolds stress can be shown as follows.With ͑r,,͒ϭ ͚ m,n m,n ͑ r ͒sin͑ mϩn ͒, ͑10͒ the Reynolds stress can be written as where ͗͘ denotes an average over toroidal angle.
The poloidal flow and flow shear are shown in Figs.16͑a͒ and 16͑b͒ at time t R ϭ0.01.Note that the dominant region of flow shear coincides with the inner region of weak magnetic shear ͑see Fig. 1͒.Because poloidal flow shear has been used to explain reduced core transport in enhanced confinement regimes, 14 the fact that small-scale MHD activity can produce large-scale, sheared poloidal flows is an important one.Recall, however, that by ignoring the viscous stress tensor and including only small levels of viscous and diffusive dissipation, we have tacitly assumed a vanishingly small damping coefficient.Realistic damping may have reduced or prevented the development of this large poloidal flow.In addition, since the parallel momentum equation is not evolved, this calculation ignores the coupling between the poloidal and toroidal flows.In order to show that this poloidal flow generation was not driven by an external source, this calculation was repeated without the equilibrium toroidal flow yielding essentially the same result.This suggests that sub-Alfve ´nic toroidal flows do not affect the development of large, poloidal flows.
Because the effects of the tearing mode were hidden by the pressure-driven instability and the poloidal flow, ␤ was reduced and the equilibrium (nϭ0) components of U and were not evolved in further calculations.In contrast to the previous case that ignored poloidal flow damping, here it is assumed that the damping coefficient is infinite.The relevant parameters are Sϭ10 5 , ␤ϭ1.5%, and v /v Hp ϭ0.005, and 0.01.The evolution of rms is shown in Fig. 17.At lower ␤ and v /v Hp ϭ0.005, the weaker, pressure-driven instability FIG. 15.Contours of at t R ϭ0.005, 0.0075, and 0.01.In the linear phase (t R ϭ0.005), appears well localized.Nonlinear coupling in the saturation phase (t R ϭ0.0075) produces a Reynolds-stress-driven poloidal flow dominated by the mϭ1, nϭ0 component (t R ϭ0.01).
saturates nonlinearly around t R ϭ0.0075 compared to t R ϭ0.005 in the previous case ͑see Fig. 14͒.Here the tearing mode at rϭ0.75 grows linearly until t R ϭ0.007 before it is overcome by the activity at the inner region of positive shear, where again the sharp increase in rms is caused by the pressure-driven mode suddenly being included in the volume average.Upon saturation of the pressure-driven instability, the resistive tearing mode emerges again around t R ϭ0.01, grows at the same growth rate, and saturates around t R ϭ0.014.Note that the growth rate of the tearing mode in this nonlinear calculation is nearly the same as the growth rate for the dominant 7/3 harmonic seen in linear calculations for the same toroidal flow shear.This suggests that the nonlocal coupling between the pressure-driven and tearing instabilities is rather weak.
In order to determine the effect that sheared toroidal flow has nonlinearly, the above calculation was repeated with v /v Hp doubled to 0.01.In accordance with linear results, Fig. 17 shows that the tearing mode is more unstable at larger flow shear.In addition, the tearing mode appears to be even less affected by the presence of the pressure-driven instability than in the case with v /v Hp ϭ0.005.This suggests that the flow shear weakly decouples the two regions of plasma over which rms is calculated.The evolution of the pressure-driven instability (rϭ0.25) was similar in both cases.
Figures 18͑a͒-18͑c͒ show the evolution of the contours at times t R ϭ0.01, 0.0125, and 0.015 for v /v Hp ϭ0.005.Here the 7/3 helicity is dominant as evidenced by the number of vortices visible in Figs.18͑b͒ and 18͑c͒.As mentioned earlier, the tearing mode prefers the 7/3 surface because of its proximity to the region of zero shear.Figure 18͑b͒ shows evidence of the resistive/interchange ͑low S͒ or ''double tearing'' ͑high S͒ structure ͑seen in linear eigenmodes at ␤ϭ1.5%͒ as the middle 7/3 surface appears to be unstable as well.At t R ϭ0.015, the 5/2 helicity appears and the mode takes on an infernal character near the outer region of weak magnetic shear. 15,16he nonlinear results presented here show that smallscale MHD activity can produce large, sheared poloidal flows in the absence of poloidal flow damping.In addition, it is seen that sub-Alfve ´nic, sheared toroidal flows have a minimal effect on weakly coupled, localized instabilities possessing a common helicity.

V. CONCLUSIONS
In this work we have used a set of reduced MHD equations to study the linear and nonlinear stability of a nonmonotomic q profile in the presence of an equilibrium, sheared toroidal flow.Linear results regarding the effects of sheared toroidal flow and finite pressure on tearing mode growth rates agree qualitatively with previous models that treated the effects separately.More importantly, however, they point to the importance of considering both effects simultaneously.Namely, the stabilizing effect attributed to finite-␤ depends upon the level of toroidal flow shear even at sub-Alfve ´nic flow speeds.Modest flow shear spoils the decrease in growth rates normally seen with increasing ␤ and finite-␤ enhances the destabilizing effect of increasing flow shear.Thus the curvature term responsible for stabilizing tearing modes is also responsible for the destabilizing effect at small but finite shear via coupling to the advective term in the pressure equation.At larger toroidal flow shear, the effect of the advective term in the vorticity equation included in previous models is weakly stabilizing.
Nonlinear calculations show the generation of a large, Reynolds-stress-driven poloidal flow due to small-scale MHD activity in the absence of flow damping.This result suggests numerically that negative shear profiles are capable of producing the sheared poloidal flows used to explain the transport barrier in enhanced core confinement regimes.In the presence of infinite flow damping, nonlinear calculations show that modest levels of toroidal flow shear have a minimal effect on weakly coupled, localized instabilities.

FIG. 9 .
FIG.9.Plot of ␥ Hp versus log s at ␤ϭ0.0% and 1.5% for Sϭ10 5 .Increasing flow shear reduces the stabilizing effect of finite ␤ making growth rates comparable around sϭ15.At large s, growth rates at zero and finite pressure coincide.

FIG. 16 .
FIG.16.Contours of the poloidal flow ͑a͒ and flow shear ͑b͒ at t R ϭ0.01.Note the dominant region of poloidal flow shear resides in the inner region of weak magnetic shear.