next up previous
Next: References

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 tex2html_wrap_inline334 points adapts to provide resolution comparable to a uniform grid of up to tex2html_wrap_inline336 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 tex2html_wrap_inline344 . Eq. (2) is from Ohm's Law and Faraday's Law and represents the conservation of magnetic flux tex2html_wrap_inline346

The MHD equations conserve magnetic flux, and energy (if the viscosity tex2html_wrap_inline348 ) Since the magnetic flux function is advected with the flow, any function of tex2html_wrap_inline350 is a constant of the motion. We assume either Dirichlet boundary conditions, with tex2html_wrap_inline352 constant on the boundary, Neumann conditions with the normal derivatives of tex2html_wrap_inline352 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 tex2html_wrap_inline350 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 tex2html_wrap_inline358 and C are time advanced, and the potentials tex2html_wrap_inline362 and tex2html_wrap_inline350 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 tex2html_wrap_inline366 , which in the present work are piecewise linear over each triangle, and satisfy tex2html_wrap_inline368

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 tex2html_wrap_inline370 , 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, tex2html_wrap_inline372 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 tex2html_wrap_inline374 is comparable to the discretization error in the magnetic flux tex2html_wrap_inline346 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 tex2html_wrap_inline380

We perform a simulation with an initial mesh, with 40 triangles on a side. Starting with the equilibrium of eq.(9), a perturbation about tex2html_wrap_inline384 smaller is inserted.

Figure 1: (a) initial magnetic flux (b) initial current

In the simulation, we take tex2html_wrap_inline386 and the simulation box has sides of length 4. Conducting boundary conditions are applied on the walls, at which tex2html_wrap_inline390 and tex2html_wrap_inline392

The previous simulations [2] were compressible, and growth rates were reported in the range tex2html_wrap_inline394 depending on the pressure. None of these cases are exactly equivalent to our strictly incompressible model. We obtain the linear growth rate tex2html_wrap_inline396

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 [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) tex2html_wrap_inline408 as a function of time(b) Log of number of mesh points as a function of time.

email: Acknowledgment

This work was supported by AFOSR Grant No. 91-0044 and U. S. Department of Energy Grant DE-FG02-86ER53223.

next up previous
Next: References

Hank Strauss
Mon Nov 22 17:03:50 EST 1999