The virtual element method for resistive magnetohydrodynamics

https://doi.org/10.1016/j.cma.2021.113815Get rights and content

Abstract

We present a virtual element method (VEM) for the numerical approximation of the electromagnetics subsystem of the resistive magnetohydrodynamics (MHD) model in two spatial dimensions. The major advantages of the virtual element method include great flexibility of polygonal meshes and automatic divergence-free constraint on the magnetic flux field. In this work, we rigorously prove the well-posedness of the method and the solenoidal nature of the discrete magnetic flux field. We also derive stability energy estimates. The design of the method includes three choices for the construction of the nodal mass matrix and criteria to more alternatives. This approach is novel in the VEM literature and allows us to preserve a commuting diagram property. We present a set of numerical experiments that independently validate theoretical results. The numerical experiments include the convergence rate study, energy estimates and verification of the divergence-free condition on the magnetic flux field. All these numerical experiments have been performed on triangular, perturbed quadrilateral and Voronoi meshes. Finally, we demonstrate the development of the VEM method on a numerical model for Hartmann flows as well as in the case of magnetic reconnection.

Introduction

Interest in the behavior of plasmas has skyrocketed in the modern age with applications ranging from fusion-based nuclear power to low power thrusters for contemporary spacecraft. Since the late 1930s, efforts have been devoted to the development of models for plasmas and discretizations that are faithful to the physics and dynamics. An approach that has proven successful and has become standard is to consider plasmas as magnetized fluids, an area called magnetohydrodynamics (MHD). Therefore, the description of these plasmas follow from a blending together of electromagnetic theory and fluid flow. The precise details of how these two theories can be coupled can be found in [1], [2], [3]. Research in MHD is driven by applications that are important to several communities including, astrophysicists that study accretion discs and the dynamics that govern evolution of stars; planetary scientists that are interested in the generation of magnetic fields at the core of planets; plasma physicists whose interest lies in the confinement of plasmas by means of external magnetic fields and engineers who have found that with external magnetic fields they can control the motion of liquid metals leading to a revolution in metallurgical techniques in industry.

The development of numerical methods for MHD is an active area of research, being developed over the last few decades. In [4], [5], two finite element methods are presented that use different techniques in order to preserve the divergence condition on the magnetic field. In [4], the condition is attained automatically, similar to how it is done in this article, whereas in [5] the scheme includes the magnetic vector potential under the temporal gauge, and the magnetic field is obtained as its curl. In [6], the convergence of finite volume methods for MHD is studied and in [7], [8] the classic upwind and Godunov methods are adapted to ideal MHD. In [9], the author presents a finite difference method based on summation by parts (SBP) to mimic the integration by parts formula in the discrete setting, in order to preserve important energy conservation properties and attain an approximate-divergence free scheme. Finally, in [10] the authors develop a MAC scheme for the fluid flow sub-system of the incompressible MHD equations, coupling it to the Yee-scheme for the electromagnetic sub-system.

Although models in MHD come about from a coupling between the equations that govern the fluid flow and Maxwell’s equations for electromagnetism, in this article we will focus on modeling the evolution of the electric and magnetic fields in a plasma for a prescribed fluid flow. Thus, we focus on the Maxwell subsystem of MHD, which combines Faraday’s law, Ampere’s law, Ohm’s law and Gauss’s law for the electromagnetic fields under a prescribed fluid flow.

The main aim of this article is to present a novel numerical discretization of Maxwell’s equations for resistive MHD in a two dimensional setting using the virtual element method (VEM). The VEM was originally proposed in [11] as a variational reformulation of the nodal high-order mimetic finite difference (MFD) method [12], [13], [14], for the numerical treatment of elliptic problems on unstructured polygonal and polyhedral meshes. The word “mimetic” reflects the nature of the method, which mimics the duality and self-adjointness of differential operators as well as identities of vector and tensor calculus. Due to such feature, mimetic methods are often dubbed as compatible methods or compatible discretizations. In particular, satisfying Gauss’s law on the divergence of the magnetic field in the discrete setting requires careful discretization of the Maxwell curl equations, i.e., Faraday’s law and Maxwell–Ampère law. This fact is in contrast to the continuous setting in which the divergence-free nature of the magnetic field is a direct consequence of the Maxwell curl equations when the initial conditions properly satisfy the Gauss’s law.

The violation of Gauss’s law is a serious source of error in the numerical discretization of Maxwell’s equations, causing the appearance of fictitious forces or magnetic monopoles, which are non-physical, thus rendering the numerical simulations unfaithful to the real physics. Over time, mimetic methods were extended from the Support Operator Method (SOM) [15], [16], [17]), which works on regular tensor grids, to the MFD method, which works on fairly general polygonal and polyhedral meshes. The MFD method is, in practice, a family of schemes depending on a set of parameters. These parameters can be optimized to satisfy additional properties such as maximum principles and low dispersion errors. This process goes by the name of mimetic adaptation or M-adaptation and it is outlined in [18]. Previous work in M-adaptation shows that the process can be implemented for problems in wave propagation, see [19], and in the study of cold plasmas, as shown in [20]. Readers interested in historical perspective on the 50-year long development of mimetic and compatible methods are referred to [21]. Development of mimetic and compatible methods are referred to a recent review in [21].

The VEM can also be interpreted as a generalization of the FEM to general polygonal and polyhedral meshes that inherits the great flexibility of the MFD method with respect to the admissible meshes used in the numerical formulation. This versatility provides a series of advantages over methods that require a degree of regularity in their meshes. For example, in order to attain a higher degree of accuracy in a region of the computational domain we tend to use adaptive mesh refinements (AMR). This approach introduces irregularities in the mesh that can be problematic for more traditional numerical methods. We explore this issue in the numerical section, Section 7. In the case of free boundaries or oddly shaped boundary layers it may be imperative to fit the mesh forcing some irregularities. We also note that the numerical dispersion can be greatly reduced on select polygonal meshes, see [22], [23]. In these works, the Finite Difference Time Domain (FDTD) method by Yee was applied to a grid of hexagonal prisms and yields much less numerical dispersion and anisotropy than on the grids that are usually considered in the formulation of Yee’s method, i.e., regular hexahedral cells. We note that another major difference when compared to a regular FEM is that, in the VEM the shape functions are defined in an implicit manner and never explicitly constructed. The name “virtual element” stems from the fact that such shape functions and the finite element space generated by their linear combinations are, in this sense, “virtual”.

The VEM was originally proposed for solving diffusion problems in [11] as a conforming FEM, and later extended to the nonconforming formulation in [24] and the mixed BDM-like and RT-like formulations in [25] and [26], respectively. Generalizations to convection–reaction–diffusion problems with variable coefficients can be found in [27], [28], [29], [30]. In a series of papers [31], [32], [33], [34], H(div)- and H(curl)-conforming virtual element spaces on general polygonal and polyhedral elements have been proposed to generalize the well known Raviart–Thomas and Nédélec finite elements to unstructured polytopal meshes. These methods, combined with the serendipity strategy that reduces the total number of degrees of freedom, see [35], [36], have successfully been applied to the numerical resolution of the magnetostatic Kikuchi’s model. In these papers, exact virtual de Rham sequences with commuting-diagram interpolation operators are built and the solenoidal nature of the discrete magnetic flux field is ensured. Finally, VEMs have also been designed for hyperbolic problems (see  [37], [38]).

In our work, we utilize the low-order spaces proposed in [31], which makes it possible to obtain the combined approximation of the H1-conforming space (0-forms) by a nodal-type virtual element space, the H(curl)-conforming space (1-forms) by an edge-type virtual element space, and the H(div)-conforming space (2-forms) by discontinuous piecewise constant polynomials.

To derive our virtual element approximation, we first reformulate the MHD equations in a variational framework, and then, approximate all L2-type integrals by using suitably defined inner products for nodal-, edge- and cell-type virtual functions. The standard way to build such inner products is through the orthogonal projection of the virtual element functions onto the subspace of linear polynomials. However, as was already noted in [31], the nodal virtual element space that we consider in this work does not provide enough information to construct such projections. Our approach in this paper is to substitute the orthogonal projector with the elliptic projector in [27], since in the low-order case we can always consider these two projection operators as equal by redefining the virtual element space appropriately. This strategy is usually referred to as the “enhanced VEM” by the VEM developers and practitioners.

A major issue occurs here, because changing the definition of the nodal virtual element space requires also changing the definition of the edge and cell virtual element spaces in order to maintain the exact de Rham commuting diagrams. This issue has led to the different virtual element space formulations that were used in the magnetostatics application mentioned above. Instead, in this work we prefer to adopt a different approach, which consists of designing a special reconstruction operator that is computable from the degrees of freedom, stable and bounded as discussed in the following sections. Applying the reconstruction operator makes it possible to recover an approximation of the nodal virtual element functions inside each mesh element and then directly integrate these reconstructed functions. The choice of the elemental reconstruction operator is not unique. In this work, we considered three different options: the elliptic projection; a Least-Squares interpolation of the nodal values; and the piecewise linear Galerkin interpolation on a triangular sub-partition of each element. Our numerical experiments show that these three options are all quite effective and the resulting scheme’s implementations have comparable accuracy.

This article is structured as follows. The rest of this section includes a brief overview of notation and some basic mathematical definitions relevant to the rest of the paper. In Section 2, we present the set of governing equations to be discretized in the continuous setting and introduce the semi-discrete and fully discrete variational formulations in the virtual element framework. Next, in Section 3, we define the virtual element spaces and detail the construction of the inner products that are used for the numerical approximation of the MHD model equations. We also discuss the exactness and commutativity properties of the De-Rham complex and prove that the divergence free condition of the numerical approximation of the magnetic flux field is preserved over time. In Section 4, we prove that the fully discrete variational formulation is well posed. In Section 5, we derive stability energy estimates for the continuous and fully discrete models. In Section 6 we discuss generalizations to arbitrary order and to three dimensions. In Section 7, we present the results of a series of numerical experiments that provide evidence regarding the convergence rate of the numerical method. Plots demonstrating that the method preserves the divergence free condition of the magnetic flux field are available as well as a numerical study of the energy estimates that are derived theoretically in Section 5. This section also includes a simulation of the solution to the Hartmann Flow problem. We conclude the numerical experiments by presenting a model for the simulation of the magnetic reconnection phenomenon. To this end, we use a mesh that is locally refined to provide higher resolution at the points were the phenomenon is taking place. This mesh includes a series of hanging nodes, thus providing an example where using the VEM facilitates the numerical simulation. Then, finally we summarize our findings in Section 8.

We use the standard definition and notation of Sobolev spaces, norms and seminorms, cf. [39]. Let k be a non-negative integer. Consider an open bounded connected subset ω of R2 with polygonal boundary ω. Subset ω can be the whole computational domain Ω, or one of the polygonal cells P of the mesh partitioning, Ωh, covering Ω.

The Sobolev space Hk(ω) consists of all square integrable functions with all square integrable weak derivatives up to order k that are defined on ω. As usual, if k=0, we prefer the notation L2(ω). Norm and seminorm in Hk(ω) are denoted by k,ω and ||k,ω, respectively. We denote the inner product in L2(ω) by (,)ω, but we omit the subscript when ω is the whole computational domain Ω. We denote the norm of an operator Π, which is a norm in the dual space, by the general notation Π, regardless of the spaces where range and image of Π are defined.

On ω, we consider the functional spaces: L2(ω)v:ωR:ω|v|2dV<,H(rot;ω)vL2(ω):rotvL2(ω)2,H(div;ω)wL2(ω)2:divwL2(ω),L(ω){w:ωR:C>0;|w|<Calmost everywhere}, where rotv=(vy,vx)T, and rotw=(wxywyx) for the vector field w=(wx,wy)T. If ω=Ω denotes the computational domain, we consider the functional spaces: VwH(div;Ω):wL2+s(Ω)2,for some real s>0,H0(rot;Ω)vH(rot;Ω):v=0 on Ω. Space V is slightly more regular than H(div;Ω) to ensure that the trace of the normal component vhn|e on each mesh edge e exists and is continuous across all the internal edges [40].

For an open bounded connected subset ωRd with d=1 or 2, we denote the linear space of polynomials of degree up to defined on ω by P(ω), with the useful conventional notation that P1(ω)={0}. We denote the space of two-dimensional vector polynomials of degree up to on ω by [P(ω)]2. Space P(ω) is the span of the finite set of scaled monomials of degree up to , that are given by M(ω)={xxωhωα with |α|},where

  • xω denotes the center of gravity of ω and hω its characteristic length, as, for instance, the edge length or the cell diameter for d=1,2;

  • α=(α1,α2) is the two-dimensional multi-index of nonnegative integers αi with degree |α|=α1+α2 and such that xα=x1α1x2α2 for any xR2.

We will also use the set of scaled monomials of degree exactly equal to , denoted by M(ω) and obtained by setting |α|= in the definition above.

Finally, we use the letter C in many inequalities to denote a strictly positive constant whose value can change at any instance. Constant C may depend on the constants of the model equations or the variational problem, like the coercivity and continuity constants, or constants that are uniformly defined for the family of meshes of the approximation while h0, such as the mesh regularity constant, the stability constants of the discrete bilinear forms, etc. However, constant C will never depend on the discretization parameters such as the mesh size h and the timestep Δt.

Section snippets

The mathematical formulation

Let Ω be an open, bounded, and polygonal subset of R2 with boundary Γ=Ω and T a positive real number. For a given fluid flow described by the velocity vector field u=(ux,uy)T[L(Ω)]2, we consider the Maxwell problem for the electric and magnetic fields, respectively denoted by E and B=(Bx,By)T, that reads as: Bt=rotEνrotBB0E0(,t)withdivB0=0in Ω×(0,T],E+u×B=νrotBrotEB0E0(,t)withdivB0=0in Ω×(0,T],B(,0)=B0withdivB0=0rotEνrotBE0(,t)in Ω,E(,t)=E0(,t)rotEνrotBB0withdivB0=0on Ω×(0,T],

Assumptions on mesh regularity

Let Ωh={P} be a mesh decomposition of Ω into polygonal element (or cell) P with boundary P, area |P| and diameter hP. As usual, h=maxPhhP is the mesh size parameter. We denote the edges of P by e and its length by |e|=he.

We assume that h belongs to H(0,+), which is a countable set of mesh sizes having 0 as its unique accumulation point. A family of meshes {Ωh}h is said to be regular if there exists a non-negative real number ρ independent of h (and, hence, of Ωh), such that

    (M1)

    (

Wellposedness of the virtual element method

Inspired by [4], in this section, we investigate the wellposedness of the virtual element method that we presented in the previous section. The major result of this section is stated by the following theorem.

Theorem 4.1

If θ>0, then solution to Problem 4.2 exists and is unique. Moreover, the map (F,g)(Bhn+1,Êhn+θ) is uniformly continuous independently of h and Δt in the norm defined in Xh.

The definition of the space Xh and its norm will be presented in the next section whereas the proof of this theorem

Stability energy estimates

In this section we show that (6) satisfies an energy estimates. We begin by finding such an estimate for the continuous system (3). The techniques used in the proof are, partially, laid out in [50].

Theorem 5.1

Let B and Ê solve (4) then ddtB0,Ω2+12σ12Ê0,Ω2E0Hσ(rot;Ω)2+2(σ)2u2+1B0,Ω2,where EHσ(rot;Ω)2=σ12E0,Ω2+rotE0,Ω2. As a consequence there exists a bounded function β:[0,T]R+ such that β(t)B(,t)0,Ω2+120tβ(τ)σ12Ê(,τ)0,Ω2dτ0tβ(τ)E0(,τ)Hσ(rot;Ω)2dτ+B0(,t)0,Ω2.

Proof

Testing

Extension of the scheme to higher orders and three dimensions

In this section, we outline how to build high order approximations as well as three dimensional approximations based on results in [31]. We also refer the reader to the papers [32], [33] for other applications in magnetostatics in 2D and 3D.

Let Ω represent a polygonal domain in R2 or a polyhedral domain in R3. Furthermore, consider Ωh to be a mesh of Ω and P a cell in Ωh. In general, for k>1 we define three local spaces Vh,kd(P),Eh,k1d(P) and Ph,k2d(P) where d=2,3 depending on the dimension

Numerical experiments

In this section we will present the results of a series of numerical experiments that sheds some light on the performance of the VEM developed and analyzed throughout this article. It is divided in three sections, the first on explores the rate of convergence and the divergence preserving nature of the numerical method. The second section studies the energy estimate that was introduced in Theorem 5.2. In the final section we introduce the Hartmann problem and use this novel discretization to

Conclusions

In this paper, we have developed the first virtual element method for the time-dependent Maxwell system of Eqs. (3) that model the evolution of the electric and magnetic fields of a magnetized fluid whose flow is prescribed. It is well documented that, in order to accurately describe the physics of resistive MHD, it is imperative for the numerical approximation of the magnetic flux field to remain divergence free. This feature is explicitly addressed in this work and Theorem 4.5 rigorously

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

S. Naranjo Alvarez’s work was supported by the National Science Foundation (NSF) grant #1545188, “NRT-DESE: Risk and uncertainty quantification in marine science and policy”, which provided a one year fellowship and internship support at Los Alamos National Laboratory. In addition, S. Naranjo Alvarez received graduate research funding from V. A. Bokil’s DMS grant #1720116, an INTERN supplemental award to Professor Bokil’s DMS grant # 1720116 for a second internship at Los Alamos National

References (54)

  • BerroneS. et al.

    SUPG stabilization for the nonconforming virtual element method for advection–diffusion–reaction equations

    Comput. Methods Appl. Mech. Engrg.

    (2018)
  • Beirão da VeigaL. et al.

    Virtual element approximation of 2D magnetostatic problems

    Comput. Methods Appl. Mech. Engrg.

    (2017)
  • Beirão da VeigaL. et al.

    Lowest order virtual element approximation of magnetostatic problems

    Comput. Methods Appl. Mech. Engrg.

    (2018)
  • Beirão da VeigaL. et al.

    Serendipity nodal VEM spaces

    Comput. & Fluids

    (2016)
  • VaccaG.

    Virtual element methods for hyperbolic problems on polygonal meshes

    Comput. Math. Appl.

    (2017)
  • AdakD. et al.

    Virtual element method for semilinear Sine-Gordon equation over polygonal mesh using product approximation technique

    Math. Comput. Simulation

    (2020)
  • DassiF. et al.

    Exploring high-order three dimensional virtual elements: bases and stabilizations

    Comput. Math. Appl.

    (2018)
  • ShadidJ.N. et al.

    Scalable implicit incompressible resistive MHD with stabilized FE and fully-coupled Newton–krylov-AMG

    Comput. Methods Appl. Mech. Engrg.

    (2016)
  • DavidsonP.A.

    An introduction to magnetohydrodynamics

    (2002)
  • SchindlerK.

    Physics of space plasma activity

    (2006)
  • MoreauR.J.

    Magnetohydrodynamics, Vol. 3

    (2013)
  • HuK. et al.

    Stable finite element methods preserving B=0 exactly for MHD models

    Numer. Math.

    (2017)
  • HiptmairR. et al.

    A fully divergence-free finite element method for magnetohydrodynamic equations

    Math. Models Methods Appl. Sci.

    (2018)
  • CortiP.

    Stable numerical scheme for the magnetic induction equation with hall effect

  • Beirão da VeigaL. et al.

    Basic principles of virtual element methods

    Math. Models Methods Appl. Sci.

    (2013)
  • Beirão da VeigaL. et al.

    Post-processing of solution and flux for the nodal mimetic finite difference method

    Numer. Methods Partial Differential Equations

    (2015)
  • Beirão da VeigaL. et al.

    Arbitrary order nodal mimetic discretizations of elliptic problems on polygonal meshes

    SIAM J. Numer. Anal.

    (2011)
  • Cited by (14)

    • The virtual element method for a 2D incompressible MHD system

      2023, Mathematics and Computers in Simulation
    • Conforming virtual element approximations of the two-dimensional Stokes problem

      2022, Applied Numerical Mathematics
      Citation Excerpt :

      The connection with de Rham diagrams and Nedelec elements and the application to the electromagnetics has been explored in [10]. A practical application of these concepts can be found in [12,77]. Other significant applications of the VEM on general meshes are found, for example, in [4,5,20–23,25–31,37,40,41,45–47,60,76,78–80,90,91] and [58,86–88].

    • Virtual elements for Maxwell's equations

      2022, Computers and Mathematics with Applications
    View all citing articles on Scopus
    View full text