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






I. Introduction


A. Overview

B. Background

C. Goals




A. List of active participants

B. Physical model

C. Numerical algorithm


IV. RFP simulations


A. Physical conditions and numerical


B. Geometry comparison

C. 2-fluid effects


V. Preliminary conclusions and future work



Addendum: Spheromak simulation





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.





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.





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.








Ahmet Aydemir IFS

Curtis Bolton OFES

Ming Chu GA

Thomas Gianakon U-WI

Alan Glasser LANL

Scott Kruger U-WI

Jean-Noel Leboeuf UCLA

Richard Nebel LANL

Steve Plimpton SNL

Nina Popova MSU

Dalton Schnack SAIC

Carl Sovinec LANL

Alfonso Tarditi SAIC





James Callen U-WI

John Finn LANL






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.





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 .




Sample Grid for D3D Simulations




More on the NIMROD spatial discretization:






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 P0 , velocity by the sound speed, and assuming an eikx dependence, the modified leap-frog advance has the form,





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








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.


• Decomposing the computational grid into blocks permits the first type of problem decomposition. Communication is performed among processors owning neighboring blocks during finite element computations and during matrix-vector multiplication (in conjugate gradient iterations for solving the matrix equations).


• The semi-implicit advance leads to separate matrix equations for each Fourier component. Thus, decomposition of the components among 'layers' of processors makes for efficient parallelization. FFTs and pseudospectral operations are handled in a scalable manner.


• Typical preconditioning for the conjugate gradient technique tends to lose effectiveness as the number of grid blocks is increased. We are presently working on a multi-level scheme to alleviate the difficulty (see Glasser, poster G4Q.11).



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)

as for n=0 and n>1. (rp1)

as , (rp2)

and as for n=1. (rp3)






where a is the linear or bilinear test function (triangular or quadrilateral elements, respectively) multiplied by the radial or azimuthal unit vector, and c is a linear or bilinear scalar function that has a uniform value at all vertices along the r=0 border of the domain and 0 at all other vertices. Thus, the modification serves to penalize any change in the radial derivative of these components. The influence of this additional integral vanishes as radial resolution is increased.


Note that a standard finite-element treatment of the (rp3) condition using natural boundary conditions is ineffective, since the 'surface integral' is identically 0, regardless of the radial derivatives.



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




Toroidal geometry study:

where i is the toroidal current, i0 is the desired current, and c0 and c1 are constants.


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











± 0.002

34.0± 1.0


± 0.006




± 0.002

32.8± 1.3


± 0.015




± 0.004

46.3± 2.8


± 0.007




± 0.004

45.2± 2.4


± 0.018




± 0.008

35.1± 5.6


± 0.022




± 0.001

33.2± 1.0


± 0.021



where A is the poloidal cross-section area, F is the toroidal flux, and dl goes around the cross section.

where ishell is the poloidal current in the shell.



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.





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.


Consider an m=1, kz=2 ideally unstable mode in periodic linear geometry with an equilibrium defined by pitch



and pressure





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.




Preliminary Conclusions and Future Work











ADDENDUM: Initial efforts to simulate spheromak formation and sustainment


Simulation parameters:



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.