## 1 Introduction

Magnetohydrodynamics (MHD) is the study of the behavior of electrically conducting fluids [1]. This field is important when studying fluids such as plasmas, liquid metals, or electrolytic solutions. Some examples of when the MHD equations would apply is in the formation of stars, behavior of cosmic dust, or plasma confinement in a nuclear fusion reactor. The fundamental equations involved in MHD are a combination of the Navier-Stokes equations from fluid dynamics and Maxwell’s equations from electromagnetism. Note that in this paper we focus only on incompressible MHD. For this we require the Mach number be sufficiently small enough to assume the density of the fluid remains constant. Even given this restriction, incompressible MHD is useful in studying problems such as coolant in a liquid metal fast-fission reactor [2] and the plasma confinement and liquid metal coolant of a fusion reactor [3, 4].

Some popular methods for solving the MHD equations include finite difference methods (FDM), finite volume methods (FVM), and finite element methods (FEM). In the category of FEM, we typically see a stabilized continuous Galerkin (CG) method [5, 6, 7, 8, 9] or a discontinuous Galerkin (DG) method [10, 11, 12, 13]

implemented. DG methods combine the advantages of FVM with FEM. DG methods are ideally suited for solving hyperbolic partial differential equations (PDEs) defined over complex geometries, and they are robust in the presence of large solution gradients, including shocks. The major drawback to DG methods is an increased number of degrees of freedom (DOF) when compared to CG methods.

In order to reduce the number of global DOFs associated with DG methods, hybridized discontinuous Galerkin (HDG) methods were introduced. With an HDG method, there are both interior DOFs that reside on the interior of elements and trace DOFs that reside on the boundaries of elements. With the use of static condensation, the interior DOFs can be written in terms of the trace DOFs, allowing the interior DOFs to be removed completely from the system of equations. This significantly reduces the number of global DOFs. This method was first introduced for symmetric elliptic problems in [14], and it has since been extended to other PDEs such as the advection-diffusion equation [15, 16, 17], the incompressible and compressible Navier-Stokes equations [18, 19], the incompressible Reynolds Averaged Navier-Stokes equations with the Spalart-Allmaras model [20], and the PDEs governing incompressible and compressible MHD [21, 22].

While HDG methods have been developed for the incompressible MHD equations, to the best of the authors’ knowledge, no HDG method presented in the literature to date yields both pointwise divergence-free velocity fields and pointwise divergence-free magnetic fields. However, divergence-conforming methods that exactly preserve divergence-free constraints harbor a number of advantageous properties over methods that satisfy divergence-free constraints in an approximate manner [23, 24]. For instance, CG and DG methods for the incompressible Navier-Stokes equations that yield pointwise divergence-free velocity fields are energy stable and properly balance momentum [25], and they further yield velocity field approximations whose error does not depend on the pressure field [26]

. This latter property is referred to as pressure robustness in the literature. A recent paper further proved velocity field error estimates for such methods that are independent of both the pressure and the Reynolds number

[27]. Divergence-conforming CG and DG methods for the incompressible MHD equations exhibit similar properties [28], though error estimates independent of the Reynolds number and magnetic Reynolds number do not exist yet for such discretizations. The above inspires us to construct an HDG method for the incompressible MHD equations that yields pointwise divergence-free velocity and magnetic fields by extending a divergence-conforming HDG method for the incompressible Navier-Stokes equations [18]. This requires special treatment of the terms that couple the velocity and magnetic fields to arrive at a method that is energy stable. It should be noted that a previously presented HDG method for the incompressible MHD equations does yield pointwise divergence-free velocity fields [29], but it yields only discretely divergence-free magnetic fields.An outline of this paper is as follows. In the next section, we recall the strong form of the incompressible MHD equations. In Section 3, we present our semi-discrete divergence-conforming HDG method for the incompressible MHD equations and show this method is consistent. In Section 4, we prove our method pointwise conserves mass, is pointwise absent of magnetic monopoles, and globally conserves momentum. In Section 5, we prove an energy stability result for our method. In Section 6, we adopt a second-order time integration scheme for our semi-discrete method, and we present a block iterative scheme for solving the nonlinear algebraic equations at each time step. In Section 7, we demonstrate how static condensation may be employed to reduce the size of the linear systems associated with our block iterative scheme. In Section 8, we assess the spatial accuracy, temporal accuracy, pressure robustness, and energy stability of our method using a suite of numerical experiments. In Section 9, we present concluding remarks.

## 2 The Strong Form of the Incompressible MHD Equations

We begin by presenting the strong form of the incompressible MHD equations. Let denote a bounded domain in , where denotes the dimensionality of the domain. We assume or . Let denote the boundary of with outward unit normal , and let be partitioned into a Dirichlet boundary and a Neumann boundary such that and . Let denote the velocity field, denote the pressure field, denote the magnetic field, and denote the magnetic pressure field. We also denote the permeability of free space as . The strong form of the incompressible MHD equations is then as follows:

Find , , , and such that:

Conservation of Momentum

(1) |

Conservation of Mass

(2) |

The Magnetic Induction Equation

(3) |

Gauss’s Law for Magnetism

(4) |

Boundary Conditions

(5) | |||||

(6) | |||||

(7) | |||||

(8) |

Initial Conditions

(9) | ||||

(10) |

In the above equations, denotes the symmetric gradient operator, denotes the antisymmetric gradient operator, denotes the outer product operator, and

denotes the identity matrix of degree

. We take the gradient of a vector

to be where is the standard basis for, and we take the divergence of a second-order tensor

to be . We later adopt the convention for two second-order tensors and . We assume constant density , variable kinematic viscosity , variable resistivity , variable body forces and , variable velocity and magnetic field specifications on the Dirichlet boundary and , variable traction specifications on the Neumann boundary and , and variable initial conditions and . The magnetic pressure is not a physical quantity but is instead included in order to enforce the condition in our HDG method. The magnetic pressure may be interpreted as a Lagrange multiplier associated with the condition and, for a suitable choice of boundary conditions, it is identically zero [30]. Inclusion of a magnetic pressure to enforce the condition is common in both continuous Galerkin [6, 9, 30, 31] and discontinuous Galerkin [21, 29, 32] finite element methods.## 3 A Semi-Discrete Divergence-Conforming HDG Method for the Incompressible MHD Equations

We are now ready to construct our divergence-conforming HDG method for the incompressible MHD equations. In this section, we discretize in space, and later, we discretize in time. As is typically done with HDG methods, we first introduce a mesh over which the incompressible MHD equations will be discretized. Let be a triangulation of the domain such that . The facet of an element is denoted by , and the outward unit normal vector on to is denoted by . Let be the set of all facets. We define the mesh skeleton as . Facets that lie on the boundary of the domain are called boundary facets, while facets that do not lie on the boundary of the domain are called interior facets. We denote the sets of interior and boundary facets as and respectively. Each interior facet is shared by two adjacent elements and . We denote the outward facing normals on to and as and respectively. A visual representation of the above objects in the two-dimensional setting is displayed in Figure 1.

Let denote the space of polynomials of degree on a domain . For , we consider the following finite element spaces for the velocity field

(11) | ||||

(12) |

the pressure field

(13) | ||||

(14) |

the magnetic field

(15) | ||||

(16) |

and the magnetic pressure field

(17) |

We approximate the velocity field, pressure field, magnetic field, and magnetic pressure field over element interiors using the spaces , , , and , and we approximate the velocity field, pressure field, magnetic field, and magnetic pressure field over the mesh skeleton using the spaces , , , and . We also introduce the spaces

(18) | ||||

(19) | ||||

(20) | ||||

(21) |

We employ the spaces and as velocity field trace trial and test spaces, and we employ the spaces and as magnetic field trace trial and test spaces.

As vector-valued functions in the spaces and are discontinuous across element boundaries, we introduce an operator to measure the jump of the normal component of these functions across element boundaries. Let be an interior facet shared by two adjacent elements and with outward facing normals and . For vector-valued functions y lying in either and , we denote the jump of the normal component of y across as where and .

With all of the above notation in hand, we are ready to present our semi-discrete HDG method for the incompressible MHD equations. Our method may be interpreted as an extension of Rhebergen and Wells’s divergence-conforming HDG method for the incompressible Navier-Stokes equations [18], wherein advective fluxes are treated using upwinding and diffusive fluxes using the symmetric interior penalty method, to the incompressible MHD equations. We further take special care in discretizing the Lorentz force appearing in the conservation of momentum equation and the coupling term appearing in the magnetic induction equation to arrive at a method that is energy stable. These two terms are responsible for the transfer of energy between the velocity field and the magnetic field.

Find for each such that:

Conservation of Momentum

(22) |

Conservation of Momentum Flux

(23) |

Conservation of Mass

(24) |

Conservation of Mass Flux

(25) |

The Magnetic Induction Equation

(26) |

Conservation of Magnetic Induction Flux

(27) |

Gauss’s Law of Magnetism

(28) |

Conservation of Magnetic Flux

(29) |

Velocity Initial Condition

(30) |

Magnetic Field Initial Condition

(31) |

In the above equations, on the th facet of the th element , is an indicator function that takes the value of one if the facet is an inflow facet and a value of zero if it is an outflow facet, that is,

(32) |

Moreover, denotes the diameter of element , and is a penalty constant associated with enforcement of the continuity of the velocity and magnetic fields across element boundary boundaries. As is typical with interior penalty methods, the penalty constant must be chosen sufficiently large to ensure the semi-discrete formulation is energy stable [33]. In our later computations, we have selected as we found this choice to be sufficient for energy stability. The above semi-discrete formulation holds for both and , though only two-dimensional numerical examples are presented in this paper. We also did not include the differential quantities and in the above integrands for the sake of conciseness and will continue to do so throughout this paper.

We now show a consistency result for our semi-discrete formulation.

###### Proposition 1 (Consistency)

Proof: Note that, for a sufficiently smooth exact solution , Equations (1)-(6) hold in a pointwise manner. Using this information, we show that (22)-(31) hold if we replace with . Equations (24) and (28) trivially hold since over , Equations (25) and (29) trivially hold since and over , and Equations (30) and (31) trivially hold since and over . Equation (22) holds since

(33) | ||||

for all by reverse integration by parts and

(34) |

in , and Equation (26) holds since

(35) | ||||

for all by reverse integration by parts and

(36) |

in . Finally, Equation (23) holds since

(37) | ||||

for all and

(38) |

on , and Equation (27) holds since

(39) | ||||