An Adaptive Finite Element Method for Magnetohydrodynamics
H. R. Strauss, NYU
D. W. Longcope, U. of Montana
Introduction
A finite element discretization for two
dimensional MHD [1] is described.
The elements are triangles with
piecewise linear basis functions.
The main computational difficulty is the accurate calculation
of the current. The most effective solution is to
employ
a current - vorticity advection formulation
of the equations. Acceptable results can also be obtained with
a two - step calculation of the current from the vector potential.
Mesh operations are
described to reconnect and refine the mesh adaptively in the
vicinity of nearly singular currents, to improve magnetic flux
conservation.
Example computations of the coalescence instability,
tilt mode and divertor tokamak equilibrium,
validating and illustrating the method,
are presented. The simulations show the formation of current
sheets, with the current density increasing exponentially in time.
During this increase, the grid of initially 
points
adapts to provide resolution comparable to a uniform grid of
up to 
grid points.
MHD
In two dimensions, with incompressible flow, the MHD equations can be written
where the two dimensional Laplacian is
and
The left hand side of
(1), along with (3) is the familiar
vorticity - stream function formulation
of two dimensional incompressible hydrodynamics. The right hand
side of (1) comes from the Lorentz force
with current density C,, and
the viscosity with coefficient 
.
Eq. (2) is from Ohm's Law and Faraday's Law and
represents the conservation of magnetic flux 
The MHD equations conserve magnetic flux, and energy (if the
viscosity 
) Since the magnetic
flux function is advected with the flow, any function of 
is a constant of the motion.
We assume either Dirichlet boundary conditions, with 
constant on the boundary, Neumann conditions with the normal
derivatives of equal to zero, or periodic boundary
conditions.
Symmetrization of Equations
The evolution of the magnetic and velocity fields are treated
in a non symmetric way in the standard formulation above.
The velocity is advanced through the vorticity, while the
magnetic field is advanced via the magnetic potential.
This can cause numerical problems when the equations are
solved on an irregular mesh, such as the unstructured, adaptive
grids we will use later.
It is desirable to formulate
the equations in a more symmetrical manner, in which
the current and vorticity are time advanced. This eliminates
numerical
problems arise from taking the Laplacian of the magnetic potential and
inserting it into the vorticity equation. Solving the equations
as they are leads effectively to third derivatives of
on the right side of (1), which is troublesome using
piecewise linear basis functions.
The best way to avoid this is by using the current advection form, in which no second derivatives appear on the right hand sides of the equations, except perhaps in dissipative terms.
Instead of solving eqs. (2) and (4), we take the Laplacian of (2) and use (3), (5), and (6). This yields an equation for the current, analogous to (1) for the vorticity,
Eqs. (7), (8)
are solved along with eqs. (1) and (3).
The equations are now symmetrical, in the sense that the source
functions 
and C are time advanced, and the
potentials 
and are obtained at each time step
by solving Poisson equations (3) and (8).
Dirichlet, Neumann, or periodic boundary equations are
applied to the potentials.
Finite Element method
We use a finite element discretization,
introducing
basis functions 
, which in the
present work are piecewise
linear over each triangle, and satisfy
The variables in the MHD equations are represented
as a sum over basis functions.
We use a zero residual Galerkin approach
in which
the equations are multiplied by a
basis function 
, and integrated over the domain.
This gives
a set of sparse matrix equations.
It is convenient to consider a diagonalized form of the
mass matrix, called
the lumped mass matrix, formed
by subtracting all the off diagonal matrix elements in each row
and adding them to the diagonal.
Using the lumped mass, the finite element discretization is equivalent to a finite volume discretization, where the control volumes are constructed by joining the barycenters of the triangles (average of the vertex positions) to the midpoints of the triangle edges.
Adaptive Refinement
We have an adaptive algorithm
to refine the mesh as the computation proceeds.
Because we want to resolve current sheets,
we monitor
the current.
If the product of current density times triangle area,
exceeds a threshold, we split the triangle in
two. The local current density in a current sheet typically
rises exponentially in time, so the refinement process has this
same behavior. We have to stop the refinement at some preset
triangle size and number of mesh points. Similarly, if the
current is too low, the triangles at that
vertex are
unrefined.
The justification for using the product of current density
and triangle area as a refinement criterion, is that
is comparable to the discretization error
in the magnetic flux
It is very important to conserve magnetic flux in the
neighborhood of a current sheet, in order to prevent numerical
reconnection.
Computational Results: Tilt Mode
We consider the two dimensional tilt instability [2]. In this calculation we have used the lumped mass matrix and the current - vorticity formulation of MHD.
The initial equilibrium state is a bipolar vortex,
When perturbed, an instability occurs, growing exponentially
as 
We perform a simulation with an initial
mesh, with 40 triangles on a side.
Starting with the
equilibrium of eq.(9),
a perturbation about 
smaller is inserted.
Figure 1: (a) initial magnetic flux (b) initial current
In the simulation, we take 
and the
simulation box has sides of length 4. Conducting boundary
conditions are applied on the walls, at which 
and 
The previous simulations
[2] were
compressible, and growth rates were reported in the
range 
depending on the pressure.
None of these cases are exactly equivalent to our
strictly incompressible model.
We obtain the linear growth rate 
Adaptive simulations were done with the current advection scheme. In the simulation, the motion is highly nonlinear by time t = 7. At this stage, the vortex has tipped over. The separatrix wraps around the two flux vortices. Current sheets are formed at the leading edges of the central vortices.
Figure 2: (a) final magnetic flux (b) final current
The mesh supporting the contours is shown in a blowup view The mesh resolution has adaptively followed the formation of the moving, curved current sheet.
The peak value of the current density grows exponentially in time [3],[4], with a large growth rate about 3 times the linear mode growth rate. This can be seen in the next figure, which shows the logarithm of the peak current density as a function of time. The logarithm of the peak current density grows approximately linearly. The two curves are from runs with different refinement parameters. As the current density increases, so does the number N of mesh points. Again, the upper curve had the lower refinement threshold, which caused more refinement both in the initial and latter stages of the simulation. The runs ended with 73,700 and 23,670 mesh points, respectively.
Figure 4: (a) 
as a function of time(b) Log of number of mesh points
as a function of time.
email: strauss@cims.nyu.edu Acknowledgment
This work was supported by AFOSR Grant No. 91-0044 and U. S. Department of Energy Grant DE-FG02-86ER53223.
