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  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.
In two dimensions, with incompressible flow, the MHD equations can be written
where the two dimensional Laplacian is
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.
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 . 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  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.
Figure 3: Part of final mesh
The peak value of the current density grows exponentially in time ,, 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: firstname.lastname@example.org Acknowledgment
This work was supported by AFOSR Grant No. 91-0044 and U. S. Department of Energy Grant DE-FG02-86ER53223.