A Modular Grad-Div Stabilization Method for Time-Dependent Thermally Coupled MHD Equations

In this paper, we consider a fully discrete modular grad-div stabilization algorithm for time-dependent thermally coupled magnetohydrodynamic (MHD) equations. The main idea of the proposed algorithm is to add an extra minimally intrusive module to penalize the divergence errors of velocity and improve the computational efficiency for increasing values of the Reynolds number and grad-div stabilization parameters. In addition, we provide the unconditional stability and optimal convergence analysis of this algorithm. Finally, several numerical experiments are performed and further indicated these advantages over the algorithm without grad-div stabilization.


Introduction
The incompressible magnetohydrodynamic (MHD) model has a wide range of applications in scientific and engineering, such as electromagnetic pumping, liquid metal, electrolyte, and so on (see [1][2][3][4]). As we know, the MHD model describes the interaction of an incompressible viscous conducting fluid and the electromagnetic field. In other words, it is a multi-physics phenomenon: the magnetic field changes the momentum of the fluid through the Lorenz force, and conversely, the conducting fluid influences the magnetic field through electric currents. When the buoyancy effects cannot be neglected in the momentum equation (owing to temperature differences in the flow), the MHD equations are usually coupled to the heat equation.
In the recent years, much effort has been spent on the development of some efficient numerical methods to investigate this problem. Meir considered the existence and uniqueness of solutions for the thermally coupled MHD flow in [6], and developed the Galerkin finite-element method (FEM). Moreover, optimal error analysis of the model was established in [7]. A stabilized finite-element method was proposed in [8]. Furthermore, a decoupled Crank-Nicolson time-stepping scheme and partitioned time-stepping scheme for the thermally coupled MHD system were considered in [9,10], respectively, and some meaningful stability and convergence results were presented. Yang and Zhang [11] gave the convergence and stability analysis of three iterative methods of the steady thermally coupled MHD equations. In addition, Ding et al. studied convergence analysis of the Crank-Nicolson-extrapolated fully discrete scheme [5] and gave a fully discrete Euler semi-implicit scheme with the magnetic equation approximated by Nédélec edge elements to capture the physical solutions [12], respectively. In addition, the modified characteristics finite-element method and projection method have been proposed in [13,14]. The unconditional stability of the fully discrete scheme and the optimal second-order convergent accuracy in both time and spatial discretizations were proved in [15]. Moreover, a linear fully decoupled velocity correction method for the thermally coupled MHD model was studied in [16].
It is worth mentioning that classical conforming finite-element discretizations for incompressible flows relax the divergence constraint, and give only a relatively weak limit. Although this enables us to construct a stable discretization of the inf-sup condition, a weak limit will lead to errors in the continuous pressure that depends on the Reynolds number and causes inaccurate computational solutions for many flow problems. In order to overcome this difficulty, a grad-div stabilization was discovered for the first time in [17], which is a simple and popular method for improving mass conservation of numerical solutions and only adds a term that equals zero in a continuous equation. The analysis of grad-div stabilization for Stokes equations and Navier-Stokes equations were proposed in [18,19]. For the time-dependent Stokes/Darcy model, two grad-div stabilization methods were proposed in [20]. In addition, a grad-div stabilized projection finite-element method for a double-diffusive natural convection model was given in [21]. In view of this, a great deal of related interesting works have been reported in the recent years [22][23][24]. Although it is easy to implement it in the program, it also has some shortcomings: on the one hand, the stabilization makes it too difficult to solve due to an increased coupling in the system. On the other hand, too large values of grad-div parameters cause a low condition number of the corresponding linear system [25].
To solve the above problems, recently, a modular grad-div stabilization method has been proposed in [26], which allows the Navier-Stokes equations to be solved in two steps. Then, they also gave a BDF2 modular grad-div stabilization method for the Navier-Stokes equations in [27]. Next, the modular grad-div stabilization method of MHD and Boussinesq equations was proposed in [28,29], respectively. In [30], Li et al. presented a rotational pressure-correction method for the Stokes/Darcy model based on the modular grad-div stabilization.
As we know, the conservation of mass plays an important role in the construction of numerical schemes for incompressible fluids. Under certain extreme situations, nonphysical phenomena may appear if discrete solutions are not mass-conservative (see [31] for comprehensive discussions). Moreover, a large Reynolds number will cause the problem of convection dominance, which makes it very difficult to solve. Recently, many researchers are interested in studying highly efficient numerical algorithms for the thermally coupled MHD problem, but there is less attention to deal with the conservation of mass and high Reynolds number simultaneously. Since the modular grad-div stabilization not only improves the conservation of mass, but has been proven useful for increasing values of the Reynolds number and grad-div parameters, the purpose of this paper is to apply the idea from [27] to the thermally coupled MHD model to improve the mass conservation of numerical solutions, and guarantee that the proposed algorithm is still effective for large Reynolds numbers. Here, we propose a first-order fully discrete modular grad-div stabilization method for the thermally coupled MHD equations, which adopts an Euler semi-implicit scheme for the time-discretization. This scheme is divided into two steps; in the first step, the intermediate velocity and other physical unknown quantities are solved. In the second step, we add two penalty terms to enforce improving the mass conservation and ensure high efficiency of the algorithm with large Reynolds numbers and grad-div stabilization parameters. Moreover, the unconditional stability and error estimation corresponding to this scheme are completed in this paper. Numerical examples further verify the reliability of the proposed algorithm.
This paper is arranged as follows. Section 2 describes the necessary notations and mathematical preliminaries. In Section 3, a fully-discrete modular grad-div stabilization method for the incompressible thermally coupled MHD equation is presented. In Sections 4 and 5, we give its complete stability and convergence analysis, respectively. We also present some numerical experiments to illustrate the reliability and effectiveness of the method in Section 6. Finally, the last section summarizes the results of the paper.

Preliminaries
In this section, we use · 0,2 and (·, ·) to denote the usual L 2 (Ω) and its inner product. For m ∈ N + , 1 ≤ p ≤ ∞, the L p (Ω) norm and W m,p (Ω) norm are denoted by · 0,p and · m,p , respectively. Particularly, H m (Ω) represents the case of p = 2. In addition, X is defined as a normed function space in Ω, L p (0, T; X) is the space of all functions defined on Ω × (0, T), and the norm is bounded These notations of Lebesgue and Sobolev spaces are used throughout this paper. We consider the classical function spaces as follows: the divergence-free subspaces of X and M are defined by: From [32], we have the following two formulas , which imply that for all Φ, Ψ ∈ M and w ∈ X, Therefore, the weak formulation of (1) reads: find (u, p, B, θ) ∈ L 2 (0, T; X) × L 2 (0, T; Q) × L 2 (0, T; M) × L 2 (0, T; W), such that for all (v, q, C, ϕ) ∈ X × Q × M × W, Furthermore, the analysis of time-discretization utilizes the following norms, for 1 ≤ m < ∞ : For u, v, w ∈ X, we have some properties of these trilinear forms (see [33]) Additionally, for v ∈ X and B, C ∈ M from [34,35], we have the following bounds: Here and after, C (with or without a subscript) denotes a general positive constant, which may represent different values in different situations. In addition, we need the following assumptions on the prescribed data for problem (1), as it is useful for our later theoretical analysis.
Assumption 3. Assume that the boundary of Ω is smooth so that the unique solution (v, q) of the Stokes problem in [34] −∆u + ∇q = f u , divv = 0 in Ω, v| ∂Ω = 0, and Maxwell's equations The following lemma is very important in convergence analysis; so, we recall it from [34].
Lemma 1 (Discrete Gronwall's Lemma). Let a n , b n and d n for the integer n ≥ 0 be nonnegative numbers such that For the spatial discretization, we define the following finite spaces, where (u, p) using the finite-element pair (P b 1 , P 1 ) and B, θ using P 1 element. Let π h = {K} be a uniformly regular family of triangulation of Ω, and define the mesh size h = max Here, P b 1 is defined as (more details see [32]) Furthermore, we need the subspace X 0h of X h which is defined as The projections satisfy the following properties (see [32]): Next, we give the numerical scheme of this paper.

A Modular Grad-Div Stabilization Method for Time-Dependent Thermally Coupled MHD Equations
Now, we construct a fully-discrete numerical scheme for solving the model system (1) and prove the corresponding unconditional stability. Divide the simulation time T into , where t n = nτ, T = Nτ. Our numerical scheme reads as follows.

Stability Analysis
Now, we focus on the stability of Algorithm 1. Our stability analysis shows that approximate solutions of Algorithm 1 are stable without any time step restriction. In order to obtain the stability result, we first present a lemma which gives a relation between solutions of Step 1 and Step 2.

Error Analysis
In this section, we show that solutions of the proposed algorithm converge to the true solutions of (1). In order to obtain the equations, we denote true solutions at time level t n+1 , the error analysis needs the following error decomposition at time level t n+1 :

Lemma 3.
Consider the second step of Algorithm 1, then it holds (see [29]) Theorem 2. Suppose that Assumption 1-3 are satisfied, then the following estimate holds Proof. For simplicity, our entire proof process is divided into three steps, as shown below.
Step 1: [The derivation of error equations] (2) with t = t n+1 , and use integration by parts to get Then, subtracting (5) from (9)-(11), respectively, we obtain Here, Using the decomposition and setting v h = 2τΛ n+1 u,h in (12), (14), and adding the three equations, we deduce that . Then, the above equation can be rewritten as follows: Step 2: [The estimation of the right-hand side terms of error equations] We now estimate each term of the right-hand sides of (15) separately. Next, using (3), the inverse inequality, the Canchy-Schwarz and Young's inequalities, we arrive at Then, using (4) with Cauchy-Schwarz and Young's inequalities, we obtain Similarly, we obtain Showing that Further, we can obtain E n+1 For the above three inequalities, we sum n = 0 to N − 1 Based on Assumptions 1 and 3, we arrive at Step 3: [The completion of the proof] Substituting all the above inequalities into the right-hand term of (15) and applying Assumptions 1-3, asserting Lemma 3, yields Then, we sum over time steps and use the results of Lemma 2 and regularity Assumptions 1-3, and arrive at Further, we apply Gronwall's lemma and use φ 0 Finally, applying the triangle inequality yields the desired results.

Numerical Experiment
In this section, we give two numerical experiments to illustrate the reliability of the thermally coupled MHD problem. One aspect is to verify the predicted stability and convergence rates of the previous section, another is to consider the flexibility of large Reynolds number and grad-div stabilization parameters. In the following tests, we use the finiteelement pair (P b 1 , P 1 , P 1 , P 1 ) for the velocity/pressure/magnetic/temperature, respectively.

An Exact Solution Problem
First, we consider an exact solution problem to verify the stability and the convergence rates of Algorithm 1 for problems (1). Let the domain Ω = [0, 1] × [0, 1] ∈ R 2 and the mesh is obtained by dividing Ω into squares and drawing a diagonal in each square. Therefore, the exact solution (u, p, B, θ) is given by then, the external force terms f 1 and f 2 , boundary conditions, and initial values in Equations (1) are selected to correspond to the exact solution. The parameters κ = S = R m = 1, end time T = 1.
In Table 1, we set β 0 = 0.2, γ 0 = 1, ν = 1 for convergence rates and vary the mesh size 1/h between 4, 8, 16, 32, 64. The expected accuracy is consistent with theoretical results. In addition, we can get the H 1 (Ω)-norm convergence rates of velocity, magnetic, and temperature fields to be O(h), the L 2 (Ω)-norm convergence rates of velocity magnetic, and temperature fields are O(h 2 ). Table 1. Errors and convergence rates of the considered scheme with τ = h 2 at T = 1. Next, our test is for Re increasing. We fix τ = h, h = 1/32 and set β 0 = 0.2, γ 0 = 1. In Table 2, errors for velocity and pressure with increasing Re of the method without grad-div and modular grad-div stabilization method are compared. The corresponding solutions are no-Stab and Modular. We observe that the error of the proposed algorithm hardly increases. However, the approximate solutions generated by the no-Stab method is getting worse and worse, especially for gradient and divergence results of velocity. Finally, we fix Re = 1, τ = h, h = 1/32, vary the grad-div parameters 0.1 ≤ β 0 ≤ 10 5 and 0.1 ≤ γ 0 ≤ 10 5 . The results are presented in Table 3. We observe that the result of velocity divergence errors of our method becomes small as γ 0 increases, but β 0 has not much impact on them. Table 3. Velocity errors and divergence of the modular grad-div methods with different β 0 , γ 0 .

Thermally Driven Cavity Flow Problem
In this experiment, we consider thermal driven cavity flow [36]  Here, we set h = 1/64, τ = 0.01, S = R m = κ = β 0 = γ 0 = 1. As shown in Figures 1 and 2, we plot streamlines of velocity and magnetic, isotherms of temperature of methods no-Stab and our proposed method at different times when the Reynolds number is Re = 1. We can observe that the results of our proposed method are consistent with no-Stab method, which verifies its correctness.
Moreover, in Figure 3, we show streamlines of velocity and magnetic, isotherms of temperature of our method at different times when the Reynolds number is Re = 10 6 . As time goes on, the vortex of the velocity streamline gradually moves from left to right, the streamline of the magnetic field gradually becomes curved in Figure 3 (Re = 10 6 ). At the same time, due to the temperature difference between the left and right walls, the isotherm also gradually becomes curved. However, these results do not change significantly when Re = 1, as can be seen from Figures 1 and 2. It is worth noting that the no-Stab method is divergent with Re = 10 6 when T = 3 s, but the modular grad-div stabilization method is still convergent in this case. So, we can find that the proposed method is efficient for the thermally driven cavity flow problem with relatively large Reynolds numbers.

Conclusions
We developed a first-order fully discrete modular grad-div stabilization algorithm for time-dependent thermally coupled MHD equations. The advantages of this scheme is to keep the conservation of mass as much as possible and its effectiveness with high Reynolds number and large grad-div stabilization parameters. Then, the scheme is proven to be stable and convergent. When compared without grad-div stabilization solutions, our algorithm exhibits a smaller divergence error of the velocity and shows how β 0 and γ 0 influence this effect. Moreover, we also confirm that the scheme still maintains the advantage with high Re and grad-div stabilization parameters. In the future, we will consider high-order schemes with modular grad-div stabilization.

Acknowledgments:
The authors would like to thank the editor and referees for their valuable comments and suggestions, which helped us to improve the results of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.