Numerical analysis of an efficient second order time filtered backward Euler method for MHD equations

TL;DR

This paper introduces a second-order time filtering method based on backward Euler for MHD equations, enhancing accuracy and stability with minimal additional complexity, and verified through numerical experiments.

Contribution

The paper presents a novel modular time filter for backward Euler that improves accuracy and stability in MHD simulations without significant computational overhead.

Findings

The method achieves second-order convergence.

It maintains long-term stability and conservation properties.

Numerical experiments confirm theoretical results and practical usefulness.

Abstract

The present work is devoted to introduce the backward Euler based modular time filter method for MHD flow. The proposed method improves the accuracy of the solution without a significant change in the complexity of the system. Since time filters for fluid variables are added as separate post processing steps, the method can be easily incorporated into an existing backward Euler code. We investigate the conservation and long time stability properties of the improved scheme. Stability and second order convergence of the method are also proven. The influences of introduced time filter method on several numerical experiments are given, which both verify the theoretical findings and illustrate its usefulness on practical problems.

Tables

Table 1. Table 1: Errors and rates of convergence for the velocity and the magnetic field.

$h=\Delta t$	${\Vert u-u_h\Vert }_{2,1}$	rate	${\Vert B-B_h\Vert }_{2,1}$	rate
1/2	0.30650	-	0.01460	-
1/4	0.08936	1.73	0.00995	0.63
1/8	0.02239	2.01	0.00282	1.81
1/16	0.00559	2.01	0.00071	1.98
1/32	0.00139	2.00	0.00017	2.00
1/64	0.00034	2.04	4.44e-5	2.01

Equations

$$u_t-R{e}^{-1}\Delta u+u\cdot \nabla u-sB\cdot \nabla B+\nabla P$$

$$\nabla \cdot u$$

$$B_t-R{e}_m^{-1}\Delta B+u\cdot \nabla B-B\cdot \nabla u-\nabla \lambda$$

$$\nabla \cdot B$$

$$X$$

$$Q$$

$$V:=\{v\in X,(\nabla \cdot v,q)=0,\forall q\in Q\}.$$

$${\Vert w\Vert }_{L^p(0,T;X)}$$

$${\Vert w\Vert }_{L^{\infty }(0,T;X)}$$

$$\Vert v\Vert \le C_p\Vert \nabla v\Vert ,$$

$$\displaystyle\begin{array}[]{rcl}(u\cdot\nabla v,w)&\leq&C\left\|\nabla u\right\|\left\|\nabla v\right\|\left\|\nabla w\right\|,\\ (u\cdot\nabla v,w)&\leq&C\left\|u\right\|^{1/2}\left\|\nabla u\right\|^{1/2}\left\|\nabla v\right\|\left\|\nabla w\right\|,\\ (u\cdot\nabla v,v)&=&0\end{array}$$

$$\nabla \cdot u_h^n$$

$$\nabla \cdot B_h^n$$

$$\underset{v_h\in X_h}{\inf }\left\{\Vert u-v_h\Vert +h\Vert \nabla (u-v_h)\Vert \right\}$$

$$\underset{B_h\in X_h}{\inf }\left\{\Vert B-B_h\Vert +h\Vert \nabla (B-B_h)\Vert \right\}$$

$$V_h=\{v_h\in X_h:(q_h,\nabla \cdot v_h)=0,\forall q_h\in Q_h\},$$

$$L_{\infty }({\mathbb{R}}_{+},{\mathbf{V}}_h^{\ast })=\{g\in {\Omega }^d\times {\mathbb{R}}_{+}\to {\mathbb{R}}^d,a.e.t>0,\exists M<\infty ,\Vert g(t){\Vert }_{{\mathbf{V}}_h^{\ast }}<M\},$$

$$\Vert w{\Vert }_{{\mathbf{V}}_h^{\ast }}=\underset{v_h\in {\mathbf{V}}_h^{\ast }}{\sup }\frac{(w,v_h)}{\Vert \nabla v_h\Vert }.$$

$${\alpha }_m+\Delta t\sum _{n=0}^m{\beta }_n\le \Delta t\sum _{n=0}^m{\delta }_n{\alpha }_n+\Delta t\sum _{n=0}^m{\xi }_n+M\text{mbox}form\ge 0.$$

$$\displaystyle\alpha_{m}+\Delta t\sum_{n=0}^{m}\beta_{n}\leq\exp\Bigg{(}\Delta t\sum_{n=0}^{m}\beta_{n}\frac{\delta_{n}}{1-\Delta t\delta_{n}}\Bigg{)}\Bigg{(}\Delta t\sum_{n=0}^{m}\xi_{n}+M\Bigg{)}\quad\mbox{for}\quad m\geq 0.$$

$$\frac{1}{\Delta t}({\widetilde{u_h}}^{n+1}-u_h^n,v_h)+R{e}^{-1}(\nabla {\tilde{u}}_h^{n+1},\nabla v_h)+({\tilde{u}}_h^{n+1}\cdot \nabla {\tilde{u}}_h^{n+1},v_h)$$

$$-s({\tilde{B}}_h^{n+1}\cdot \nabla {\tilde{B}}_h^{n+1},v_h)-({\tilde{P}}_h^{n+1},\nabla \cdot v_h)=(f(t^{n+1}),v_h),$$

$$(\nabla \cdot {\tilde{u}}_h^{n+1},q_h)=0,$$

$$\frac{1}{\Delta t}({\tilde{B}}_h^{n+1}-B_h^n,{\chi }_h)+R{e}_m^{-1}(\nabla {\tilde{B}}_h^{n+1},\nabla {\chi }_h)-({\tilde{B}}_h^{n+1}\cdot \nabla {\tilde{u}}_h^{n+1},{\chi }_h)$$

$$+({\tilde{u}}_h^{n+1}\cdot \nabla {\tilde{B}}_h^{n+1},{\chi }_h)+({\tilde{\lambda }}_h^{n+1},\nabla \cdot {\chi }_h)=(\nabla \times g(t^{n+1}),{\chi }_h),$$

$$(\nabla \cdot {\widetilde{B_h}}^{n+1},r_h)=0,$$

$$u_h^{n+1}$$

$$B_h^{n+1}$$

$$P_h^{n+1}$$

$${\lambda }_h^{n+1}$$

$$\mathcal{F}[w_h^{n+1}]=\frac{3}{2}{w}_h^{n+1}-w_h^n+\frac{1}{2}{w}_h^{n-1}$$

$$\frac{1}{\Delta t}(\frac{3}{2}{u}_h^{n+1}-2u_h^n+\frac{1}{2}{u}_h^{n-1},v_h)+R{e}^{-1}(\nabla (\mathcal{F}[u_h^{n+1}]),\nabla v_h)$$

$$(\nabla \cdot (\mathcal{F}[u_h^{n+1}]),q_h)=0,$$

$$\frac{1}{\Delta t}(\frac{3}{2}{B}_h^{n+1}-2B_h^n+\frac{1}{2}{B}_h^{n-1},{\chi }_h)+R{e}_m^{-1}(\nabla \mathcal{F}[B_h^{n+1}],\nabla {\chi }_h)$$

$$(\nabla \cdot (\mathcal{F}[B_h^{n+1}]),r_h)=0,$$

$$G=\left(\begin{array}{cc}\frac{3}{2}& -\frac{3}{4}\\ -\frac{3}{4}& \frac{1}{2}\end{array}\right),$$

$$\displaystyle\left\|\begin{bmatrix}{u}\\ {v}\end{bmatrix}\right\|_{G}^{2}=\big{(}\begin{bmatrix}{u}\\ {v}\end{bmatrix},G\begin{bmatrix}{u}\\ {v}\end{bmatrix}\big{)},$$

$${\Vert u\Vert }_F=(u,Fu).$$

$$C_1{\Vert \left[\begin{array}{c}u\\ v\end{array}\right]\Vert }_G^2\le {\Vert \left[\begin{array}{c}u\\ v\end{array}\right]\Vert }^2\le C_2{\Vert \left[\begin{array}{c}u\\ v\end{array}\right]\Vert }_G^2.$$

Full text

Keywords: time filter, backward Euler, MHD equations

1 Introduction

This paper considers a modular time filter method combined with the backward Euler method for the magnetohydrodynamics (MHD) flow problems. A simple method of incorporating this time filter into an existing code is to add extra lines for each fluid variables, thus it can be considered as a post processing step. As discussed in [14], for ODEs adding such time filter to backward Euler not only increases accuracy from first order to second order, but also reduces spurious oscillations of numerical solutions, preserves A-stability of the method and yields a useful error estimator. Recently, the time filter of [14] was considered for Navier-Stokes equations for constant and variable time steps by DeCaria, Layton and Zhao in [7], resulting a stable, second order time accurate adaptive method with a low complexity.

The goal of this paper is to extend this novel idea from [7] of time accurate flow approximation to the MHD system for constant time steps, which describes the mutual interaction between the magnetic field and electrically conductive fluids. These flows have diverse applications in, e.g., hydrology, geophysics, astrophysics and cooling system designs [6, 18, 27, 28]. It was first presented by Ladyzhenskaya and has been developed in [5, 11, 12, 13, 25, 26]. Using Navier Stokes equations (NSE) and Maxwell equations, the governing equations of MHD system are given by

[TABLE]

in a bounded polyhedral domain $\Omega\subset\mathbb{R}^{d},d\in\{2,3\}$. Here, $u$, $P:=p+\frac{s}{2}|B|^{2}$, $p$ and $B$ denote the unknown velocity, modified pressure, pressure and magnetic field, respectively. The body forces ${f}$ and $\nabla\times g$ are forcing on the velocity and magnetic field, respectively. Also, $Re$ is the Reynolds number, $Re_{m}$ is the magnetic Reynolds number, and $s$ is the coupling number. The Lagrange multiplier (dummy variable) $\lambda$ corresponds to the solenoidal constraint on the magnetic field. In the continuous case, provided the initial condition $B_{0}$ is solenoidal, then the use of $\lambda$ is unnecessary, see [5]. However, when discretizing with the finite element method, this solenoidal constraint is needed to be enforced explicitly and thus the additional variable is required. We also assume that the system (1.1)-(1.4) is equipped with homogeneous Dirichlet boundary conditions for velocity and the magnetic field.

Due to the coupling of the equations of the velocity and the magnetic field, developing efficient, accurate numerical methods for solving MHD system (1.1)-(1.4) remains a great challenge in computational fluid dynamics community. It is well known that time filter methods combined with leapfrog scheme are commonly used in geophysical fluid dynamics to reduce spurious oscillations to improve predictions, see e.g.[3], but these methods degrade the numerical accuracy and over damps the physical mode. A successfully tuned model was developed by Williams [29] reducing undesired numerical damping of [3] with higher order accuracy, see [2, 22, 24, 30] and references therein. On the other hand, in practice, the use of the backward Euler method is often preferred to extend a code for the steady state problem and this yields stable but inefficient time accurate transient solutions, see [10]. To improve this behavior, time filters are used to stabilize the backward Euler discretizations in [14] for the classical numerical ODE theory.

The present work extends the method of [7] tailored to MHD flows for constant time step. As it is mentioned in this study, the constant time step method is equivalent to a general second order, two step and A-stable method given in [9] and [19]. The scheme we consider is the time filtered backward Euler method, which is efficient, $\mathcal{O}(\Delta t^{2})$ and amenable to implementation in existing legacy codes. In addition, we also consider the numerical conservation of physically conserved quantities such as the energy and the helicity. It is worth noting that for ideal MHD with periodic boundary conditions, we prove both analytically and numerically the time filtered backward Euler method preserves the exact conservation of energy and helicity with the strong enforcement of the solenoidal constraints on the velocity and magnetic field. In addition, we prove the method’s velocity and magnetic field are both stable and long time stable without any time step restriction.

This paper is arranged as follows. Section 2 gathers notations and preliminary results which will be used for the analysis. In Section 3, the time filtered backward Euler method is described along with the proof of conservation properties. Section 4 presents stability and convergence analysis of for the fully discrete scheme. Numerical experiments are presented to verify theoretical results in Section 5. Finally, conclusions of the paper are given in Section 6.

2 Notation and Preliminaries

Standard notations of Lebesgue and Sobolev spaces are used throughout this paper. The inner product of $(L^{2}(\Omega))^{d}$, will be denoted by $(\cdot,\cdot)$, the norm in $(L^{2}(\Omega))^{d}$ by $\|\cdot\|$ and the norm in the Hilbert space $(H^{k}(\Omega))^{d}$ by $\|\cdot\|_{k}$. For clarity of presentation, we assume no-slip boundary conditions. We consider the classical function spaces

[TABLE]

The norm of the dual space $H^{-1}$ of $X$ is denoted by $\|\cdot\|_{-1}$. As usual, one has $X\subset L^{2}(\Omega)\subset H^{-1}(\Omega)$ with compact injection. The divergence free velocity space is given by

[TABLE]

We define the following norms for all Lebesgue measurable $w:[0,T]\to X$:

[TABLE]

In the error analysis, we use the Poincaré-Friedrichs’ inequality,

[TABLE]

for all $v\in X$, where $C_{p}$ is a constant depending only on the size of $\Omega$. The following properties of for the skew symmetric form are necessary in the analysis.

Lemma 2.1**.**

The trilinear skew-symmetric form $(u\cdot\nabla v,w)$ satisfies

[TABLE]

for all $u,v,w\in X$.

Proof.

Utilizing Hölder’s inequality, interpolation theorem, Sobolev embedding theorem and Poincaré inequality gives the stated results, see [21]. ∎

We use conforming finite element spaces based on edge to edge triangulations of $\Omega$ (with maximium element diameter $h$) by $X_{h}\subset X$ and $Q_{h}\subset Q$. In the computations, we consider the Scott-Vogelius finite element spaces for velocity-pressure and magnetic field-Lagrange multiplier pairs. It is well known that on a barycenter refinement of regular mesh, this element satisfies the discrete inf-sup condition, see [9] and the optimal approximation properties, [31]. Since Scott-Vogelius elements enforce mass conservation pointwisely for both velocity and magnetic field, e.g.

[TABLE]

it has been successfully used for multiphysics problems, see e.g.[4, 5].

Following [20], one admits the optimal approximation properties for the velocity and magnetic field.

[TABLE]

The discretely divergence-free space is defined by

[TABLE]

which is also the divergence-free subspace of $X_{h}$ when using Scott-Vogelius pair.

We also use the following space in the analysis:

[TABLE]

where ${\bf V}_{h}^{*}$ and $\|\cdot\|_{{\bf V}_{h}^{*}}$ are dual spaces of ${\bf V}_{h}$ and its norm which is given by

[TABLE]

The following discrete Gronwall lemma, stated in [17] plays an important role in the error analysis.

Lemma 2.2**.**

[Discrete Gronwall Lemma] Let $\Delta t$, M, and $\alpha_{n},\beta_{n},\xi_{n},\delta_{n}$ (for integers $n\geq 0$) be finite nonnegative numbers such that

[TABLE]

Suppose $\Delta t\delta_{n}<1$ for all $n$, then

[TABLE]

3 Time Filtered MHD Equations

Time filtered finite element scheme we study consists of two steps. In the first step the usual backward Euler method is applied to MHD equations and the second step includes the linear combination of solutions at previous time levels. As we will prove later, while the second step doesn’t require additional function evaluations, it has a profound impact on the solution quality such that it increases time accuracy. By assuming the prescribed values are nodal interpolants of the fluid variables, we now present the finite element approximation of (1.1)-(1.4) for constant time step method. Let $T$ denote the final time, $M$ denote the number of time steps to take and define the time step ${\Delta T}=\frac{T}{M}$. The fully discrete solution at time $t_{n}=n\Delta t$, $n=0,1,2...M$ will be denoted by $u^{n}_{h}$ and $B^{n}_{h}$. The scheme applied to the problem (1.1)-(1.4) reads as follows:

Algorithm 3.1**.**

*Given $u_{h}^{n-1},u_{h}^{n},B_{h}^{n-1},B_{h}^{n},P_{h}^{n-1},P_{h}^{n},\lambda_{h}^{n-1},\lambda_{h}^{n}$, find $(u_{h}^{n+1},B_{h}^{n+1},P_{h}^{n+1},\\ \lambda_{h}^{n+1})\in(X_{h},X_{h},Q_{h},Q_{h})$ satisfying

Step 1:*

[TABLE]

Step 2:**

[TABLE]

*for all $({v}_{h},\chi_{h},q_{h},r_{h})\in({X}_{h},{X}_{h},Q_{h},Q_{h})$. *

Step 1, without Step 2, is the classical backward Euler scheme for MHD equations analyzed in [5]. The numerical efficiency of the method is obvious. Step 2 is just an application of time filters as a modular step and its implementation is easy.

By using the following operator, Step 2 can be embedded into Step 1 in the following way. Define the interpolation operator $\mathcal{F}$ as

[TABLE]

which is formally $\mathcal{F}[{w}_{h}^{n+1}]={w}_{h}^{n+1}+O(\Delta t^{2})$. Note that reorganizing (3.5) gives $\tilde{{u}}_{h}^{n+1}=\frac{3}{2}{{u}}_{h}^{n+1}-{u}_{h}^{n}+\frac{1}{2}{u}_{h}^{n-1}$. If one repeats same calculations for the other variables, inserting all of them in (3.1)-(3.4) along with (3.9) gives

[TABLE]

for all $({v}_{h},\chi_{h},q_{h},r_{h})\in({X}_{h},{X}_{h},Q_{h},Q_{h})$. Naturally, the formulations (3.1)-(3.4) and (3.10)-(3.13) are equivalent. For simplicity of analysis, the equivalent formulation (3.10)-(3.13) of the method will be used for the complete stability and convergence analysis. However, the utilization of the method for computer simulations will be based on (3.1)-(3.4).

In the analysis, we use the following $G$-norm and $F$-norm. In general since $G$- stability implies $A$-stability, the use of $G$-matrix is very common in BDF2 analysis, see e.g.,[15] and references therein. These norms and properties are already have been given in [19]. With respect to notation of [19],(see page 392), analysis of the described method here uses the choices of $\theta=1$ and $\nu=2\epsilon$,

[TABLE]

and the $G$-norm is given by

[TABLE]

which can be negative. Here $\begin{bmatrix}{u}\\ {v}\end{bmatrix}$ is a $2n$ vector.

We also consider $F=3I^{n}$ $\in$ $\mathbb{R}^{n\times n}$ symmetric positive matrix in general case, see [19] for details. For any ${u}\in\mathbb{R}^{n}$, define $F$ norm of the $n$ vector ${u}$ by

[TABLE]

The following properties of $G$-norm are well known and for a detailed derivation of these estimations, the reader is referred to [15, 19].

Lemma 3.1**.**

$L^{2}$* norm and $G$-norm are equivalent in the following sense: there exist constants $C_{1}>C_{2}>0$ such that*

[TABLE]

Lemma 3.2**.**

The symmetric positive matrix $F$ $\in$ $\mathbb{R}^{n\times n}$ and the symmetric matrix $G$ $\in$ $\mathbb{R}^{2n\times 2n}$ satisfy the following equality:

[TABLE]

Lemma 3.3**.**

For any ${u},{v}\in\mathbb{R}^{n}$, we have

[TABLE]

[TABLE]

Proof.

Letting $\theta=1$ and $\nu=2\epsilon$ in Lemma 3.1 on p. 392 of [19] gives the stated result. ∎

The following consistency error estimations are required in the analysis.

Lemma 3.4**.**

There exists $C>0$ such that

[TABLE]

Proof.

Utilizing integral version of Taylor’s theorem, we have

[TABLE]

Hence, we get

[TABLE]

In a similar manner, one gets

[TABLE]

Integrating (3.23) and (3.24) with respect to $x$ produces

[TABLE]

Multiplying by $\Delta t$ and summing from $1$ to $N-1$ gives (3.20) and (3.21). ∎

3.1 Conservation Laws

We study conservation properties of the scheme (3.10)-(3.13). Energy and helicity are very important flow quantities and play an important role in flow’s structures, [23]. It is well known that, an accurate model must predict these quantities correctly to verify the physical fidelity of the model. We now show the time filtered backward Euler (3.10)-(3.13) is an energy and helicity preserving scheme.

Lemma 3.5**.**

(Global Energy Conservation) Scheme (3.10)-(3.13) satisfies the energy equality:

[TABLE]

Proof.

Set $v_{h}=\mathcal{F}[{u}_{h}^{n+1}]$ in (3.10), $q_{h}=\mathcal{F}[P_{h}^{n+1}]$ in (3.11) , $\chi_{h}=\mathcal{F}[{B}_{h}^{n+1}]$ in (3.12) and $r_{h}=\mathcal{F}[\lambda_{h}^{n+1}]$ in (3.13), then the trilinear terms $(\mathcal{F}[{u}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{u}_{h}^{n+1}]),\mathcal{F}[{u}_{h}^{n+1}])$ and $(\mathcal{F}[{u}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{B}_{h}^{n+1}]),\mathcal{F}[{B}_{h}^{n+1}])$, the pressure term and the $\lambda$ term vanish by the use of (2). Then, one gets

[TABLE]

Note that since $\nabla\cdot\mathcal{F}[B_{h}^{n+1}]=0$, we get $(\mathcal{F}[{B}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{B}_{h}^{n+1}]),\mathcal{F}[{u}_{h}^{n+1}])=-(\mathcal{F}[{B}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{u}_{h}^{n+1}]),\mathcal{F}[{B}_{h}^{n+1}])$. Multiplying (3.29) by $s$ and adding (3.28) to (3.29) produces

[TABLE]

Reorganizing (3.30) by using Lemma 3.2 and multiplying with $\Delta t$ yields

[TABLE]

Summing (3.31) from $n=1$ to $n=N-1$ gives the stated energy result. ∎

Lemma 3.6**.**

(Global Cross Helicity Conservation) Scheme (3.10)-(3.13) satisfies the cross helicity equality:

[TABLE]

where $I[{w}_{h}^{n+1}]={w}_{h}^{n+1}-2{w}_{h}^{n}+{w}_{h}^{n-1}$

Proof.

To prove the global cross helicity conservation, set $v_{h}=\mathcal{F}[{B}_{h}^{n+1}]$, $q_{h}=\mathcal{F}[\lambda_{h}^{n+1}]$, $\chi_{h}=\mathcal{F}[{u}_{h}^{n+1}]$, $r_{h}=\mathcal{F}[P_{h}^{n+1}]$ in (3.10)- (3.13), respectively. Since the trilinear terms $(\mathcal{F}[{B}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{B}_{h}^{n+1}]),\mathcal{F}[{B}_{h}^{n+1}])$, $(\mathcal{F}[{B}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{u}_{h}^{n+1}]),\mathcal{F}[{u}_{h}^{n+1}])$, the pressure term and the $\lambda$ term vanish by the use of (2) , one has

[TABLE]

Note that since $(\mathcal{F}[{u}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{u}_{h}^{n+1}]),\mathcal{F}[{B}_{h}^{n+1}])=-(\mathcal{F}[{u}_{h}^{n+1}]\cdot\nabla(\mathcal{F}[{B}_{h}^{n+1}]),\mathcal{F}[{u}_{h}^{n+1}])$ and

[TABLE]

where $I[{w}_{h}^{n+1}]={w}_{h}^{n+1}-2{w}_{h}^{n}+{w}_{h}^{n-1}$, adding (3.33) and (3.34) yields

[TABLE]

Summing (3.36) from $n=1$ to $n=N-1$ and multiplying by $\Delta t$ produces the cross helicity conservation result. ∎

4 Convergence Analysis

4.1 Stability and Long Time Stability

This section presents unconditional stability, long time stability and convergence analysis of the proposed method.

Lemma 4.1**.**

Let $f\in L^{2}(0,T;H^{-1}(\Omega))$ and $g\in L^{2}(0,T;L^{2}(\Omega))$. Then, solutions to the scheme (3.10)-(3.13) are unconditionally stable, and satisfy the following bounds at $t_{M}=M\Delta t$

[TABLE]

Proof.

The proof starts with using the global energy conservation (3.27).

[TABLE]

The forcing terms can be bounded by using Cauchy-Schwarz and Young’s inequalities as

[TABLE]

Inserting (4.3) and (4.4) in (4.2) gives

[TABLE]

Using Lemma 3.3, we get

[TABLE]

Multiplying (4.6) by $\frac{4}{3}$ and using induction finishes the proof. ∎

We also show that the scheme is unconditionally long time stable.

Lemma 4.2**.**

Let ${f},{g}\in L^{\infty}(\mathbb{R}_{+},{\bf V}_{h}^{*})$, then the approximation (3.10)-(3.13) is long time stable in the following sense: for any $\Delta t>0$

[TABLE]

where $\omega=(1+\alpha)(1+\beta)$, $\alpha=\min\{\dfrac{C_{1}^{2}Re^{-1}\Delta t}{8C_{p}^{2}},2\}$,$\beta=\min\{\dfrac{C_{1}^{2}Re_{m}^{-1}\Delta t}{8C_{p}^{2}},2\}$, $C_{1}$ is given by Lemma 3.1 and $C_{p}$ is given by (2.2).

Proof.

Applying Cauchy-Schwarz and Young’s inequalities for the global energy conservation equation (3.27) yields

[TABLE]

Dropping the fourth term in the left hand side of (4.8), and adding both sides $\dfrac{Re^{-1}\Delta t}{8}\|\nabla\mathcal{F}[{u}_{h}^{n}]\|^{2}$ and $\dfrac{sRe_{m}^{-1}\Delta t}{8}\|\nabla\mathcal{F}[{B}_{h}^{n}]\|^{2}$ results in

[TABLE]

The third and fourth terms can be bounded by using Poincaré’s-Friedrichs’ inequality and Lemma 3.1 as

[TABLE]

where $\alpha=\min\{\dfrac{C_{1}^{2}Re^{-1}\Delta t}{8C_{p}^{2}},2\}$. Using similar techniques for the fifth and the sixth terms, we get

[TABLE]

where $\beta=\min\{\dfrac{C_{1}^{2}Re_{m}^{-1}\Delta t}{8C_{p}^{2}},2\}$. Inserting (4.10)-(4.11) in (4.9) and using induction yields the stated result.

∎

4.2 A-priori Error Estimate

In this section, we present a detailed convergence analysis of the proposed time filtered method for MHD equations. We define the discrete norms as

[TABLE]

For the optimal asymptotic error estimation, we assume the following regularity assumptions for the exact solution of (1.1)-(1.4):

[TABLE]

The mesh and velocity approximating polynomial degree $k$ is chosen so that the Scott-Vogelius pair is inf-sup stable and the properties (2.8)-(2.9) hold.

Theorem 4.1**.**

Suppose regularity assumptions (4.13) hold. Under the following time step condition

[TABLE]

there exists a positive constant $C^{*}$ independent of $h$ and $\Delta t$ such that

[TABLE]

Proof.

The proof starts by deriving the error equations. We consider continuous variational formulations of (1.1)-(1.4) at the time level $t=t^{n+1}$. Adding and subtracting terms yields the following variational formulations for the velocity,

[TABLE]

for all $v_{h}\in V_{h}$ and for the magnetic field

[TABLE]

for all $\chi_{h}\in V_{h}$ where

[TABLE]

Denote the error between finite element solution and continuous solution by $e_{u}^{n}:={u}^{n}-u_{h}^{n}$ and $e^{n}_{B}:={B}^{n}-B_{h}^{n}$. The error equations are obtained by subtracting (3.10)-(3.12) from (4.16)-(4.17), respectively:

[TABLE]

and

[TABLE]

We split the errors as follows

[TABLE]

where $I_{{u}}^{n}$ and $I_{{B}}^{n}$ are the interpolations of $u^{n}$ and $B^{n}$ in $V_{h}$, respectively. Substituting (4.22) into (4.20) and (4.23) into (4.21), choosing $v_{h}=\mathcal{F}[\phi_{{u},h}^{n+1}]$ and using Lemma 3.2 and (2), leads to

[TABLE]

Then, we now bound the terms in the right hand side of (4.24) and obtain

[TABLE]

for the first two terms along with the Cauchy-Schwarz and Young’s inequalities. Also, with Lemma 2.1, Cauchy-Schwarz and Young’s inequalities, we get estimations for the nonlinear terms:

[TABLE]

In addition, the terms in consistency error $|E_{1}({u},{B},\mathcal{F}[\phi_{{u},h}^{n+1}])|$ are bounded by using Cauchy-Schwarz, Poincarè and Young’s inequalities as follows:

[TABLE]

Inserting (4.25)-(4.36) into (4.24) yields

[TABLE]

Multiplying (4.37) by $\Delta t$ and summing from $t=1$ to $t=N-1$, we have

[TABLE]

Using Lemma 3.4 and approximation properties (2.8)-(2.9), we have

[TABLE]

Substituting (4.39)-(4.44) in (4.38) and utilizing Lemma 4.1, one gets

[TABLE]

Reorganizing equation (4.45), we have

[TABLE]

In a similar manner, substituting (4.23) into (4.21), and setting $\chi_{h}=\mathcal{F}[\phi_{{B},h}^{n+1}]$ gives

[TABLE]

Multiplying (4.47) by $s$, adding it to (4.46) and using that

[TABLE]

we get

[TABLE]

Application of the discrete Gronwall inequality with

[TABLE]

and utilization of Lemma 3.3 yields

[TABLE]

Multiplying (4.50) with $\frac{4}{3}$ and applying induction produces

[TABLE]

The proof is completed by applying the triangle inequality. ∎

5 Numerical Studies

In this section, Algorithm 3.1 presented in Section 3 will be studied at examples given in a two-dimensional domain $\Omega$. We perform three different numerical tests in order to expose the promise of proposed method. The first example has been designed to confirm the theoretically predicted results of Theorem 4.1. In the second test, we check the energy and cross-helicity conservation properties of the scheme for an ideal MHD case. In the final test, we investigate the flow behavior in a channel over a step under the effect of magnetic field. The initial velocity, the initial pressure and the initial magnetic field were computed as nodal interpolants if not stated otherwise. For all simulations, the Scott-Vogelious pair of finite elements $((P_{2})^{2},P_{1}^{disc})$ on barycenter refined triangular meshes is used. The computations were performed with the public license finite element software FreeFem++ [16].

5.1 Convergence Rate Verification

We consider the MHD equation (1.1)-(1.4) in the unit square and in the time interval $[0,1]$ where the right hand side and the boundary conditions are chosen such that

[TABLE]

is the solution. Other problem parameters are chosen as $Re=Re_{m}=s=1$. Since we are studying convergence, the spatial meshwidth $h$ and the time step $\Delta t$ are set to be same in order to see the errors and rates at once. We measure the errors in the discrete norm $L^{2}(0,T;{H}^{1}(\Omega))$ for the velocity and the magnetic field which could be written for the velocity for example:

[TABLE]

Table 1 reports the order of convergence for Algorithm 3.1. One can observe the predicted second order convergence for the errors estimated in Theorem 4.1.

We note that this test was also carried out for both with filtering and not filtering the pressure. More precisely, $\tilde{P}_{h}^{n+1}$ in (3.1) and $\tilde{\lambda}_{h}^{n+1}$ in (3.3) are chosen as $P_{h}^{n+1}$ and $\lambda_{h}^{n+1}$ and not updated in Step $2$ of Algorithm 3.1. In both ways, we obtain the same error rates showing that the pressure filtering does not affect the velocity and magnetic field solution, exactly the same situation for Navier-Stokes equations, [7].

5.2 Orszag-Tang Vortex Test

As a second numerical test, we solve Orszag-Tang vortex problem which is a well-known model for testing MHD codes. Due to the complex interaction between various shock waves traveling at different speed regimes, this problem tests robustness of the code in the formation of shocks and shock-shock interactions in the ideal MHD case, (see [8, 26] and references therein). In addition, since the numerical solution of Orzag-Tang vortex system does not necessarily preserve the incompressible constraint $\nabla\cdot B=0$, this problem also provides some quantitative estimations for the effect of significant magnetic monopoles on the numerical solutions. In this test problem, our goal is to show the confirmation of the conserved quantities and compare the results with unfiltered case in order to see the effect of the time filter explained in detail in Section 3.1. The problem is solved in $[0,2\pi]\times[0,2\pi]$ using the meshwidth $h=1/32$, the time step $\Delta t=0.01$ and the final time $2.7$. For an ideal MHD case, the selected parameter choices are $Re=Re_{m}=\infty$, $s=1$ and $f=\nabla\times g=0$. Consider the following initial conditions

[TABLE]

along with the periodic boundary conditions. Since an ideal MHD case is assumed, the global energy and the cross helicity defined by

[TABLE]

should be conserved through the solutions obtained by Algorithm 3.1. As depicted in Figure 1, the quantities of interest are exactly conserved, while the backward Euler scheme which consists of only discarding the filters fails to preserve them. We can deduce that the classical backward Euler method ruins the energy and cross helicity properties and time filters correct this behavior. Thus, the results for conserved quantities are consistent with the theory.

It is worth noting that the simulations are ran using the coarse mesh which already provides very similar results as the finest grid $4096\times 4096$ of [8] and $1024\times 1024$ of [26].

5.3 MHD Channel Flow Over a Step

Our final numerical example is to test Algorithm 3.1 for the benchmark MHD channel flow over past a step. The problem geometry consists of a rectangular $40\times 10$ channel with a $1\times 1$ step places $5$ units into the channel at the bottom. We pick $Re=1000$ and $Re_{m}=1$ along with varying $s$ and Dirichlet boundary conditions corresponding to no slip velocity on the walls. We impose ${u}=\left(y(10-y)/25,0\right)^{T}$ for the velocity on the inlet and outlet and ${u}=0$ for the rest. For the magnetic field boundary condition, we take $B=\left(0,1\right)^{T}$ on all boundaries. As initial conditions, we take ${u}=\left(y(10-y)/25,0\right)^{T}$ and $B=0$. The computations are carried out with $\Delta t=0.025$ up to an end time of $40$ that provides $328,148$ total degrees of freedom. The development of the flow is depicted in Figure 2. Our interest is only flow behaviour behind the step, thus we present the figures up to $30\times 10$ part of the channel. Note that since there is no magnetic force in the case of $s=0$, we only give velocity streamlines over speed contours. For $s=0.01$, two eddies start to develop behind the step and the eddies separate from the step between $t=5$ and $t=10$. Due to the effect of the Lorentz force, the peeling of the eddies behind the step is suppressed for $s=0.05$. As result, the solution captures the correct eddy formation and detachment behind the step. We note that the initial parabolic profile of the initial velocity is changed and the results shown in Figure 2 are compatible with [1].

6 Conclusions

An efficient time filtered method as a post processing step is introduced to MHD equations in a given backward Euler code. We have shown that the time filtered algorithm increases accuracy from first order to second order without any extra programming effort. We have provided a complete numerical analysis of the method, including unconditional, long time stability and optimal convergence rates. Moreover, time filtered backward Euler method conserves energy and cross helicity when the solenoidal constraints on the velocity and magnetic field enforced strongly. Results of several numerical tests have been presented in order to verify all theoretical findings. The numerical investigations have shown the time filtered backward Euler method to be very effective and to predict the energy and helicity very well in comparison with the backward Euler method.

Several research directions will be pursued in future. For instance, we will study variable time step methods for MHD, which require only one BDF solve at each time level followed by addition of the solution at previous time steps to develop embedded family of higher order accuracy. In addition, the extension of time filtering to more complex coupled flow problems such as MHD convection and double diffusive convection will be topics of future research.