Toroidal Geometry Effects in the Reversed-field Pinch MHD Dynamo

Carl R. Sovinec, Richard A. Nebel

Los Alamos National Laboratory

Dalton D. Schnack

Science Applications International Corporation

presented at the

40th Annual Meeting of the APS Division of Plasma Physics

New Orleans, Louisiana

November 16-20, 1998

Outline

I. Introduction

A. Overview

B. Background

C. Goals

II. The NIMROD Code

A. List of active participants

B. Physical model

C. Numerical algorithm

IV. RFP simulations

A. Physical conditions and numerical

parameters

B. Geometry comparison

C. 2-fluid effects

V. Preliminary conclusions and future work

Addendum: Spheromak simulation

Overview

Nonlinear, three-dimensional fluid computations are used to investigate magnetohydrodynamic (MHD) activity in reversed-field pinches (RFP). Results pertaining to toroidal geometry effects at small aspect ratio and two-fluid modifications to quasilinear interactions in a periodic cylinder are presented.

The NIMROD (Non-Ideal MHD with Rotation) code used for this study features a finite-element representation of the poloidal plane and a semi-implicit temporal advance. Simulations on massively parallel computers are made possible through two types of domain decomposition and message-passing communication.

We find that toroidal curvature plays only a moderate role in RFP MHD activity, even in simulations with R/a<2. Simulations in linear geometry show that Hall and electron inertia terms tend to enhance quasilinear coupling.

As an addendum to the main topic, initial work on studying nonlinear MHD activity in spheromaks is also presented.

Background

The MHD dynamo is an important feature of RFPs. Ohmic current drive with weak (relative to tokamaks) toroidal magnetic field leads to a peaked current profile that is unstable to current-gradient driven modes. The modes saturate due to nonlinear coupling to stable modes and due to quasilinear feedback on the mean profile. The quasilinear feedback involves a mean electric field, , that sustains parallel current where the magnetic field is purely poloidal. This allows reversed toroidal field to persist near the wall in the presence of finite resistivity.

While previous numerical studies have examined the MHD activity in detail, all (?) have been performed in periodic, linear geometry. This is usually a good approximation, since the safety factor is less than unity across the profile. However, linear coupling associated with toroidal curvature can affect the MHD activity at low aspect ratio.

An earlier numerical investigation of aspect ratio [Ho, Phys. Plasmas 2, 3407] in periodic, linear geometry shows that the fluctuation energy tends to concentrate into a single Fourier component at low aspect ratio. This has led to speculation that such a configuration may be less stochastic than conventional RFPs.

Goals

The purpose of this work is to determine what role toroidal curvature plays in the MHD activity in RFPs and the extent of its significance. Besides continuing the development of our understanding of RFPs, this will also provide bounds of validity for the linear geometry approximation.

Another goal of our study is to provide quantitative information regarding magnetic stochasticity and magnetic diffusivity at small aspect ratio in the full toroidal geometry. Numerical tools for computing the appropriate measures are being developed.

The NIMROD Code

LIST OF PRESENTLY ACTIVE TEAM MEMBERS

ADVISORS

PHYSICAL MODEL

The physical foundation of NIMROD lies in the two-fluid/ Maxwell system of equations.

Quasineutrality and fluid moments:

Maxwell's (no displacement current):

where a=i,e,

and primes indicate the center of mass reference frame,

 The separate electron and ion momentum equations are combined into a single-fluid form for numerical convenience. Only the net charge forces are approximated away, assuming low-frequency behavior.

Additional notes regarding the equations advanced in NIMROD:

Closure terms may be modified to suit a variety of physical situations, e.g., a linearized Hirshman-Sigmar closure has also been added to Ohm's law as an initial step towards neoclassical physics in tokamaks.

In its present form, NIMROD solves all of the low-order fluid moments, except continuity; however, the current advection terms in Ohm's law and all nonadiabatic terms in the pressure equations have not been implemented, and the stress tensor term in the momentum equation is presently an isotropic kinematic viscosity. The Hall, electron inertia, and stress tensor terms may be omitted in a simulation to solve the resistive MHD equations.

NUMERICAL ALGORITHM

NIMROD represents functions of space with 2-D finite elements for the poloidal plane and Fourier series for the perpendicular direction (toroidal or periodic linear).

The poloidal plane is decomposed into an unstructured collection of blocks. This permits geometric flexibility and creates a natural decomposition for running on parallel machines, where MPI is used for communication between processors.

The blocks containing triangles are themselves unstructured and may be used for complex boundary shapes or for singular points in the grid.

Triangulation of the 'tblocks' is done with the TRIANGLE code written by Jonathan R. Shewchuck at Carnege Mellon University. TRIANGLE may be downloaded from http://www.cs.cmu.edu/~quake/triangle.html .

Sample Grid for D3D Simulations

More on the NIMROD spatial discretization:

• Blocks of structured grid (quadrilaterals) are aligned with flux surfaces in the core region to resolve tearing layers, while blocks of unstructured grid permit decomposition near complex wall shapes.
• • the grid and equilibrium fields are represented by cubic splines in structured blocks.
  • dependent variables are treated as bilinear within quadrilateral cells.
  • all quantities are linear in triangular cells.

• Finite elements allow a straightforward procedure for decomposing complex domains.
• • discretization is based on functional form instead of difference operators.
  • the order of spatial convergence depends solely on the chosen functional form.
  • symmetry of differential operators is preserved.

• The disadvantages of finite elements include
• • consistency of differential operators is not guaranteed.
  • added coding complexity

• Nonlinear terms are treated in a pseudospectral fashion
• • FFTs are used to transform to configuration space where nonlinear products are computed

A semi-implicit algorithm [Schnack, J. Comput. Phys. 70, 330, for example] is used to advance the solution from initial conditions.

Simple example of a 1D linear sound wave:

A semi-implicit scheme replaces mass density in the momentum equation with

After normalizing pressure by P₀ , velocity by the sound speed, and assuming an e^ikx dependence, the modified leap-frog advance has the form,

where

The complex eigenvalue representing the linear operations of the time step,

satisfies

where

The magnitude of l is unity (advance is numerically stable without dissipation) if

Setting the coefficient

ensures stability for all Dt.

The choice of semi-implicit operator has a strong influence on accuracy at large Dt (large compared with propagation times of normal modes).

In NIMROD, our semi-implicit scheme for the MHD waves leads to the following equation for advancing velocity:

After integration by parts in the finite element procedure, the operator has the form of the linear potential energy of the ideal MHD system.

This allows time steps determined by the phenomena of interest, even if they are large with respect to the Alfven time (important for large S conditions).

Caveat: all advection is treated in a predictor/corrector fashion, and large flows impose a CFL restriction on Dt.

Hall terms, when used, are time-split from the MHD advance of B, and require a separate semi-implicit operator. [Harned, J. Comput. Phys. 83, 1]

NIMROD uses two types of problem decomposition and message-passing communication to allow efficient computation on hardware ranging from massively parallel supercomputers to multi-processor workstations.

This scaling has been performed on the NERSC T3E for a few time steps of an RFP simulation similar to those reported below, albeit with twice the number of Fourier components (86 vs. 43 below). Simulations in the scaling run on 86 processors or less use a single grid block. Simulations with more than 86 processors use 86 processor layers and multiple grid blocks.

The RFP simulations reported below have been run with either 43 or 86 processors.

Our equations are written in cylindrical components (about the geometric axis), so simulations of compact toroids (not RFPs) require regularity conditions.

From a 2D Taylor series in and of an arbitrary scalar or vector function near the geometric axis, the following must be preserved.

Scalars: as for n>0. (s)

Vectors: as for n>0. (z)

• The homogeneous Dirichlet conditions (s, z, and rp1) are enforced as essential boundary conditions.

• The n=1 rotation condition (rp3) is enforced in each matrix equation by summing the equations for and and for and , setting the difference of each pair to zero, and solving for the average.

• The Neumann conditions (rp2) are satisfied through a volume constraint integral, added to each time-advance equation for a radial and azimuthal n=1 vector component.

The eigenfunction of a nonresonant, ideally unstable kink mode in a periodic cylinder, where the r-z plane is gridded, demonstrates the effectiveness of the penalty function.

without penalty<----------->with

RFP Simulations

PHYSICAL CONDITIONS AND NUMERICAL PARAMETERS

Toroidal geometry study:

• resistive MHD terms in Ohm's law (no Hall or electron inertia terms)
• 0-b conditions
• S=2500 (t_{r set to 1s)}
• P_m=m₀n/h=1
• circular cross section
• initial conditions for the RFP simulations are provided by separate n=0 simulations of a paramagnetic pinch.
• 48 X 48 poloidal grid resolution (based on a convergence study of the linear modes), and n<43. For toroidal cases, the center of the grid is adjusted to match the Shafranov shift of the paramagnetic pinch.
• conducting wall boundary conditions, except for an effective gap for applying a loop voltage.
• the loop voltage is adjusted to approximately maintain a constant current:

Unstable modes grow until nonlinear and quasilinear interactions cause saturation. Reversal of the toroidal magnetic field results from quasilinear feedback.

Temporal evolution of the reversal parameter and the magnetic energies in the different toroidal Fourier components from the toroidal simulation with R/a=1.5 and pinch parameter of 1.6.

The MHD activity flattens the parallel current profile and drives q(a) below 0. The green traces show the initial paramagnetic state; other colors show subsequent times. [These plots show flux-surface averaged quantities that are based on the n=0 Fourier component, only.]

Comparison of Global Results

• Results are averaged over 0.12-0.2 t_{r after saturation, except for simulation L3, which is being continued.}
• 'T' and 'L' refer to toroidal and periodic linear geometry, respectively.
• The pinch parameter is a measure of the toroidal current

• The reversal parameter is

For the R/a=1.5 cases at low Q (T1 and L1), the linear simulation displays more MHD activity and more reversal than the toroidal simulation. At high Q, the toroidal cases generate more reversal. However, these toroidal cases (T2 and T3) tend to a near-steady, single-helicity state.

Temporal averages and ± one standard deviation are plotted for the R/a=1.5 cases.

Case T2 slowly emerges from its single-helicity state, and by the end of the simulation, the field is stochastic--possibly more so than case L2.

We plan to compute Liapunov exponents and magnetic diffusivities to compare the different geometries and aspect ratios quantitatively.

TWO-FLUID EFFECTS

Hall terms have been incorporated in the generalized Ohm's law, though gyro-viscous terms have not. We are presently testing the incorporated two-fluid effects in linear and quasilinear simulations.

•The Hall parameter is large:

•Rotation of the mode near the diamagnetic frequency is observed.

•Rotation decreases substantially in the 0-b limit.

•From a reduced MHD framework, most of the rotation likely results from the magnetic curvature correction, , which is large in this reversed-field equilibrium.

Resistive MHD shows only island formation, no rotation.

With Hall and electron inertia, the island rotates as it develops.

• Low R/a simulations of RFPs in toroidal geometry appear to have a greater tendency to develop single-helicity states than simulations in periodic linear geometry.

• Where single helicity states are not formed, toroidal curvature appears to diminish MHD activity, including field reversal.

• Where single-helicity states are formed in toroidal geometry, there is more reversal than in the comparable (multi-helicity) linear geometry simulation.

• More simulations are needed to verify these findings. Simulations at larger S values are also necessary.

• Numerical means of computing quantitative measures of field stochasticity are being developed.

• The Hall and electron inertia terms in Ohm's law will be applied in full RFP simulations.

ADDENDUM: Initial efforts to simulate spheromak formation and sustainment

Simulation parameters:

• resistive MHD
• 0-b
• S=1000 (t_{r set to 1s)}
• P_m=m₀n/h=1
• the geometry is a cylinder 1m in height and radius.
• the ends of the cylinder form electrodes through which magnetic field and current pass. The field is line-tied to the electrodes.
• a stabilized pinch with large current, l(0)=10m^-1, is used as an equilibrium.
• the equilibrium is unstable to an ideal kink-type mode (gt_A~20%).
• the poloidal grid is 40x40, and Fourier components with n<6 are included in the nonlinear simulation.

The n=1 eigenfunction shows the traits of a kink mode.

The unstable mode saturates in a steady-state (at least at this low value of S), coupling with the larger n components and the mean field.

This results in the formation of closed contours of poloidal flux (based on the toroidal average field).

The poloidal flux amplification resulting from the MHD activity is ~100% of the initial open poloidal flux in this case.