**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.
**

**
**

Mon Nov 22 17:03:50 EST 1999