A modified characteristics projection finite element method for unsteady incompressible Magnetohydrodynamics equations

1.   Introduction

Magnetohydrodynamics(MHD) is a subject that studies the motion of conductive fluids in magnetic fields. It is mainly used in astrophysics, controlled thermonuclear reactions and industry. The system of MHD equations is a coupled equation system of Navier-Stokes equations of fluid mechanics and Maxwell's equations of electrodynamics (see^[1,2]). In this paper, we consider the unsteady incompressible MHD equations as follows:

where $u$ is velocity, $\nu$ kinematic viscosity, $\mu$ magnetic permeability, $p$ hydrodynamic pressure, $B$ magnetic field, $f$ body force, and $g$ applied current. The coefficients are the magnetic Reynolds number $Rm$. As in ^[3], we assume that $\Omega$ is a convex domain in $R^d$, $d = 2, 3$ and $t\in [0, T]$. Here, $u_t = \frac{\partial u}{\partial t}$, $B_t = \frac{\partial B}{\partial t}$. The term $u\times B$ explains the effect of hydrodynamic flow on current. The magnetic Reynolds number $Rm$ is moderate or low (from $10^{-2}$ to 1) (see^[4]). For the purposes of simplicity, we choose $Rm = 1$ in our numerical analysis.

System (1.1) is considered in conjunction with the following initial boundary conditions

where ${\bf{n}}$ denotes the outer unit normal of $\partial \Omega$.

Over that last few decades, there has been a lot of works devoted on the numerical solution of MHD flows. In ^[5], some mathematical equations related to the MHD equations have been studied. In ^[6], a weighted regularization approach was applied to incompressible MHD problems. Reference ^[7,8] decoupled the linear MHD problem involving conducting and insulating regions. A mixed finite element method for stationary incompressible MHD equations was given in ^[9]. In order to avoid in-sup condition, Gerbeau ^[10] presented a stabilized finite element method and Salah et al. ^[11] proposed a Galerkin least-square method. In ^[12], a two-level finite element method with a stabilizing sub-grid for the incompressible MHD equations has been shown.

Recently, a convergence analysis of three finite element iterative methods for 2D/3D stationary incompressible MHD equations has been given by Dong and He (see^[13]). In order to continue the in-depth explore of three-dimensional incompressible MHD system, an unconditional convergence of the Euler semi-implicit scheme for the three-dimensional incompressible MHD equations was studied by He ^[14]. In ^[15], a two-level Newton iteration method for 2D/3D steady incompressible MHD equation was proposed by Dong and He. In addition, Yuksel and Ingram have discussed the full discretization of Grank-Nicolson scheme at small magnetic Reynolds numbers (see^[16]). The defect correction finite element method for the MHD equations was shown by Si et. al. ^[17,18,19]. In ^[20], a decoupling penalty finite element method for the stationary incompressible MHD equation was presented. To further study this method, we can refer to ^[21,22,23].

In this paper, a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations will be given. In this method, modified characteristics finite element method and the projection method will be combined for solving the incompressible MHD equations. Both the stability and the optimal error estimates both in $L^2$ and $H^1$ norms for the modified characteristics projection finite element method will be derived. In order to demonstrate the effectiveness of our method, we will present some numerical results at the end of the manuscript. The contents of this paper are divided into sections as follows. To obtain the unconditional stability and optimal error estimates, in section 2, we introduce some notations and present a modified characteristics projection algorithm for the MHD equations. Then, in section 3, we give stability and error analysis and show that this method has optimal convergence order by mathematical induction. In order to demonstrate the effectiveness of our method, some numerical results are presented in the last section.

2.   Notation and preliminaries

In this section, we introduce some notations. We employ the standard Sobolev spaces $H^{k}(\Omega) = W^{k, 2}$, for nonnegative $k$ with norm $\|v\|_{k} = (\sum_{|\beta| = 0}^{k}\|D^{\beta}v\|_{0}^{2})^{1/2}$. For vector-value function, we use the Sobolev spaces $\mbox{H}^{k}(\Omega) = (H^{k}(\Omega))^{d}$ with norm $\|\mbox{v}\|_{k} = (\sum_{i = 1}^{d}\|v_{i}\|_{k}^{2})^{1/2}$(see ^[24,25]). For simplicity, we set:

which are equipped with the norms

In this paper, we will use the following equality

Then, we recall the following Poincaré-Friedrichs inequality in $W$: there holds

with a constant $C > 0$ solely depending on the domain $\Omega$; (see^[26]). Here, we define the following trilinear form $a_1(\cdot, \cdot, \cdot)$ as follows, for all $u\in V$, $w, v\in X$,

and the properties of the trilinear form $a_1(\cdot, \cdot, \cdot)$ ^[27] are helpful in our analysis

where $N > 0$ is a constant, $\gamma_0$ (only dependent on $\Omega$) is a positive constant, and $C_0$ (only dependent on $\Omega$) is an embedding constant of $H^{1}(\Omega)\hookrightarrow L^{4}(\Omega)$ ($\hookrightarrow$ denotes the continuous embedding).

Furthermore, we define the Stokes operator $\mathcal{A}:D(\mathcal{A})\longrightarrow \tilde{H}$ such that $\mathcal{A} = -P\triangle$, where $D(\mathcal{A}) = H^{2}(\Omega)^{d}\cap X_{0}$ and $P$ is an $L^2$-projection from $L^2(\Omega)^d$ to $\tilde{H} = \{v\in L^2(\Omega)^d; \nabla\cdot v = 0, v\cdot n|_{\partial\Omega} = 0\}$.

Next, we recall a discrete version of Gronwall's inequality established in ^[14].

Lemma 2.1 Let $C_{0}, a_n, b_n$ and $d_n$, for integers $n\geq0$, be the non-negative numbers such that

Then

Throughout this paper, we make the following assumptions on the prescribed date for the problem (1.1), which satisfies the regularity of the data needed for our main results.

Assumption (A1)(see^[27,28]): Assume that the boundary of $\Omega$ is smooth so that the unique solution $(v, q)$ of the steady Stokes problem

for prescribed $f\in L^{2}(\Omega)^{3}$ satisfies

and Maxwell's equations

for the prescribed $g\in L^{2}(\Omega)^{3}$ admits a unique solution $B\in W_{0}$ which satisfies

Assumption (A2): The initial data $u_0, B_0, f$ satisfy

Assumption (A3): We assume that the MHD problem admits a unique local strong solution $(u, p, B)$ on $(0, T)$ such that

Let $0 = t_{0} < t_{1} < \cdots < t_{N} = T$ be a uniform partition of the time interval $[0, T]$ with $t_{n} = n\tau$, and $N$ being a positive integer. Set $u^{n} = (\cdot, t_{n})$, For any sequence of functions $\{f^{n}\}_{n = 0}^{N}$, we define

At the initial time step, we start with $U^{0} = \tilde{U}^{0} = u_{0}, B^{0} = \tilde{B}^{0} = b_{0}$ and an arbitrary $P^{0}\in M\cap H^{1}(\Omega)$. Then, we can give

Then, we give the algorithm as follows

Algorithm 2.1 (Time-discrete modified characteristics projection finite element method)

Step 1. Solve $B^{n+1}$, for instance

where

Step 2. Find $\tilde{U}^{n+1}$, such that

Step 3. Calculate $(U^{n+1}, P^{n+1})$, for example

The corresponding weak form is

Step 1. Solve $B^{n+1}$, for instance

Step 2. Find $\tilde{U}^{n+1}$, such that

Step 3. Calculate $(U^{n+1}, P^{n+1})$, for example

Remark 2.1 As we all know $(B^{n}\times \mbox{curl}B^{n+1}, v) = (B^{n}\cdot\nabla v, B^{n+1})-(v\cdot\nabla B^{n}, B^{n+1})$. In order to improve our algorithm, we change it as $(B^{n}\cdot\nabla v, B^{n+1})-(v\cdot\nabla B^{n}, B^{n})$.

We denote $T_{h}$ be a regular and quasi-uniform partition of the domain $\Omega$ into the triangles for $d = 2$ or tetrahedra for $d = 3$ with diameters by a real positive parameter $h(h\rightarrow 0)$. The finite element pair $(X_{h}, M_{h}, W_{h})$ is constructed based on $T_{h}$.

Let $W_{0n}^h = X_h\times W_h$. There exits a constant $\beta_0 > 0$ (only dependent on $\Omega$) such that

There exits a mapping $r_h\in \mathcal{L}(H^2(\Omega)^2\cap V, X_h)$ satisfying

and a $L^2$-projection operator $\rho_h:M\rightarrow M_h$ satisfying

and a mapping $R_h\in \mathcal{L}(H^2(\Omega)^2\cap V_n, W_h)$ satisfying

Here, we define the discrete Stokes operator $A_{h} = -P_h\triangle_h$ and $\triangle_h$ (see ^[15] and the references therein) defined by

Meanwhile, we define the discrete operator $A_{1h} B_h = R_{1h}(\nabla_h\times \nabla\times B_h+\nabla_h\nabla\cdot B_h)\in W_h$ (see ^[15]) as follows

We define the Stokes projection $(R_{h}(u, p), Q_{h}(u, P)):(X, M)\rightarrow (X_{h}, M_{h})$ by

By classical FEM theory (^[25,29,30]), we have the following results.

Lemma 2.2 Assume that $u \in H_{0}^{1}(\Omega)^{d}\cap H^{r+1}(\Omega)^{d}$ and $p \in L_{0}^{2}(\Omega)\cap H^{r}(\Omega)$. Then,

and

Lemma 2.3 If $(u, p)\in W^{2, k}(\Omega)^{d}\times W^{1, k}(\Omega)$ for $k > d$,

Algorithm 2.2 (The fully discrete modified characteristics projection finite element method)

Step 1. Solve $B_{h}^{n+1}$, for instance

where

Step 2. Find $\tilde{U}^{n+1}$, such that

Step 3. Calculate $(u_{h}^{n+1}, p_{h}^{n+1})$, for example

3.   Unconditional stability and convergence

3.1.   Error estimation for the time-discrete method

Here, we make use of the following notation:

Lemma 3.1 (see^[31]). Assume that $g_{1}, g_{2}, \rho$ are three functions defined in $\Omega$ and $\rho|_{\partial\Omega} = 0$. If

then

where $1/q_{1}+1/q_{2} = 1/q, 1 < q < \infty$.

The following theorem gives insight on the error estimate between the time-discrete solution and the solution of the time-dependent MHD system (1.1).

Theorem 3.1 Suppose that assumption A2–A3 are satisfied, there exists a positive $C > 0$ such that

Proof. Firstly, we prove the following inequality

by mathematical induction for $m = 0, 1, ..., N$.

Since $U^{0} = u^{0}, B^{0} = b^{0}$ the above inequality hold for $m = 0$.

We assume that (3.1) and (3.2) hold for $m\leq n$ for some integer $n\geq0$. By Sobolev embedding theory, we get

and

when $\tau\leq\tau_{1}$ for some positive constant $\tau_{1}$.

To prove (3.1) and (3.2) for $m = n+1$, we rewrite (1.1) by

where $u^{n+1} = U(t_{n+1}), b^{n+1} = B(t_{n+1})$, and

and

define the truncation error. We can refer to (^[31]), there holds that

The corresponding weak form is

Combining (2.13) and (2.14), we obtain

Then, we rewrite (3.14) and (2.12) as

and

Subtracting (3.15) and (3.16) from (3.11) and (3.13), respectively, leads to

and

Testing (2.9) by $\tau e^{n+1}$, since $\nabla\cdot e^{n+1} = 0$, we arrive at

Testing (2.9) by $\tau e^{n}$, we get

Subtracting (3.20) from (3.19), we can obtain

Let $v = 2\tau \tilde{e}^{n+1}$ in (3.17), we obtain

Let $\psi = 2\mu\tau \varepsilon^{n+1}$ in (3.18), we get

Taking sum of (3.22) and (3.23) yields

From (2.9), we have

Then, we can obtain

On the other hand, using Lemma 3.1, we deduce that

Substituting these above inequality into (3.24), we obtain

Taking sum of the (3.25) for $n$ from $0$ to $m\leq N$ and using discrete Gronwall's inequality Lemma 2.1, we get

Owing to (3.21), we have

To acquire the $H^{1}$ estimates, we take $v = 2\tau D_{\tau}e^{n+1}$ in (3.17) to get

Then, we take $\psi = 2\mu\tau D_{\tau}\varepsilon^{n+1}$ in (3.18) to get

Combining (3.28) and (3.29), we obtain

Similarly, by Lemma 3.1, we have

Substituting these above inequality into (3.30), we obtain

Taking sum of the (3.31) for $n$ from $0$ to $m\leq N$ and using discrete Gronwall's inequality Lemma 2.1, we arrive at

Furthermore, the standard result for (3.17) and (3.18) with $p = 2$, respectively, leads to

Thanks to (3.19), we can obtain

Then, we have

From (3.35), we can arrive at

From (2.7) and (2.5), we have

and

By (2.9), we have

Considering the identity $\nabla^2\tilde{U}^{n+1} = \nabla\nabla\cdot\tilde{U}^{n+1}-\nabla\times\nabla\times\tilde{U}^{n+1}$ and (3.46), we can obtain

Now, the standard result for the Stokes system (3.17) and (3.18) with $p = d^{*} > 2$, respectively, leads to

By (3.48) and (3.49), we get

Thus, we can obtain

The proof is complete.

3.2.   Error estimation for the fully discrete method

For the sake of simplicity, we denote $R_{h}^{n} = R_{h}(U^{n}, P^{n}), Q_{h}^{n} = Q_{h}(U^{n}, P^{n}), R_{0h} = R_{0h}B^{n}, n = 1, 2, ..., N$, and $e_{h}^{n} = u_{h}^{n}-R_{h}^{n}, \tilde{e}_{h}^{n} = \tilde{u}_{h}^{n}-R_{h}^{n}, \varepsilon_{h}^{n} = B_{h}^{n}-R_{0h}^{n}, n = 1, 2, ..., N$, where $R_{0h}:L^{2}(\Omega)^{3} \rightarrow W_{h}$ is the $L^{2}$-orthogonal projection, where

Theorem 3.2 Suppose that assumption A2–A3 are satisfied, there exists a positive $C > 0$ such that

Proof. We prove following estimates

by mathematical induction for $m = 0, 1, ..., N$.

When $m = 0, u_{h}^{0} = u_{0}, B_{h}^{0} = B_{0}$. It is obvious (3.56) holds at the initial time step. We assume that (3.56) holds for $0\leq m\leq n$ for some integer $n\geq0$. By (2.20), (2.21), (3.54), (3.55), Theorem 3.1 and inverse inequality, we infer

and

Combining (2.23) and (2.24), we obtain

When $m = n+1$, letting $v_{h} = 2\tau \tilde{u}_{h}^{n+1}$ in (3.61), we can obtain

Letting $v_{h} = \tau u_{h}^{n+1}$ in (2.24), since $\nabla\cdot u_{h}^{n+1} = 0$, it follows that

Taking $v_{h} = \tau u_{h}^{n}$ in (2.24), we get

Subtracting (3.64) from (3.63), we can obtain

Using $2(a-b, a) = \|a\|_{0}^{2}-\|b\|_{0}^{2}+\|a-b\|_{0}^{2}$, leads to

Taking $v_{h} = \nabla p_{h}^{n+1}$ in (2.24), we deduce

When $m = n+1$, letting $\psi_{h} = 2\mu\tau B_{h}^{n+1}$ in (2.22), we can obtain

Using $2(a-b, a) = \|a\|_{0}^{2}-\|b\|_{0}^{2}+\|a-b\|_{0}^{2}$, we have

Taking sum of (3.66) and (3.69) yields

Then, we can obtain

Substituting these above estimates into (3.70), we obtain

Taking sum of the (3.71) for $n$ from $0$ to $m\leq N$ and using discrete Gronwall's inequality Lemma 2.1, we get

Letting $v_{h} = 2\tau\mathcal{A}u_{h}^{n+1}$ in (3.61), we obtain

From (3.63), we arrive at

Using $2(a-b, a) = \|a\|_{0}^{2}-\|b\|_{0}^{2}+\|a-b\|_{0}^{2}$, leads to

Letting $\psi_{h} = 2\mu\tau \mbox{curlcurl}B_{h}^{n+1}$ in (2.22), we can obtain

Taking sum of (3.73) and (3.74) yields

Then, we get

Substituting these above estimates into (3.75), we obtain

Taking sum of the (3.76) for $n$ from $0$ to $m\leq N$ and using discrete Gronwall's inequality Lemma 2.1, we get

Then, we rewrite (3.61) and (2.22) as

and

Subtracting (3.15) and (3.16) from (3.77) and (3.78), respectively, leads to

and

Taking $v_{h} = 2\tau e_{h}^{n+1}$ in (2.24), we can obtain

Let $\theta^{n} = U^{n}-R_{h}^{n}, \eta^{n} = B^{n}-R_{0h}^{n}$. By taking $v_{h} = 2\tau e_{h}^{n+1}$ in (3.79), we deduce

By taking $\psi_{h} = 2\mu\tau \varepsilon_{h}^{n+1}$ in (3.80), we have

Combining (3.81) and (3.82), we obtain

Due to (2.19) and Theorem 3.1, we can get $\|\theta^{n}\|_{0}\leq Ch^{2}(\|U^{n}\|_{2}+\|P^{n}\|_{1})\leq Ch^{2}$, $\|\eta^{n}\|_{0}\leq Ch^{2}\|B^{n}\|_{2}$. By Lemma 3.1, there holds that

Substituting these above estimates into (3.83), we obtain

Taking the sum of the (3.76) for $n$ from $0$ to $m\leq N$ and using the discrete Gronwall's inequality Lemma 2.1, we get

Then, we can deduce

Thus, all the results in Theorem 3.2 have been proved.

4.   Numerical results

In order to show the effect of our method, we present some numerical results. Here, we consider the

The boundary conditions and the forcing terms are given by the exact solution. The finite element spaces are chosen as P1b-P1-P1b finite element spaces. Here, we choose $\tau = h^2$ and $h = 1/n$, $n = 8, 16, 24, 32, 40, 48$. Different Reynolds number $Re$ and magnetic Reynolds number $Rm$ are chosen to show the effect of our method, the numerical results were shown in Tables 1–6. Here, we use the software package FreeFEM++ ^[32] for our program.

Tables 1 and 2 give the numerical results for $Re = 1.0$ and $Rm = 1.0$. In order to show the effect of our method for high Reynolds number, we give the numerical results for $Re = 100.0, \ 1000.0$ and $Rm = 100.0, \ 1000.0$. Tables 3 and 4 give the numerical results for $Re = 100.0$ and $Rm = 100.0$. Tables 5 and 6 give the numerical results for $Re = 1000.0$ and $Rm = 1000.0$. It shows that our method is effect for high Reynolds numbers. It shows that the errors are goes small as the space step goes small and the convergence orders are optimal. We can see that $|\int_{\Omega} \nabla\cdot B dx|$ and $|\int_{\Omega} \nabla\cdot {{\bf u}} dx|$ are small, which means that our method can conserve the Gauss's law very well.

5.   Conclusions

In this paper, we have given a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations. In this method, the modified characteristics finite element method and the projection method were combined for solving the incompressible MHD equations. Both the stability and the optimal error estimates in $L^2$ and $H^1$ norms for the modified characteristics projection finite element method have been given. In order to demonstrate the effectiveness of our method, some numerical results were given at the end of the manuscript.

Acknowledgments

The authors would like to thank the anonymous referees for their valuable suggestions and comments, which helped to improve the quality of the paper. This work is supported in part by National Natural Science Foundation of China (No. 11971152), China Postdoctoral Science Foundation (No. 2018M630907) and the Key scientific research projects Henan Colleges and Universities (No. 19B110007).

Conflict of interest

The authors declare there is no conflict of interest in this paper.