23P128-B
Concept for Numerical Calculation of
3D MHD Equilibria with Flow and FLR Effects
流れと有効ラーモア半径効果を取り入れた
3次元MHD平衡の数値解析の進展
Daniel Raburn, Atsushi Fukuyama
レイバン ダニエル, 福山 淳
Kyoto University, Department of Nuclear Engineering
Sakyo-ku, Yoshida-Honmachi, Kyoto University, Kyoto 606-8501, Japan
京都大学 606-8501 京都市左京区吉田本町
One of the simplest models for describing plasma equilibrium is the single-fluid magnetohydrydynanimc(MHD) model under axisymmetry with no fluid flow. However, both non-axisymmetric effects and fluid
flow can significantly alter plasma equilibrium. A code which can rapidly calculate 3D MHD equilibriumwith fluid flow is important for future modeling of fusion plasmas. Presently, there are many codes which
can look at a subset of the physics of interest. Based on those codes, a new code for calculating 3D MHD
equilibrium with fluid flow is being developed.

1. Overview of MHD Equilibrium
Single-fluid MHD equilibrium is governed byconservation of mass, Ohm's law, force-balance,and Maxwell's equations with quasi-neutrality.
Closure requires an equation of state, often taken tobe an adiabatic law. Consider a cylindricalcoordinate system (R, , Z), where R is the major
radius,  is the toroidal angle, and Z is the distanceabove the mid-plane. For the case of anaxisymmetric plasma (/ = 0) with no bulk fluid
flow, equilibrium is governed by the well-known
Grad-Shafranov equation [1]:∇⋅1
1 dI 2
dp
∇

=02
2
R
2R d d(1)where  is the poloidal magnetic flux, I is thepoloidal current, and p is the pressure. Equilibrium
depends on the boundary conditions for , as wellas the two free functions p() and I().
With the inclusion of fluid flow, equilibrium isgoverned by a generalized Grad-Shafranov
equation and a Bernoulli equation [2]. Equilibriumin this case depends on five free functions of thepoloidal flux, which can be taken to be the fluxsurface averages of the following quantities: (1) the
poloidal current I; (2) the pressure p; (3) the densityn; (4) the electrostatic scalar potential ; and, (5)the stream function , which measures the poloidal
flow [3].
However, the inclusion of fluid flow can produce
a sharp profile pedestal [2-6]. In the single-fluid
MHD model, the pedestal occurs where the poloidal
flow transitions from sub- to super-poloidal sonic.
However, because of the small scale-length of thepedestal, the inclusion of finite Larmor radius(FLR) effects are generally necessary to accurately
model the plasma. Ito and Nakajima havedeveloped a formulation for MHD equilibrium withflow and FLR effects using an expansion in the
inverse aspect-ratio (), which includes FLR effects
by the ordering /a ~ , where  is the ion Larmor
radius and a is the plasma minor radius [6]. In thisformulation, equilibrium depends on five freeprofiles of the lowest-order poloidal flux, which arerelated to the lowest non-constant order of thefollowing quantities: (1) the poloidal current I; (2)the electron pressure pe; (3) the ion pressure pi; (4)the density n; and, (5) the electrostatic scalarpotential  [7]. However, because the effects offlow are included as a linear perturbation, thepedestal becomes a singularity, and the formulation
cannot be directly applied to trans-poloidal-sonicflows.
In 3D cases, even with no flow, equilibrium iscomplicated by the absence of a guarantee of nested
flux-surfaces: islands and stochastic regions arepossible. Even in nominally axisymmetric devices,magnetic islands may form through tearing modes.
3D MHD equilibrium without flow depends on two
free functions, which must be constant on eachmagnetic field line: (1) the field-line average of theparallel electrical current density divided by themagnetic field strength <J||/B>B; and, (2) thepressure p [8].
2. Overview of Existing Codes
There presently exist many codes for calculating
MHD equilibrium in a variety of cases. Here, webriefly review three codes: (1) the FLOW code; (2) acode for calculations under the formulation by Itoand Nakajima; and, (3) the PIES code.
The FLOW code is capable of calculating
axisymmetric single-fluid MHD equilibrium withflow [2]. FLOW handles the Bernoulli equation
explicitly, identifying the location of the transpoloidal-sonic surface, and separates thecomputational domain in to two regions, with theouter region assumed to have super-poloidal-sonicflow.
We have written a solver for reduced two-fluid
MHD equilibrium in the large aspect-ratio
axisymmetric limit, as per the formulation by Itoand Nakajima [7]. Because of the limitations of theformulation, the code cannot handle trans-poloidalsonic flows.
The PIES code is capable of calculating 3D
single-fluid MHD equilibrium without flow [8].
The algorithm is based on that originally suggested
by Spitzer [9] and by Grad and Rubin [10], andinvolves iterated direct integration of theequilibrium equations. However, there aresubstantial complexities involved with determining
the topology of the magnetic field lines, and PIES
uses a system of magnetic coordinates to effectively
treat magnetic islands and stochastic field lines. As
a consequence of the complications, PIES is
prohibitively slow for some applications.
Recently, an external wrapper for the PIES code
has been developed to improve the speed ofconvergence [11]. The external wrapper makes useof a custom Jacobian-free Newton-Krylov (JFNK)
numerical solver with adaptive preconditioning anda novel subspace-restricted Levenberg-Marquardt
backtracking algorithm. The details of the solver are
beyond the scope of this document, but for
information on JFNK, see Ref. [12], and, for
information on adaptive preconditioning, see Ref.
[13].
3. Concept for a New CodeWe have begun work on a code for calculating
3D MHD equilibrium with fluid flow and FLReffects. We refer to the code as KITES – Kyoto
ITerative Equilibrium Solver. Based on an analysisof the governing equations, we believe thatequilibrium is governed by the magnetic field-line
averages of six quantities, which can taken to be thefollowing: (1) the parallel current density J||; (2) theelectron temperature Te; (3) the ion temperature Ti;(4) the density n; (5) the electrostatic scalarpotential ; and, (6) the parallel fluid velocity v||.
For a numerical solver, we plan to use the custom
JFNK solver which was used with PIES. Unlike

PIES, we hope to work entirely in physicalcoordinates without the use of magnetic
coordinates: this is expected to reduce the speed of
the solver but increase the robustness andmaintainability.
There appear to be three major hurdles to thedevelopment of the code. The first is to have a way
to accurately and rapidly determine the topology ofthe magnetic field lines. We have already developed
an algorithm for field line tracing in physical space,
and the results appear to be competitive with PIES.
The second major hurdle is accurately solving thedifferential equation along each field line; these
determine the variation of several quantities alongeach field line. We hope to be able to do this bybrute-force by following the field lines for a
sufficient length. The third hurdle is handling thetrans-poloidal-sonic surface. The presence of FLR
terms may eliminate the pedestal, in which case thishurdle would be moot. Otherwise, we would
attempt a multi-region approach similar to that usedin the FLOW code.
Acknowledgments
This work was supported by Grant-in-Aid for
Scientific Research (S) (20226017) from JSPS,Japan.
References[1] J. Wesson: "Tokamaks," (Oxford University Press,New York, NY, 2004) p. 110, 3rd ed.
[2] L. Guazzotto, R. Betti, J. Manickam, and S. Kaye:
Phys. Plas. 11 (2004) 604.
[3] D. Raburn and A. Fukuyama: Plasma Fusion Res. 6(2011) 2403044.
[4] T. Taylor, H. S. John, et al: Plasma Phys. Control.
Fusion, 36 (1994) B229.
[5] R. Betti and J. P. Freidberg: Phys. Plas., 7, (2000)
2439.
[6] A. Ito, J.J. Ramos, and N. Nakajima: Plasma Fusion
Res. 3 (2008) 034.
[7] D. Raburn and A. Fukuyama: Phys. Plas. 17(2)(2010) 122504.
[8] H. S. Greenside, A. H. Reiman, and A. Salas: J.
Comput. Phys. 81(1) (1989) 102-136.
[9] L. Spitzer: Phys. Fluids, 1 (1958) 253.
[10] H. Grad and H. Rubin, in Proceedings, Second
United Nations International Conf. on the PeacefulUses of Atomic Energy (United Nations, Geneva)
Vol. 31, p. 190.
[11] D. Raburn: PhD Thesis (2011).
[12] D. Knoll and D. Keyes: J. Comput. Phys. 193(2004) 357.
[13] D. Loghin, D. Ruiz, and A. Touhami: J. Comput. &Applied Math, 189 (2006) 362-374.