Divergence-free finite-difference method for 2D ideal MHD

The divergence-free finite-difference scheme for 2d ideal MHD using triangular unstructured staggered grids is described. In this approach the density, pressure and velocity fields are attributed to grid cells and magnetic field is defined by its normal components in cells edges. The HLLD approximate Riemann solver is used to compute the numerical flux of gas-dynamics variables at edges. It also provides the electrical field values used to compute the magnetic field. The magnetic field is interpolated into the cells using Raviart — Thomas basis functions. The algorythm is implemented in parallel numerical code. The method is tested on several well-known two-dimensional MHD problems.


Introduction
Models based on ideal magnetohydrodynamics equations are high-demanded in geophysical and astrophysical investigations [1]. They are widely applied in modelling of macroscale plasma behavior [2,3]. The magnetic field divergence-free condition is the crucial for such models [4,5]. Artificial numerical non-zero divergence of magnetic field leads to non-physical oscillations of solution and the magnetic charge generation. This problem solution demands special numerical techniques to be implemented during computational research [6,7,8].
We consider new finite-difference divergence-free numerical scheme for ideal two-dimentional magnetohydrodynamics (MHD). We use the unstructured staggered triangular grids. The grid cells are used as the control volumes for gas variables such as density, velocity and pressure. The plane magnetic field components are associated with cells edges and are defined by normal to edge vector component (like [8]). Under proposed division the gas variables are computed in every time layer using Godunov-type scheme with reduced HLLD numerical flux [9]. Then the magnetic field normal components are computed directly from the Faraday equation using the Stokes theorem by calculating the circulation of electrical field [10] obtained as the components of HLLD flux.
The Stokes theorem exploitation leads to multiple results. On the one hand it guarantees the divergence-free of numerical scheme without any additional steps [6]. On the other hand the direct approximation of the Faraday law using the desribed technique is the general approach applicable for any grid type. The triangular grid is used for example.

Governing equations
Ideal MHD equations system includes the conservation laws for gas mass, pulse, energy and Faraday law for magnetic field evolution: Here ρ is gas density, v = [u, v, w] T is velocity, p is pressure,Î is unity tensor, B = [B x , B y , B z ] T magnetic induction vector, e is the total gas energy excluding the magnetic field energy and including internal and kinetic energy of unity gas volume: γ is adiabatic index. Let all variables depend only on time t and plane coordinates x and y, all derivatives with respect to z are equal to zero. The system (1) -(4) should be closed by setting the initial and boundary conditions for density, velocity, energy and magnetic field. The conservation of magnetic field is prescribed by divergence-free condition Applying the divergence operator to (4) one can obtain It means that the ∇ · B will not increase in time and magnetis field comply with condition (5) if initial distribution is divergence-free.

Numerical scheme
To construct the numerical scheme we use the explicit first order time integration and the physical processes splitting technique. Under this paradigm on every time step to obtain new variables values one should compute • the gas fluxes (mass, pulse and energy densities fluxes) on every grid edge using HLLD approximate Riemann solver; • gas variables inside cells implementing gas fluxes from-edges-to-cell balance condition prescribed in equations (1) -(3) like it usually done in Godunov-type schemes; • normal components of magnetic field on grid edges using the Stokes theorem; • magnetic field inside the grid cells by divergence-free interpolation procedure.
Obtained in such approach scheme is monotonic, conservative and divergence-free. x' y' Figure 1. Base grid cell. (x, y) is basic system of axes, (x , y ) is local rotated system of axes.

Numerical fluxes
Consider the triangular cell and its edge with outer normal depicted in Fig. 1. Let x axis be normal to edge and y axis be tangential. In standart approach [11] to obtain numerical flux we can locally for edge central point reduce model to 1D approximation and rescribe equations (1)(2)(3)(4) in new system of axes (x , y ) as: Here U = [ρ, ρu , ρv , ρw, e , B x , B y , B z ] T is conservative variables vector including magnetic T is magnetic field vector in rotated system of axis (x , y ). The dependance of total pressure p T and total energy e on gas pressure p is given by formulas This hyperbolic system has seven eigenvalues that correspond to entropy, two Alfvenic and four magnetosonic (two fast and two slow) waves [1]: Here All eigenvalues are real and ordered: Equations system (6) is usually [11] used to construct approximate solver for Riemann problem with two states prescribed by base and neighbouring cells. We use HLLD numerical flux which is described in detail in [9]. It considers six states rising from discontinuity point in cell's edge and divided by five waves depicted in Fig. 2. States U * L , U * * L , U * * R and U * R are diveded by Alfvenic waves travelling with velocities S * L and S * R to the left and to the right direction respectively and by contact wave travelling with velocity S M . Index L designates ,,left" side or the base cell and index R designates ,,right" side or neighbouring cell for considered edge.  Figure 2. Plasma states after breakdown of a discontinuity for HLLD solver. Dark blue is for fast magnetosonic waves, red is for Alfvenic waves and green is for enthropy wave.
Let normal velocity component and total pressurt are unchanged by Alfvenic and contact waves: Here After S M and p * T calculation and using Rankine -Hugoniot condition for discontinuity [1], one can compute states U * α with α = L or R in following way: Alfvenic waves velocities are computed using formulas Here In formula for energies e * α and e * * α plus sign is suitable for α = L and minus one is for α = R. Finally the HLLD numerical flux F at the center of edge is computed as Values S L and S R one can compute as the totally minimal and maximal eigenvalues (7) for MHD equations system

Faraday equation approximation
During the numerical simulations based on system (6) the divergence-free condition (5) may fail. In this case, non-physical transfer of matter orthogonally to mathbf B and violation of the energy conservation law are possible. To avoid this trouble, we use the direct approximation of Faraday equation (4) with Stokes theorem.
Let's write the approximation for the Faraday equation on base triangular cell depicted in Fig. 1. Consider the unity height triangular prism stretched out of this cell in z direction. For the prism face S based on one of the cell edges we produce the equation (4) by outer normal and integrate it over this face:  Figure 3. Triangular grid containing equilateral triangles.
Using Stokes formula one can modify right hand side of the equation obtained: Since we consider pure 2D model, B z = 0 and we write the following formula for magnetic field normal components on cell edges in new time layer computation: Here B n is the normal component of magnetic field in center of cell edge, τ is time step, is edge length, [·] i z is designation for vertical component of the vector defined in ith vertex of triangular cell. Consider this component in the expression (8) in more detail: The right side of the formula corresponds to the flux component for B y in the system of equations (6), let's denote it as F B y . Note that previously we have computed fluxes F in centers of cells edges. Using these values we compute the average F B y by formula where N being number of edges adjacent to ith vertex. The magnetic field values B inside the grid cells now can be calculated aided by the B n interpolation using Raviart -Thomas vector basis functions [12] with constant divergence: Since the B n are calculated using the Faraday law, summary value of div B is equal to zero in whole the cell.

Testing
The presented algorithm for solving the ideal MHD system is implemented as a program in C++ using the OpenMP library for its parallelization. It is tested and verified on several 2D MHD problems described in detail in the literature [6]. In all calculations the discretization of the computational domain is carried out by equilateral triangles (Fig. 3).

Brio -Wu test
First test is Brio -Wu shock tube [6]. It is the variant of Riemann problem and illustrates process of complex discontinuities evolution. Initial conditions are Boundary conditions (BC) for magnetic field are defined using the field B normal components values in the computational domain boundary edges. Both magnetic and gas BC implement free flow condition. Fig. 4 presents the solution distribution along the x axis. The grid contains n x = 500 triangles along the x axis (along the waves travelling direction). The results show good resolution of If r < r 0 values of gas variables are If r < r 1 = 0.115 values are The plasma between rotating and ambient zones r 0 < r < r 1 has linear distribution of variables: At the initial moment centrifugal forces are not balanced, so the system set in this way is not an equilibrium configuration. The rotating substance will be gradually redistributed over the computational domain, invading its ambient part, while the magnetic field will keep the rotating substance in a flattened form.
We used following values: The solution obtained on a grid with 400 computational cells along one of the axes is shown in Fig. 5. The distribution of density and pressure in the section along the diagonal x = y is shown in Fig. 6.    Fig. 8 shows the pressure distribution along line y = 0.3125. Notice, that finite-difference used doesn't lead to spurious oscillations and has good resolution of discontinuities (only 2 cells on shock wave front).

Conclusions
The numerical scheme for 2D ideal MHD based on physical processes splitting is considered. The unstructured triangular staggered grids technique is used. The matter fluxes (mass, pulse, energy) are calculated using Godunov-type finite-difference scheme and HLLD numerical flux. Evolution of magnetic field is modelled in divergence-free way using direct Faraday law approxiation and Stokes theorem application. The algorythm is implemented in parallel program package. Program is verified using standart test problems. The tests performed showed the effectiveness of the developed technique, which provides suppression of the growth of magnetic field divergence. They also showed that proposed method gives a good resolution for both smooth and discontinuous solutions.