A convergent linearized Lagrange finite element method for the magneto-hydrodynamic equations in 2D nonsmooth and nonconvex domains

TL;DR

This paper introduces a new linearized finite element method for 2D magneto-hydrodynamics equations that converges reliably in complex, nonconvex domains without requiring high regularity assumptions.

Contribution

It develops a fully discrete, linearized $H^1$-conforming Lagrange finite element method that guarantees convergence in general domains with minimal regularity assumptions.

Findings

Proves convergence of numerical solutions in nonconvex, nonsmooth domains.

Establishes strong convergence in $L^2$-norms for velocity and magnetic field.

No mesh restrictions needed for convergence.

Abstract

A new fully discrete linearized $H^1$-conforming Lagrange finite element method is proposed for solving the two-dimensional magneto-hydrodynamics equations based on a magnetic potential formulation. The proposed method yields numerical solutions that converge in general domains that may be nonconvex, nonsmooth and multi-connected. The convergence of subsequences of the numerical solutions is proved only based on the regularity of the initial conditions and source terms, without extra assumptions on the regularity of the solution. Strong convergence in $L^2(0,T;{\bf L}^2(\Omega))$ was proved for the numerical solutions of both ${\bm u}$ and ${\bm H}$ without any mesh restriction.

Tables

Table 1. Table 1: Numerical results of the proposed numerical method ( 38 )-( 42 ) with τ = h 𝜏 ℎ \tau=h

$h$	${\Vert {𝑯}_h^N-𝑯(\cdot ,t_N)\Vert }_{L^2}$	${\Vert {𝒖}_h^N-𝒖(\cdot ,t_N)\Vert }_{L^2}$
1/16	4.295E-03	6.540E-05
1/32	1.894E-03	2.073E-05
1/64	1.092E-03	8.492E-06
1/128	7.036E-04	3.845E-06
convergence rate	$O\left(h^{0.63}\right)$	$O\left(h^{1.14}\right)$

Table 2. Table 2: Numerical results of the H 1 superscript 𝐻 1 H^{1} -conforming FEM ( 5.1 )-( 128 ) with τ = h 𝜏 ℎ \tau=h

$h$	${\Vert {𝑯}_h^N-𝑯(\cdot ,t_N)\Vert }_{L^2}$	${\Vert {𝒖}_h^N-𝒖(\cdot ,t_N)\Vert }_{L^2}$
1/16	1.180E-02	6.516E-05
1/32	7.419E-03	2.046E-05
1/64	4.958E-03	8.297E-06
1/128	3.206E-03	3.723E-06
convergence rate	$O\left(h^{0.63}\right)$	$O\left(h^{1.16}\right)$

Table 3. Table 3: Numerical results of the proposed numerical method ( 38 )-( 42 ) with τ = h 𝜏 ℎ \tau=h

$h$	${\Vert {𝑯}_h^N-𝑯(\cdot ,t_N)\Vert }_{L^2}$	${\Vert {𝒖}_h^N-𝒖(\cdot ,t_N)\Vert }_{L^2}$
1/16	5.625E-02	6.376E-05
1/32	3.472E-02	2.232E-05
1/64	2.179E-02	8.678E-06
1/128	1.366E-02	3.767E-06
convergence rate	$O\left(h^{0.67}\right)$	$O\left(h^{1.20}\right)$

Table 4. Table 4: Numerical results of the H 1 superscript 𝐻 1 H^{1} -conforming FEM ( 5.1 )-( 128 ) with τ = h 𝜏 ℎ \tau=h

$h$	${\Vert {𝑯}_h^N-𝑯(\cdot ,t_N)\Vert }_{L^2}$	${\Vert {𝒖}_h^N-𝒖(\cdot ,t_N)\Vert }_{L^2}$
1/16	4.681E-01	1.290E-04
1/32	2.912E-01	4.582E-05
1/64	1.819E-01	1.798E-05
1/128	1.153E-01	7.388E-06
convergence rate	$O\left(h^{0.66}\right)$	$O\left(h^{1.28}\right)$

Equations

$$\mu {\partial }_t\mathbf{H}+{\sigma }^{-1}\nabla \times (\nabla \times \mathbf{H})-\mu \nabla \times (\mathbf{u}\times \mathbf{H})={\sigma }^{-1}\nabla \times J$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}-\nu \Delta \mathbf{u}+\nabla p=\mathbf{f}-\mu \mathbf{H}\times (\nabla \times \mathbf{H})$$

$$\nabla \cdot \mathbf{u}=0$$

$$\nabla \times \mathbf{B}=\frac{\partial B_2}{\partial x_1}-\frac{\partial B_1}{\partial x_2}\nabla \cdot \mathbf{B}=\frac{\partial B_1}{\partial x_1}+\frac{\partial B_2}{\partial x_2}$$

$$\displaystyle\nabla\times\psi=\bigg{(}\frac{\partial\psi}{\partial x_{2}}\,,\,-\frac{\partial\psi}{\partial x_{1}}\bigg{)}\quad\nabla\psi=\bigg{(}\frac{\partial\psi}{\partial x_{1}}\,,\,\frac{\partial\psi}{\partial x_{2}}\bigg{)}$$

$$\mathbf{B}\times \psi =(B_2\psi ,-B_1\psi ).$$

$$\mathbf{H}\cdot \mathbf{n}=0,\nabla \times \mathbf{H}=J,\text{mbox}and\mathbf{u}=0\text{mbox}on\partial \Omega \times (0,T],$$

$$\mathbf{H}{|}_{t=0}={\mathbf{H}}_0\text{mbox}and\mathbf{u}{|}_{t=0}=u_0\text{mbox}in\Omega .$$

$$\begin{array}{rl}& J\in C([0,T];L^2(\Omega ))\mathbf{f}\in C([0,T];{\mathbf{L}}^2(\Omega ))\\ & {\mathbf{H}}_0,{\mathbf{u}}_0\in {\mathbf{L}}^2(\Omega )\nabla \cdot {\mathbf{H}}_0=\nabla \cdot {\mathbf{u}}_0=0.\end{array}$$

$$({\bm{H}},{\bm{u}})\in L^{2}(0,T;{\bf H}^{1}(\Omega))^{2}\cap L^{\infty}(0,T;{\bf L}^{2}(\Omega))^{2}\quad\mbox{such that $\nabla\cdot{\bm{H}}=\nabla\cdot{\bm{u}}=0$}.$$

$$(\mathbf{u},p)\in {\mathbf{H}}^{s+1}(\Omega )\times H^s(\Omega )\text{mbox}and(\mathbf{H},\nabla \times \mathbf{H})\in {\mathbf{H}}^s(\Omega )\times H^s(\Omega )\text{mbox}forsomes>\frac{1}{2}$$

$$\mu {\partial }_{t}A-{\sigma }^{-1}\Delta A-\mu \mathbf{u}\times (\nabla \times A)={\sigma }^{-1}J$$

$$\displaystyle\partial_{t}{\bm{u}}+{\bm{u}}\cdot\nabla{\bm{u}}-\nu\Delta{\bm{u}}+\nabla p={\bm{f}}+\mu\nabla\cdot\big{[}(\nabla\times A)\otimes(\nabla\times A)-\frac{1}{2}|\nabla\times A|^{2}\mathbb{I}\big{]}$$

$$\nabla \cdot \mathbf{u}=0$$

$${\partial }_t\mathbf{H}=-\nabla \times \mathbf{E}$$

$$\mathbf{E}=-\mathbf{u}\times \mathbf{H}+{\sigma }^{-1}\nabla \times \mathbf{H}.$$

$$\sigma {\partial }_t\mathbf{A}+\sigma (\nabla \times \mathbf{A})\times \mathbf{u}+{\mu }^{-1}\nabla \times \nabla \times \mathbf{A}=0$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}-\nu \Delta \mathbf{u}+\nabla p=\mathbf{f}-\sigma ({\partial }_t\mathbf{A}+(\nabla \times \mathbf{A})\times \mathbf{u})\times \nabla \times \mathbf{A}$$

$$\nabla \cdot \mathbf{u}=0$$

$$\underset{k\to \infty }{\lim }\Vert {\mathbf{u}}_{\tau ,h_k}^n-{\mathbf{u}}_{\tau }^n{\Vert }_{L^4(\Omega )}\to 0$$

$$\displaystyle\,\,{\bm{A}}_{\tau,h_{k}}^{n}\rightharpoonup{\bm{A}}_{\tau}^{n}\quad\mbox{weakly in $L^{\frac{3}{2}}(\Omega)$ as $k\rightarrow\infty$}$$

$$\underset{k\to \infty }{\lim }\Vert {\mathbf{H}}_{\tau ,h_k}^n-{\mathbf{H}}_{\tau }^n{\Vert }_{L^2(\Omega )}\to 0$$

$$\nabla \cdot \mathbf{H}=0.$$

$$\mathbf{H}=\nabla \times A+\sum _{j=1}^m{\beta }_j\nabla \times {\phi }_j,$$

$$\{\begin{array}{rlrl}& -\Delta A=\nabla \times \mathbf{H}& & \text{mbox}in\Omega \\ & A=0& & \text{mbox}on\partial \Omega ,\end{array}$$

$$\{\begin{array}{rlrl}& \Delta {\phi }_j=0& & \text{mbox}in\Omega \\ & {\phi }_j=1& & \text{mbox}on{\Gamma }_j\\ & {\phi }_j=0& & \text{mbox}on\partial \Omega \backslash {\Gamma }_j.\end{array}$$

$$\nabla\times\nabla\times\phi=\nabla\times\bigg{(}\frac{\partial\phi}{\partial x_{2}}\,,\,-\frac{\partial\phi}{\partial x_{1}}\bigg{)}=\frac{\partial}{\partial x_{1}}\bigg{(}-\frac{\partial\phi}{\partial x_{1}}\bigg{)}-\frac{\partial}{\partial x_{2}}\bigg{(}\frac{\partial\phi}{\partial x_{2}}\bigg{)}=-\Delta\phi.$$

$$({\partial }_t\mathbf{H},\nabla \times {\phi }_i)=\sum _{j=1}^m\frac{\mathrm{d}{\beta }_j}{\mathrm{d}t}(\nabla \times {\phi }_j,\nabla \times {\phi }_i),$$

$$(\nabla \times A,\nabla \times {\phi }_i)=-(A,\nabla \times (\nabla \times {\phi }_i))=(A,\Delta {\phi }_i)=0.$$

$$\displaystyle(\partial_{t}{\bm{H}},\nabla\times\varphi_{i})=0\quad\mbox{for $i=1,\dots,m$}.$$

$$\frac{\mathrm{d}{\beta }_j}{\mathrm{d}t}=0\text{mbox}fori=1,\dots ,m.$$

$$\displaystyle\nabla\times\big{(}\mu\partial_{t}A-\sigma^{-1}\Delta A-\sigma^{-1}J-\mu{\bm{u}}\times(\nabla\times A+\mbox{$\sum_{j=1}^{m}$}\beta_{j}\nabla\times\varphi_{j})\big{)}={\bm{0}}$$

$$\displaystyle\mu\partial_{t}A-\sigma^{-1}\Delta A-\sigma^{-1}J-\mu{\bm{u}}\times\big{(}\nabla\times A+\mbox{$\sum_{j=1}^{m}$}\beta_{j}\nabla\times\varphi_{j}\big{)}={\rm const}.$$

$$\displaystyle\mu\partial_{t}A-\sigma^{-1}\Delta A-\mu{\bm{u}}\times\big{(}\nabla\times A+\mbox{$\sum_{j=1}^{m}$}\beta_{j}\nabla\times\varphi_{j}\big{)}=\sigma^{-1}J$$

$$\displaystyle\partial_{t}{\bm{u}}+{\bm{u}}\cdot\nabla{\bm{u}}-\nu\Delta{\bm{u}}+\nabla p={\bm{f}}+\mu\big{(}\nabla\times A+\mbox{$\sum_{j=1}^{m}$}\beta_{j}\nabla\times\varphi_{j}\big{)}\times\Delta A$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\sum _{j=1}^m{\beta }_j(\nabla \times {\phi }_j,\nabla \times {\phi }_i)=({\mathbf{H}}_0,\nabla \times {\phi }_i)$$

Full text

1 Introduction

This article is concerned with numerical approximation of the incompressible magneto-hydrodynamics (MHD) equations

[TABLE]

in a two-dimensional polygonal type domain $\Omega=\Omega_{0}\backslash(\cup_{j=1}^{m}\Omega_{j})\subset{\mathbb{R}}^{2}$, where both $\Omega_{0}$ and $\Omega_{j}\subset\Omega_{0}$, $j=1,\dots,m$, are polygons (thus the domain is possibly nonconvex and multi-connected), and the following 2D notations for the curl, divergence, and gradient operators are used for a vector field ${\bm{B}}=(B_{1},B_{2})$ and a scalar field $\psi$:

[TABLE]

In the system (1)-(3), ${\bm{u}}$ denotes the velocity field, ${\bm{H}}$ the magnetic field, $p$ the pressure, ${\bm{f}}$ and $J$ the given source terms, $\nu$ the viscosity of the fluid, $\sigma$ the magnetic Reynolds number, and $\mu=M^{2}\nu\sigma^{-1}$, where $M$ denotes the Hartman number.

We consider (1)-(3) under the perfectly conducting and no-slip boundary conditions

[TABLE]

and the initial conditions

[TABLE]

The given source terms and the initial data for ${\bm{H}}$ and ${\bm{u}}$ are assumed to satisfy

[TABLE]

Under the assumption $\nabla\cdot{\bm{H}}_{0}=0$ for the initial value ${\bm{H}}_{0}$, a solution of (1) will automatically satisfy $\nabla\cdot{\bm{H}}=0$ for all time.

The MHD equations (1)-(3) describe the interaction between a magnetic field and a viscous incompressible conducting fluid flow. The mathematical theory of existence and uniqueness of weak solutions for the initial-boundary value problem (1)-(5) was established in [40, Theorem 3.1] for a smooth or convex domain $\Omega$, where the weak formulation of the equations (1)-(5) does not involve the variable $p$; see [40, Problem 2.1]. In particular, under the regularity assumption (6) for the source terms and initial data, the problem has a unique weak solution (for any given $T>0$)

[TABLE]

Numerical methods and analysis for the MHD equations have been done from many different point of views. For the stationary MHD equations, existence, uniqueness, and finite element approximations were studied in [21] for small data. To overcome the numerical instability caused by possibly small hydrodynamic diffusion, a stabilized finite element method (FEM) was introduced in [16]. These articles are concerned with $H^{1}$-conforming FEMs for the magnetic field ${\bm{H}}$ and the error estimates are based on the ${\bf H}^{1}(\Omega)$-regularity of the magnetic field ${\bm{H}}$. Such regularity holds in convex or smooth domains. However, in nonconvex and nonsmooth domains, the solution of the magnetic field is generally in ${\bf H}({\rm curl},\Omega)$ instead of ${\bf H}^{1}(\Omega)$.

In more general domains, possibly nonconvex and nonsmooth, a mixed FEM with curl-conforming Nédélec edge elements was proposed for solving the magnetic field in the stationary MHD equations in [39], where an additional gradient term $\nabla q$ was added to the magnetic potential equation to enforce the divergence-free condition for the magnetic field in a weak sense; the $H^{1}$-conforming FEM was used for the velocity. An error estimate for this numerical method was proved under the regularity assumptions

[TABLE]

which hold when the source term of the stationary MHD equations is sufficiently small. The same method for the magnetic field was also used in [18] for solving the MHD equations, where a divergence-conforming FEM was used for the velocity.

In the case of low magnetic Reynolds numbers, the MHD model usually consists of a time-dependent Navier-Stokes equation and a stationary electric potential equation with given magnetic field; see [21, 34]. In [27], two implicit-explicit methods (of first and second order) decoupling the velocity from electric potential were proposed and analyzed. In [37], the method proposed in [27] was further combined with the technique of [42] to decouple the pressure from velocity. For the model with low magnetic Reynolds numbers considered in [27, 37], the $H^{1}$-conforming FEMs were proved to be convergent.

For time-dependent MHD equations, many different numerical methods have been developed and analyzed:

1. $\bullet$

$H^{1}$-conforming FEMs for the velocity and magnetic fields were used in [2, 19, 23, 24, 45] with several different time discretization methods. These methods work very well in smooth or convex domains. In this type of domains, error estimates have been established for several of these $H^{1}$-conforming FEMs. Similarly as for the stationary problem, in a nonconvex and nonsmooth domain, the solution of the magnetic field of the time-dependent MHD equations is generally in ${\bf H}({\rm curl},\Omega)$ in the spatial direction instead of ${\bf H}^{1}(\Omega)$. In this case the ${\bf H}^{1}(\Omega)$ or ${\bf H}^{2}(\Omega)$ regularity assumptions used in the existing analyses for the $H^{1}$-conforming FEMs generally do not hold. 2. $\bullet$

A Galerkin least square FEM was proposed for solving the augmented MHD equations by adding an additional gradient term $\nabla q$ to the equation of magnetic field [38], where numerical simulations were shown to illustrate the performance of the numerical methods. Rigorous proof for the convergence of numerical solutions remains open. 3. $\bullet$

The magnetic potential formulation was used in [41], where the equivalent formulation of the two-dimensional MHD equations

[TABLE]

was solved by a fully implicit $H^{1}$-conforming FEM, where $\mathbb{I}$ denotes the identity matrix. Numerical results were given without proof for the convergence of numerical solutions. 4. $\bullet$

Divergence-free preserving methods for MHD and ideal MHD have been developed in many articles. In particular, the locally divergence-free subspace was used in the discontinuous Galerkin (DG) methods for the MHD equations [28]. The divergence-free subspace of the Brezzi-Douglas-Marini (BDM) finite element space was used for the magnetic field in [29, 30]. An additional bubble function was added to each element in [8] in order to have additional degree of freedoms to enforce the divergence-free condition. More recently, the equation of magnetic field was rewritten as a first-order system in [26], i.e.,

[TABLE]

Then, the divergence-conforming and curl-conforming Nédélec edge elements were used for ${\bm{H}}$ and ${\bm{E}}$, respectively. In this method no constraint was enforced on the magnetic field, but the numerical solution automatically satisfied the divergence-free condition provided that the initial data for the magnetic field was divergence free. For these methods, numerical results were shown to illustrate the performance of the numerical methods, while convergence proofs remain open. 5. $\bullet$

A finite difference method for the 2D incompressible MHD equations in rectangular domains was proposed in [31], where the MAC and Yee’s scheme were used for the fluid and magnetic equations, respectively. The convergence of this method was proved in [36] by using discrete energy estimates in $L^{\infty}(0,T;H^{1}(\Omega))$ and $L^{2}(0,T;H^{2}(\Omega))$ for both velocity and magnetic fields. An energy and helicity preserving finite difference method was proposed in [32] for the 3D axisymmetric incompressible MHD equations. The convergence of this method was demonstrated numerically but remains open theoretically. An error estimate of the corresponding energy and helicity preserving method for a single fluid equation was proved in [33] based on some regularity assumption. 6. $\bullet$

As far as we know, the only existing proofs for the (subsequence) convergence of numerical solutions in possibly nonsmooth domains were given in [35] and [25]. In [35], Prohl studied several fully discrete linearized FEMs (with different time discretization and decoupling methods) with curl-conforming Nédélec edge elements for the magnetic field, and proved the convergence of two numerical schemes to weak solutions under the mesh restrictions $\tau=O(h^{4})$ and $\tau=O(h^{3})$, respectively, where $\tau$ denotes the time-step size and $h$ the spacial mesh size. Without such mesh restrictions, the weak∗ convergence of numerical solutions in $L^{\infty}(0,T;{\bf L}^{2}(\Omega))$ was proved. In [25], Hiptmair, Li, Mao and Zheng discretized a magnetic potential formulation of the three-dimensional MHD equations:

[TABLE]

which is different from the magnetic potential formulation in [41]. The curl-conforming Nédélec elements were used for the discretization of ${\bm{A}}$. For every fixed time-step size $\tau$, it was proved that a subsequence $({\bm{u}}_{\tau,h_{k}}^{n},{\bm{H}}_{\tau,h_{k}}^{n})$, $k=1,2,\dots$, of numerical solutions will converge to a semi-discrete solution $({\bm{u}}_{\tau}^{n},{\bm{H}}_{\tau}^{n})$, i.e.,

[TABLE]

It was also shown that a subsequence $({\bm{u}}_{\tau_{m}},{\bm{A}}_{\tau_{m}})$, $m=1,2,\dots$, of the semi-discrete solutions will converge weakly in $L^{2}(0,T;{\bf H}^{1}_{0}(\Omega))\times W^{1,p}(0,T;{\bf H}({\rm curl}))$ for some $p>1$ to a weak solution of the MHD equations. Strong convergence in $L^{2}(0,T;{\bf L}^{2}(\Omega))$ of the numerical solutions of ${\bm{u}}$ and ${\bf H}$ remains open.

Overall, the existing $H^{1}$-conforming FEMs for the MHD equations were proved to be convergent only in convex or smooth domains. In more general domains the mixed FEM with curl-conforming Nédélec elements is more suitable for approximating the magnetic field directly, while existing proofs for the convergence of numerical solutions either require mesh restriction $\tau=O(h^{3})$ or yield only weak- or weak∗-convergence of the numerical solutions.

The most important applications of the MHD model considered in the present paper occur in metallurgy and liquid-metal processing [3, 43]. In that context, consideration of general polygonal yet nonconvex and multi-connected domains is very relevant. The aim of this article is to develop a new fully discrete linearized $H^{1}$-conforming Lagrange FEM for the two-dimensional MHD equations based on a magnetic potential formulation such that the numerical solutions would converge not only in convex and smooth domains but also in nonconvex and nonsmooth domains.

In the magnetic potential formulation, the magnetic potential $A$ would naturally has $H^{1}$ regularity and therefore can be solved correctly by using $H^{1}$-conforming finite element methods. Then the magnetic field is obtained through taking partial derivatives, i.e., ${\bm{H}}=\nabla\times A+\sum_{j=1}^{m}\beta_{j}\nabla\times\varphi_{j}$, where $\beta_{j}$ and $\varphi_{j}$ are some time-independent constants and functions, depending on the geometry of the computational domain. This approach is different from the singular basis methods [4, 22], weighted regularization method [9], $L^{2}$ projected method [14] and the residual-based stabilized augmented formulation [5]. These methods focus on approximating the magnetic field ${\bm{H}}$ in the Maxwell equations directly by using $H^{1}$-conforming finite element methods. Similarly as [25, 35], the proof of convergence in this paper is only based on the regularity of the initial conditions and source terms, without any extra assumptions on the regularity of the solution. Strong convergence of subsequences in $L^{2}(0,T;{\bf L}^{2}(\Omega))$ as $(\tau_{n},h_{n})\rightarrow(0,0)$ is proved for the numerical solutions of ${\bm{u}}$ and ${\bm{H}}$ without mesh restrictions.

The rest of the paper is organized as follows. In Section 2, we present an equivalent magnetic potential formulation of the the two-dimensional MHD equations (1)-(3). In Section 3, we propose a fully discrete linearized $H^{1}$-conforming Lagrange finite element method for solving the problem, and present the main theoretical result about the convergence the numerical solutions. Rigorous proof of the main theoretical result is presented in Section 4. Numerical experiments are given in Section 5 to support the theoretical analyses.

2 Equivalent formulation

In section 2.1 we first formally derive an equivalent formulation of the two-dimensional MHD equations (1)-(3) in terms of the magnetic potential. Then we define weak solutions of the problem in section 2.2. It is easy to verify that a weak solution of the reformulated problem is also a weak solution of the original MHD equations (see Remark 2.1).

2.1 Formal derivation

By taking the divergence of (1) we obtain $\mu\partial_{t}\nabla\cdot{\bm{H}}=0,$ which together with the divergence-free initial condition $\nabla\cdot{\bm{H}}_{0}=0$ in (6) implies

[TABLE]

Let $m$ denote the number of holes of the domain $\Omega$, and let $\Gamma_{j}=\partial\Omega_{j}$ denote the boundary of the $j$th hole. Then the divergence-free vector field ${\bm{H}}$ can be decomposed as (cf. [7], with slightly different boundary conditions)

[TABLE]

where $\beta_{j}$ for $j=1,\dots,m$, are constants independent of the spatial variables, $A$ is the solution of

[TABLE]

and $\varphi_{j}$ is the solution of

[TABLE]

Note that for the scalar-valued function $\phi$, there holds

[TABLE]

Therefore, integrating the time derivative of (15) against $\nabla\times\varphi_{i}$ yields

[TABLE]

where we have used (16)-(17) and

[TABLE]

Furthermore, integrating (1) against $\nabla\times\varphi_{i}$ and using integration by parts, and with the boundary condition (4), we obtain

[TABLE]

Since the matrix $M_{ij}=(\nabla\times\varphi_{j},\nabla\times\varphi_{i})$ is positive definite, the two identities above imply

[TABLE]

Therefore, $\beta_{j}$ for $j=1,\dots,m$, are constants independent of time.

Now we substitute ${\bm{H}}=\nabla\times A+\sum_{j=1}^{m}\beta_{j}\nabla\times\varphi_{j}$ into (1). This yields

[TABLE]

which implies

[TABLE]

With the boundary condition (4), it is easy to derive that the constant on the right-hand side of the above equation equals zero.

Thus, instead of solving (1)-(2) directly, we propose to solve (17) and the following equations:

[TABLE]

where $\beta_{j}$, $j=1,\dots,m$, are determined by

[TABLE]

which can be obtained by integrating (15) against $\nabla\times\varphi_{i}$ at time $t=0$.

The boundary and initial conditions for (20)-(22) are given by

[TABLE]

and

[TABLE]

where $A_{0}$ is the solution of

[TABLE]

After solving of (17) and (20)-(25), we can obtain the magnetic field {\bm{H}}=\nabla\times A+\mbox{\sum_{j=1}^{m}}\beta_{j}\nabla\times\varphi_{j}.

2.2 Weak solution

For $k\geq 0$ and $p\in[1,\infty]$, we denote by $W^{k,p}(\Omega)$ the conventional Sobolev space of functions defined on $\Omega$, with abbreviations $L^{p}(\Omega)=W^{0,p}(\Omega)$ and $H^{k}(\Omega)=W^{k,2}(\Omega)$. Let $W^{1,p}_{0}(\Omega)$ be the space of functions in $W^{1,p}(\Omega)$ with zero traces on the boundary $\partial\Omega$, and denote $H^{1}_{0}(\Omega)=W^{1,2}_{0}(\Omega)$. The corresponding vector-valued spaces are

[TABLE]

The inner product of $L^{2}(\Omega)$ is denoted by $(\cdot,\cdot)$.

A quadruple $(A,{\bm{u}},(\varphi_{j})_{j=1}^{m},(\beta_{j})_{j=1}^{m})$ is called a weak solution of (17) and (20)-(26) if

[TABLE]

with $A|_{t=0}=A_{0}$ and ${\bm{u}}|_{t=0}={\bm{u}}_{0}$, and the following equations hold for all test functions $a\in L^{\infty}(0,T;H_{0}^{1}(\Omega))$ and ${\bm{v}}\in L^{\infty}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$,

[TABLE]

Remark 2.1.

1. (1)

The pressure $p$ does not appear in the definition of weak solutions as we have restricted both ${\bm{u}}$ and the test function ${\bm{v}}$ to ${\bf H}^{1}_{0,{\rm div}}(\Omega)$. 2. (2)

If the domain is simply connected (without any holes) then $m=0$. In this case, a pair $(A,{\bm{u}})$ is called a weak solution if (27)-(28) and (31)-(32) hold. 3. (3)

By substituting $a=\nabla\times{\bm{b}}$ into (31), we see that if $(A,{\bm{u}},(\varphi_{j})_{j=1}^{m},(\beta_{j})_{j=1}^{m})$ is a weak solution of (17) and (20)-(26) then $({\bf H},{\bm{u}})$ is a weak solution of (1)-(5) in the sense that

[TABLE]

and the following equations hold for all test functions ${\bm{b}}\in L^{\infty}(0,T;{\bf H}(\Omega;{\rm curl}))$ and ${\bm{v}}\in L^{\infty}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$,

[TABLE]

3 Numerical method

In this section, we introduce a fully discrete numerical method for solving (17) and (20)-(26), and then present the main theoretical results about the convergence of the numerical solutions.

Let $\Im_{h}$, $0<h<1$, be a family of quasi-uniform partition of $\Omega$ into triangles $K_{j}$, $j=1,\ldots,M$, such that no vertex of any triangle lies on the interior of an edge of another triangle. Here, quasi-uniformity means that the diameter $h_{K_{j}}$ of a triangle $K_{j}$ and the radius $\rho_{K_{j}}$ of the inscribed circle satisfy

[TABLE]

where $\displaystyle h=\max_{1\leq j\leq M}\{\mbox{diam}(K_{j})\}$ is the mesh size, and $c$ is a constant independent of $h$.

For any integer $r\geq 1$, we define the Taylor–Hood finite element space $\mathring{\bf S}_{h}^{r+1}\times S_{h}^{r}/{\mathbb{R}}$ with

[TABLE]

where $P_{r}(K_{j})$ is the space of polynomials of degree $r$ on the triangle $K_{j}$.

For any given $j=1,\dots,m$, let $\varphi_{j,h}\in S_{h}^{r+1}$ be the finite element solution of (17), i.e.,

[TABLE]

such that $\varphi_{j,h}=1$ on $\Gamma_{j}$ and $\varphi_{j,h}=0$ on $\partial\Omega\backslash\Gamma_{j}$. Let $\beta_{j,h}$, $j=1,\dots,m$, be the constants (independent of space and time) determined by the equations

[TABLE]

Let $\{t_{n}=n\tau\}_{n=0}^{N}$ denote a uniform partition of the time interval $[0,T]$, with a step size $\tau=T/N$, and $u^{n}=u(x,t_{n})$. A fully discrete numerical scheme for the system (20)-(22) is to find $A_{h}^{n}\in\mathring{S}_{h}^{r+1}$, ${\bm{u}}_{h}^{n}\in\mathring{\bf S}_{h}^{r+1}$, and $p_{h}^{n}\in S_{h}^{r}/{\mathbb{R}}$ such that

[TABLE]

hold for all test functions $a_{h}\in\mathring{S}_{h}^{r+1}$, ${\bm{v}}_{h}\in\mathring{\bf S}_{h}^{r+1}$, $q_{h}\in S_{h}^{r}/{\mathbb{R}}$, and $n=1,2,\dots,N$.

Then the discrete magnetic field can be obtained by

[TABLE]

As a result, ${\bf H}_{h}^{n}$ satisfies $\nabla\cdot{\bf H}_{h}^{n}=0$ automatically.

Here, the operator $\Delta_{h}:\mathring{S}_{h}^{r+1}\rightarrow\mathring{S}_{h}^{r+1}$ is defined via the duality:

[TABLE]

The initial condition $A^{0}_{h}\in S_{h}^{r+1}$ can be determined by

[TABLE]

The initial condition for velocity is given by ${\bm{u}}_{h}^{0}=P_{h}{\bm{u}}_{0}$, where $P_{h}$ denotes the $L^{2}$ projection from $L^{2}(\Omega)$ onto $\mathring{S}_{h}^{r+1}$, defined by

[TABLE]

which automatically extends to the vector-valued space $\mathring{\bf S}_{h}^{r+1}$, i.e.,

[TABLE]

For any sequence $\omega_{h}^{n}$, $n=1,2,\dots,$ we define the piecewise constant functions $\omega_{h,\tau}^{+}$ and $\omega_{h,\tau}^{-}$ by

[TABLE]

for $t\in(t_{n-1},t_{n}]$ and $n=1,2,\dots,N$. Then we have the following result on the convergence of the numerical solutions.

Theorem 1**.**

Under assumption (6), the fully discrete finite element scheme (38)-(42) has a unique solution. For any sequence $(\tau_{n},h_{n})\rightarrow(0,0)$, there exists a subsequence $(\tau_{n_{k}},h_{n_{k}})\rightarrow(0,0)$ such that the corresponding numerical solutions converge to a weak solution $(A,{\bm{u}},(\varphi_{j})_{j=1}^{m},(\beta_{j})_{j=1}^{m})$ of (17) and (20)-(26) in the following sense:

[TABLE]

which also imply that {\bm{H}}_{h_{n_{k}},\tau_{n_{k}}}^{+}=\nabla\times A_{h_{n_{k}},\tau_{n_{k}}}^{+}+\mbox{\sum_{j=1}^{m}}\beta_{j,h}\nabla\times\varphi_{j,h} converges to {\bm{H}}=\nabla\times A+\mbox{\sum_{j=1}^{m}}\beta_{j}\nabla\times\varphi_{j} in $L^{2}(0,T;{\bf L}^{2}(\Omega))$.

Remark 3.1**.**

The uniqueness of weak solutions was proved in [40] for convex and smooth domains, but remains open for nonconvex and nonsmooth domains. Our proof shows that every sequence of numerical solutions contains a subsequence converging to a weak solution of the PDE problem. If the weak solution is unique, then Theorem 1 implies that the numerical solutions converge to the unique weak solution as $(\tau,h)\rightarrow(0,0)$ (without taking a subsequence).

Theorem 1 implies the existence of weak solutions for (20)-(26) with regularity (27)-(30). However, the regularity (27) only implies that the right-hand side of (21) is in $L^{2}(0,T;{\bf L}^{1}(\Omega))\hookrightarrow L^{2}(0,T;{\bf W}^{-1,p}(\Omega))$ for all $1\leq p<2$. Whether the right-hand side of (21) is in $L^{2}(0,T;{\bf H}^{-1}(\Omega))$ is unknown and requires further investigation.

4 Proof of Theorem 1

In this section, we prove Theorem 1 by using a compactness argument. We first introduce some standard notations of finite element spaces in Section 4.1, and then present energy estimates for the numerical solutions in Section 4.2. In Section 4.3 we utilize the compactness of the numerical solutions to prove the existence of a subsequence (in every sequence of numerical solutions) that converges to a weak solution of the PDE problem.

Throughout this paper, we denote by $C$ a generic positive constant which could be different at different places but would be independent of $n$, $h$, and $\tau$. To simplify notation, we use the abbreviations $W^{k,p}=W^{k,p}(\Omega)$, $L^{p}=L^{p}(\Omega)$ and $H^{k}=H^{k}(\Omega)$ for $k\geq 0$ and $1\leq p\leq\infty$.

4.1 Prelimiaries

It is known that the Taylor–Hood finite element space $\mathring{\bf S}_{h}^{r+1}\times S_{h}^{r}/{\mathbb{R}}$ ($r\geq 1$) satisfies the following discrete inf-sup condition for some constant $\gamma>0$:

[TABLE]

Over the finite element spaces $\mathring{S}_{h}^{r+1}$ and $\mathring{\bf S}_{h}^{r+1}$, we define the $L^{2}$ projection $P_{h}$ in (45)-(46). Besides, we also define the $L^{2}$ projection $\widetilde{P}_{h}:L^{2}(\Omega)\rightarrow S_{h}^{r}$ (without enforcing boundary conditions), i.e.,

[TABLE]

We denote by $W^{-1,p}(\Omega)$ the dual space of $W^{-1,p^{\prime}}_{0}(\Omega)$. Under the assumptions on the triangulation in Section 3, the $L^{2}$ projections defined above satisfy the following standard estimates for $1\leq p\leq\infty$:

[TABLE]

The $L^{p}$ and $W^{1,p}$ estimates above were proved in [10, Section 2] for meshes more general than quasi-uniformity assumed in this paper. Since $W^{-1,p}(\Omega)$ is the dual space of $W^{-1,p^{\prime}}_{0}(\Omega)$, the self-adjointness of the $L^{2}$-projection and (49) together imply (51), which is used in the proof for the convergence of numerical solutions; see the last inequality of (4.2).

For any ${\bm{v}}\in{\bf H}^{1}_{0,{\rm div}}(\Omega)$ we define $Q_{h}{\bm{v}}$ to be the Fortin projection of ${\bm{v}}$ onto $S_{h}^{r+1}$, satisfying

[TABLE]

which has the following property (cf. [17] and [20, Lemma 3.4])

[TABLE]

Furthermore, over the finite element spaces $\mathring{S}_{h}^{r+1}$, the following inverse inequality holds; see [15, 44].

[TABLE]

for all $\chi_{h}\in\mathring{S}_{h}^{r+1}$, and $1\leq\tilde{q}\leq q\leq\infty$, $m=0,1$.

4.2 Energy estimate

For any sequence of functions $\omega_{h}^{n}$, $n=0,1\dots,N$, we let $\omega_{h,\tau}$ denote a piecewise linear function in time defined by

[TABLE]

for $t\in(t_{n-1},t_{n}]$ and $n=1,2,\dots,N$. Then we have the following energy estimate for the numerical solutions.

Proposition 2**.**

The numerical scheme (38)-(42) admits a unique solution $(A_{h}^{n},{\bm{u}}_{h}^{n},p_{h}^{n})$, which satisfies the following estimates for $s\in(1,2)$:

[TABLE]

Furthermore, the finite element solution $(A_{h,\tau}^{\pm},{\bm{u}}_{h,\tau}^{\pm})$ satisfies the estimate

[TABLE]

To prove Proposition 2 we need to use the following two lemmas.

Lemma 3**.**

There exists $\alpha>0$ (depending only on the domain $\Omega$) such that

[TABLE]

for any $\chi_{h}\in\mathring{S}_{h}^{r+1}$. Furthermore, if $h_{n}\rightarrow 0$ as $n\rightarrow\infty$, and $\Delta_{h}\chi_{h_{n}}$ is bounded in $L^{2}(\Omega)$, then $\chi_{h_{n}}$ is compact in $W^{1,4+\alpha}(\Omega)$.

*Proof.   * For any given $\chi_{h_{n}}\in\mathring{S}_{h_{n}}^{r+1}$, let $\chi_{n}\in H^{1}_{0}(\Omega)$ be the weak solution of the PDE problem

[TABLE]

with homogeneous Dirichlet boundary condition. Thus $\chi_{n}$ satisfies the weak formulation

[TABLE]

for any $\chi\in\mathring{S}_{h_{n}}^{r+1}$, which implies $\chi_{h_{n}}$ is the Ritz projection of $\chi_{n}$.

On the one hand, as the solution of the PDE problem (65), $\chi_{n}$ satisfies the standard PDE estimate

[TABLE]

for some constant $\delta\in(0,\frac{1}{2}]$ depending on the maximal interior angle of the domain $\Omega$. The estimate (66) is a consequence of [12, Corollary 3.9, with fractional $k$]; also see [13, p. 30] and [11, (23.3)]. Since $H^{\frac{3}{2}+\delta}(\Omega)$ is compactly embedded into $W^{1,4+\alpha}(\Omega)$ for $\alpha\in(0,\frac{8\delta}{1-2\delta})$ in a polygon (cf. [1, Theorem 7.34], with $p=q=2$ and $s=\frac{1}{2}+\delta$), it follows that

[TABLE]

and the set of functions $\{\chi_{n}:n=1,2,\dots\}$ is compact in $W^{1,4+\alpha}(\Omega)$.

On the other hand, as the Ritz projection of $\chi_{n}$, the finite element function $\chi_{h_{n}}$ satisfies the standard error estimate

[TABLE]

as $n\rightarrow\infty$. By using the triangle inequality and the inverse inequality of the finite element space, we have

[TABLE]

where we have used the fact $\frac{1}{2}+\delta-\frac{2+\alpha}{4+\alpha}>0$ for $\alpha\in(0,\frac{8\delta}{1-2\delta})$ in the last inequality. Thus, by (67), we obtain

[TABLE]

Since $\{\chi_{n}:n=1,2,\dots\}$ is compact in $W^{1,4+\alpha}(\Omega)$ and with the above result, the set of functions $\{\chi_{h_{n}}:\,n=1,2,\dots\}$ is also compact in $W^{1,4+\alpha}(\Omega)$. This completes the proof of Lemma 3.

Lemma 4**.**

There exists $\alpha>0$ (depending only on the domain $\Omega$) such that the function $\varphi_{j,h}$ deteremind by (38) converges to $\varphi_{j}$ in the following sense:

[TABLE]

Furthermore,

[TABLE]

*Proof.  * Let $\chi$ be a smooth cut-off function such that $\chi=1$ in a neighborhood of $\Omega_{j}$ and $\chi=0$ on $\cup_{i\neq j}\partial\Omega_{i}$. Then $\varphi_{j}-\chi$ is the solution of

[TABLE]

which implies (similar as (66))

[TABLE]

Therefore, $\|\varphi_{j}\|_{H^{\frac{3}{2}+\delta}}\leq C$. In view of the embedding $H^{\frac{3}{2}+\delta}(\Omega)\hookrightarrow W^{1,4+\alpha}(\Omega)$ for $\alpha\in(0,\frac{8\delta}{1-2\delta})$ (we denote by “$\hookrightarrow$” continuously embedding), we have $\nabla\times\varphi_{j}\in L^{4+\alpha}(\Omega)$. Since $\varphi_{j,h}$ is the finite element solution of $\varphi_{j}$, it follows that

[TABLE]

where $P_{h}$ is the $L^{2}$-projection operator, satisfying (49). By using the inverse inequality we obtain

[TABLE]

For $\alpha\in(0,\frac{8\delta}{1-2\delta})$ we have $\frac{1}{2}+\delta-\frac{2+\alpha}{4+\alpha}>0$. Since $H^{\frac{3}{2}+\delta}(\Omega)\hookrightarrow W^{\frac{3}{2}+\delta-\frac{2+\alpha}{4+\alpha},4+\alpha}(\Omega)$, it follows that

[TABLE]

Threfore we have

[TABLE]

which implies (68).

Let $M=(M_{ij})$ and $M_{h}=(M_{ij,h})$ be two $m\times m$ matrices, with $M_{ij}=(\nabla\times\varphi_{j},\nabla\times\varphi_{i})$ and $M_{ij,h}=(\nabla\times\varphi_{j,h},\nabla\times\varphi_{i,h})$. Since $M=(M_{ij})$ is positive definite and (68) implies $\lim_{h\rightarrow 0}M_{h}=M$, it follows that $\lim_{h\rightarrow 0}M_{h}^{-1}=M^{-1}$. Then (39) and (23) imply

[TABLE]

With above results, we start to prove Proposition 2.

*Proof of Proposition 2.   * First, we start with estimating $\|\nabla A_{h}^{n}\|_{L^{2}}$ and $\|{\bm{u}}_{h}^{n}\|_{L^{2}}$. To this end, we substitute $a_{h}=-\Delta_{h}A_{h}^{n}$, ${\bm{v}}_{h}={\bm{u}}_{h}^{n}$ and $q_{h}=p_{h}^{n}$ into (40)-(41). Then we obtain

[TABLE]

Since

[TABLE]

summing up (75)-(76) and dropping the square terms in the above two equations yield

[TABLE]

By choosing a sufficiently small $\epsilon$, the inequality above implies

[TABLE]

By summing up the inequality above for $n=1,\dots,m$, with $1\leq m\leq N$, we obtain

[TABLE]

where we have used the boundedness of $\|\nabla A_{h}^{0}\|_{L^{2}}$ and $\|{\bm{u}}_{h}^{0}\|_{L^{2}}$, which are direct consequences of the definitions of ${\bm{u}}_{h}^{0}=P_{h}{\bm{u}}_{0}$ and $A_{h}^{0}$ in (44). The estimate (4.2) also implies that, by setting $A_{h}^{n-1}=J^{n}=0$ and ${\bm{u}}_{h}^{n-1}={\bm{f}}^{n}={\bm{0}}$, the homogeneous linear system given by (38)-(42) has only zero solution $A_{h}^{n}=0$ and ${\bm{u}}_{h}^{n}=0$.

In this case, (41) reduces to

[TABLE]

By using the inf-sup condition (48), we further derive $p_{h}^{n}=0$. This implies the existence and uniqueness of solutions of the linear system (38)-(42).

Second, we estimate \big{\|}\frac{A_{h}^{n}-A_{h}^{n-1}}{\tau}\big{\|}_{L^{2}}. To this end, we use (40) and note that for any $q,\bar{q}>2$ satisfying $\frac{1}{\bar{q}}+\frac{1}{q}=\frac{1}{2}$

[TABLE]

The inequality above implies (via the duality argument)

[TABLE]

which holds for all $q\in(2,4)$ with $\bar{q}=\frac{2q}{q-2}$. By using Lemma 3 and (79), we have

[TABLE]

Therefore,

[TABLE]

where we have used Lemma 4 in the last inequality. For any given $1<s<2$ one can choose $q=\frac{4s}{3s-2}\in(2,4)$ so that $(2-\frac{4}{q})\frac{2s}{2-s}=2$ and

[TABLE]

where we have used (79) again. Substituting the estimate above into (80) yields

[TABLE]

Third, we present the estimate for \big{\|}\frac{{\bm{u}}_{h}^{n}-{\bm{u}}_{h}^{n-1}}{\tau}\big{\|}_{{\bf H}^{1}_{0,{\rm div}}(\Omega)^{\prime}}. From (41) we can derive that

[TABLE]

where

[TABLE]

Substituting ${\bm{v}}_{h}=\frac{{\bm{u}}_{h}^{n}-{\bm{u}}_{h}^{n-1}}{\tau}$ into (82) yields

[TABLE]

where we have used (42). Then,

[TABLE]

From (56), we have $(\nabla\cdot{\bm{v}},p_{h}^{n})=(\nabla\cdot Q_{h}{\bm{v}},p_{h}^{n})=0$ for any ${\bm{v}}\in{\bf H}^{1}_{0,{\rm div}}(\Omega)$ and thus

[TABLE]

By using Hölder’s inequality and (56), we have

[TABLE]

where we have used the inverse inequality in the second to last inequality, and (51) in the last inequality. The two estimates above imply

[TABLE]

From (82) we obtain for any $q>1$

[TABLE]

where we have used Sobolev embeddding $L^{q}(\Omega)\hookrightarrow H^{-1}(\Omega)$ ($\forall\,q>1$) in the second to last inequality, and Hölder’s inequality with $\frac{1}{q}=\frac{1}{2}+\frac{1}{2q/(2-q)}$ in the last inequality. In particular, for $1<q<\frac{4}{3}$ we have $\frac{2q}{2-q}<4$ and so

[TABLE]

Substituting the three inequalities above into (4.2) yields

[TABLE]

for all $q\in(1,4/3)$, where we have used (79) in the second to last inequality and Young’s inequality at the last step. For any $s\in(1,2)$ we can choose $q$ to be sufficiently close to $1$ so that $s(3-\frac{2}{q})\leq 2$ and $s(5-\frac{4}{q})\leq 2$, and therefore

[TABLE]

Substituting (4.2) into (86) yields

[TABLE]

Thus, (79), (81) and (90) imply the estimates (2) and (63). The proof of Proposition 2 is complete.

It is easy to see that the numerical scheme (38)-(42) implies that the following equations hold

[TABLE]

for all $a_{h}\in L^{2}(0,T;S_{h}^{r+1})$, ${\bm{v}}_{h}\in L^{2}(0,T;{\bf S}_{h}^{r+1})$ and $q_{h}\in L^{2}(0,T;S_{h}^{r})$.

In the next subsection we pass to the limit in (91)-(93) by using the energy estimates given in Proposition 2.

4.3 Compactness and convergence

Let $s\in(1,2)$ be fixed. The estimate (2) implies the existence of functions

[TABLE]

such that the following convergence results hold: for any sequence $(h_{n},\tau_{n})\rightarrow 0$, there exists a subsequence $(h_{n_{m}},\tau_{n_{m}})$, $m=1,2,\dots$, such that

[TABLE]

where (95) is an obvious consequence of (94). Furthermore, we have $F=\Delta A\in L^{2}(0,T;L^{2}(\Omega))$. To see this, we test $F$ by a function $v\in L^{2}(0,T;H^{1}_{0}(\Omega))$. Then we have

[TABLE]

where we have used (95). Since

[TABLE]

it follows that the function $F\in L^{2}(0,T;L^{2}(\Omega))$ satisfies

[TABLE]

This implies $\Delta A=F\in L^{2}(0,T;L^{2}(\Omega))$ and therefore

[TABLE]

Moreover, the initial value of the limit function $A$ must be equal to $A_{0}$. This can be proved in the following way: on the one hand, for any given smooth function $a$ such that $a|_{t=T}=0$,

[TABLE]

on the other hand,

[TABLE]

where we have used the convergence of $A_{h_{n_{m}},\tau_{n_{m}}}$ in (94) and the convergence of $A_{h_{n_{m}}}^{0}$ to $A_{0}$ in $H^{1}_{0}(\Omega)$. The latter is a simple consequence of the convergence theory for the elliptic problem (26) and (44). Therefore,

[TABLE]

which implies $A|_{t=0}=A_{0}$ in $L^{2}(\Omega)$. Similarly, we also have ${\bm{u}}|_{t=0}={\bm{u}}_{0}$ in ${\bf H}^{1}_{0,{\rm div}}(\Omega)^{\prime}$. This proves

[TABLE]

In Proposition 2 we have proved

[TABLE]

Then Lemma 3 and the Aubin–Lions–Simon lemma (cf. [6, Theorem II.5.16]) imply that $A_{h_{n_{m}},\tau_{n_{m}}}$ is compact in $L^{2}(0,T;W^{1,4+\alpha}(\Omega))$, and thus one can choose the subsequence to have the following property:

[TABLE]

Since $A_{h_{n_{m}},\tau_{n_{m}}}-A$ is bounded in $L^{\infty}(0,T;H^{1}(\Omega))$ and convergent to zero in $L^{2}(0,T;W^{1,4+\alpha}(\Omega))$, and (interpolation inequality)

[TABLE]

for $\frac{1-\theta}{2}+\frac{\theta}{4+\alpha}=\frac{1}{4}$ and $\frac{1}{q_{\theta}}=\frac{1-\theta}{\infty}+\frac{\theta}{2}$, it follows that

[TABLE]

Since $H^{1}(\Omega)\hookrightarrow\hookrightarrow L^{q}(\Omega)\hookrightarrow{\bf H}^{1}_{0,{\rm div}}(\Omega)^{\prime}$ for all $2\leq q<\infty$ and

[TABLE]

the Aubin–Lions–Simon lemma (cf. [6, Theorem II.5.16]) implies that ${\bm{u}}_{h_{n_{m}},\tau_{n_{m}}}$ is compact in $L^{2}(0,T;L^{q}(\Omega))$ for all $2\leq q<\infty$, and thus one can choose the subsequence to have the following property:

[TABLE]

The estimates above are for the piecewise linear functions $A_{h_{n_{m}},\tau_{n_{m}}}$ and ${\bm{u}}_{h_{n_{m}},\tau_{n_{m}}}$. For the piecewise constant functions $A_{h_{n_{m}},\tau_{n_{m}}}^{\pm}$ and ${\bm{u}}_{h_{n_{m}},\tau_{n_{m}}}^{\pm}$ we have similar estimates:

[TABLE]

Since ${\bm{u}}_{h_{n_{m}},\tau_{n_{m}}}-{\bm{u}}$ is bounded in $L^{\infty}(0,T;L^{2}(\Omega))$ and convergent to zero in $L^{2}(0,T;L^{4}(\Omega))$, and (interpolation inequality)

[TABLE]

it follows that

[TABLE]

Then Lemma 4 and (110)-(111) imply

[TABLE]

(110) and (113) imply

[TABLE]

(109) and (112) imply

[TABLE]

Lemma 4, (108) and (111) imply

[TABLE]

For any given $a\in L^{\infty}(0,T;L^{2}(\Omega))$, its $L^{2}$ projection $a_{h}=P_{h}a$ converges to $a$ strongly in $L^{p^{\prime}}(0,T;L^{2}(\Omega))$ for the number $p>1$ in (4.3).

Similarly, for any given ${\bm{v}}\in L^{\infty}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))\hookrightarrow L^{\frac{4}{3}}(0,T;L^{\frac{8}{7}}(\Omega)^{2})^{\prime}\cap L^{2}(0,T;L^{\frac{4}{3}}(\Omega)^{2})^{\prime}$, its Fortin projection ${\bm{v}}_{h}=Q_{h}{\bm{v}}$ converges to ${\bm{v}}$ strongly in $L^{s^{\prime}}(0,T;{\bf H}^{1}_{0}(\Omega))$ and satisfies $(q_{h},\nabla\cdot{\bm{v}}_{h})=0$ for all $q_{h}\in S_{h}^{r}$.

Then, by taking the limit of a subsequence in (91)-(92), we obtain

[TABLE]

for all $a\in L^{\infty}(0,T;L^{2}(\Omega))$ and ${\bm{v}}\in L^{\infty}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$.

The convergence results (94)-(100) show that $A$ and ${\bm{u}}$ possess the regularity in (27)-(28), except the regularity ${\bm{u}}\in L^{2}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$.

Finally, we show that ${\bm{u}}\in L^{2}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$. In fact, for any given $q\in L^{2}(0,T;L^{2}(\Omega))$, its $L^{2}$ projection $q_{h}=\widetilde{P}_{h}q$ converges to $q$ strongly in $L^{2}(0,T;L^{2}(\Omega))$. By taking limit in (93) we obtain

[TABLE]

which implies $\nabla\cdot{\bm{u}}=0$. Since ${\bm{u}}\in L^{2}(0,T;{\bf H}^{1}_{0}(\Omega))$ and $\nabla\cdot{\bm{u}}=0$, it follows that ${\bm{u}}\in L^{2}(0,T;{\bf H}^{1}_{0,{\rm div}}(\Omega))$.

Therefore, (102) and (118)-(119) imply that $(A,{\bm{u}},(\varphi_{j})_{j=1}^{m},(\beta_{j})_{j=1}^{m})$ is a weak solution of (20)-(26). Furthermore, (103), (105) and Lemma 4 imply the convergence results in Theorem 1.

5 Numerical examples

In this section we present two numerical examples to illustrate the convergence of the proposed numerical method in nonconvex and nonsmooth domains. All the numerical tests are done by using FreeFEM++ with double precision.

Example 5.1**.**

In the first example, we consider the MHD equations

[TABLE]

in a simply connected L-shape domain $\Omega$ whose longest side has unit length, centered at the origin; see Figure 1 (left). The source terms ${\bm{f}}$, $g$ and

[TABLE]

are calculated by substituting the following exact solution into (120)-(122):

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$ and $\theta={\rm arg}(x+iy)$, with $(x,y)\in\Omega$, and $\Phi(r)$ is a $C^{3}(\Omega)$ cut-off function. Here,

[TABLE]

and $\Upsilon(r)$ is the unique $7^{\rm th}$ order polynomial satisfying the conditions $\Upsilon^{\prime}(0.1)=\Upsilon^{\prime\prime}(0.1)=\Upsilon^{\prime\prime\prime}(0.1)=\Upsilon(0.4)=\Upsilon^{\prime}(0.4)=\Upsilon^{\prime\prime}(0.4)=\Upsilon^{\prime\prime\prime}(0.4)=0$ and $\Upsilon(0.1)=0.1$. The constructed solutions in (123)-(124) satisfy the boundary conditions

[TABLE]

We solve the equations (120)-(122) up to time $T=1$ by the proposed numerical method (38)-(42) with $r=1$, i.e., with the P2 element for $A$ and P2-P1 elements for $({\bm{u}},p)$. The numerical solution of magnetic field is given by ${\bm{H}}_{h}^{n}=\nabla\times A_{h}^{n}.$ Since the L-shape domain is simply connected, it follows that $m=0$ (the constants $\beta_{j}$, $j=1,\dots,m$, are not needed). The L-shape domain is triangulated quasi-uniformly with $M$ nodes per unit length on each side, and we denote by $h=1/M$ for simplicity.

We compare the numerical solutions with the exact solution given by (123)-(124) and present the errors of the numerical solutions in Table 1. For comparison, we also present the numerical results of a “direct $H^{1}$-conforming FEM” in Table 2, with the same time-stepping scheme as the proposed method (38)-(42). In particular, the direct $H^{1}$-conforming FEM seeks ${\bm{H}}_{h}^{n}\in{\bf S}_{h}^{2}$, ${\bm{u}}_{h}^{n}\in\mathring{\bf S}_{h}^{2}$ and $p_{h}^{n}\in S_{h}^{1}$, with ${\bm{H}}_{h}^{n}\cdot{\bm{n}}=0$ on $\partial\Omega$, such that the following equations hold for all test functions ${\bf a}_{h}\in{\bf S}_{h}^{2}$, ${\bm{v}}_{h}\in\mathring{\bf S}_{h}^{2}$ and $q_{h}\in S_{h}^{1}$ with ${\bf a}_{h}\cdot{\bm{n}}=0$ on $\partial\Omega$:

[TABLE]

The numerical results in Tables 1 and 2 show that the proposed method has slightly higher accuracy than the $H^{1}$-conforming FEM in computing the magnetic field ${\bm{H}}$ (with the same degree of finite elements and similar computational complexity). The convergence of the proposed numerical method is proved in Theorem 1 (also see Remark 3.1), while the convergence of the direct $H^{1}$-conforming FEM remains open in nonconvex and nonsmooth domains.

Example 5.2**.**

In the second example, we consider the MHD equations (120)-(122) in a multi-connected domain shown in Figure 1 (right), with the source terms ${\bm{f}}$, $g$ and

[TABLE]

determined by the exact solution

[TABLE]

where $\varphi$ is the harmonic function satisfying

[TABLE]

In this case $m=1$ and $\beta_{1}=1$.

We solve the MHD equations up to time $T=1$ by the proposed numerical method (38)-(42) with the P2 element for $A$ and P2-P1 elements for $({\bm{u}},p)$, and compare the numerical solutions with the exact solution given by (129)-(130) (where $\varphi$ is approximated numerically). The domain is triangulated quasi-uniformly, with $M$ nodes per unit length on each side, and we denote by $h=1/M$ for simplicity. The errors of the numerical solutions are presented in Table 3. For comparison, we also present the numerical results of the direct $H^{1}$-conforming FEM (5.1)-(128) with the P2 element for ${\bm{H}}$ and P2-P1 elements for $({\bm{u}},p)$ in Table 4. Numerical results in Tables 3–4 show that the proposed method is much more accurate than the direct $H^{1}$-conforming FEM in such a multi-connected nonsmooth domain. The reason may be that the proposed method approximates the harmonic part $\nabla\times\varphi$ at the initial time in a more accurate way than the direct $H^{1}$-conforming FEM, which cannot separate $\nabla\times\varphi$ from ${\bm{H}}$.

6 Conclusion and remarks

We have proposed a fully discrete and linearized $H^{1}$-conforming Lagrange finite element method for solving a magnetic potential formulation of the two-dimensional MHD equations in general domains that can be non-convex, nonsmooth and multi-connected. The theoretical analysis shows that, for given source terms and initial data with regularity (6), a subsequence of numerical solutions converges in $L^{2}(0,T;{\bf L}^{2}(\Omega))$ to a weak solution of the MHD equations without any mesh restriction. The uniqueness of weak solutions for the two-dimensional MHD equations in general domains remains open. If the weak solution is unique, then the numerical solutions converge to the unique weak solution as $(\tau,h)\rightarrow(0,0)$ (without taking a subsequence).

The results proved in this paper can be generalized to more general $H^{1}$-conforming finite element spaces. For example, the finite element space $\mathring{\bf S}_{h}^{r+1}\times S_{h}^{r}/{\mathbb{R}}$ for $({\bm{u}},p)$ can be replaced by other $H^{1}$-conforming finite element spaces $X_{h}\times M_{h}$ satisfying the inf-sup condition (4.1), meanwhile with the following approximation property:

[TABLE]

The finite element space $\mathring{S}_{h}^{r+1}$ for $A$ can be replaced by any $H^{1}$-conforming finite element space $V_{h}\subset H^{1}_{0}(\Omega)$ such that

[TABLE]

The assumption of perfect-conductor-type boundary conditions for the magnetic field is an ideal situation. In realistic problems, the magnetic field is often not confined to the fluid region, but extends throughout space, satisfying interface conditions (rather than boundary conditions) on the surface of the fluid region. In this case, the numerical approximation of solutions to the full incompressible MHD equations in the whole space, satisfying interface conditions on the surface of the fluid region, is still an open problem. The present paper, focused on two-dimensional domains with perfect-conductor-type boundary conditions, could be an incremental step towards the solution.

Funding. The work of Buyang Li was supported in part by the Hong Kong RGC grant 15301818. The work of Jilu Wang was supported in part by the USA National Science Foundation grant DMS-1315259 and by the USA Air Force Office of Scientific Research grant FA9550-15-1-0001. The work of Liwei Xu was supported in part by a Key Project of the Major Research Plan of NSFC grant no. 91630205 and the NSFC grant no. 11771068.