Linear Full Decoupling, Velocity Correction Method for Unsteady Thermally Coupled Incompressible Magneto-Hydrodynamic Equations

We propose and analyze an effective decoupling algorithm for unsteady thermally coupled magneto-hydrodynamic equations in this paper. The proposed method is a first-order velocity correction projection algorithms in time marching, including standard velocity correction and rotation velocity correction, which can completely decouple all variables in the model. Meanwhile, the schemes are not only linear and only need to solve a series of linear partial differential equations with constant coefficients at each time step, but also the standard velocity correction algorithm can produce the Neumann boundary condition for the pressure field, but the rotational velocity correction algorithm can produce the consistent boundary which improve the accuracy of the pressure field. Thus, improving our computational efficiency. Then, we give the energy stability of the algorithms and give a detailed proofs. The key idea to establish the stability results of the rotation velocity correction algorithm is to transform the rotation term into a telescopic symmetric form by means of the Gauge–Uzawa formula. Finally, numerical experiments show that the rotation velocity correction projection algorithm is efficient to solve the thermally coupled magneto-hydrodynamic equations.


Introduction
The unsteady incompressible MHD equation is studied in this paper. The buoyancy effect cannot be ignored in the momentum equation due to the temperature difference of fluid flow [1][2][3]. Therefore, we consider unsteady incompressible thermally coupled MHD system. This is a strongly coupled model through the famous Boussinesq approximation [1,2,4,5]. The model is a multi-physics phenomenon: first of all, the movement of the conductor under the presence of a magnetic field generates an electric current, which changes the existing electromagnetic field; then, current and magnetic field produce Lorenz forces, which accelerate fluid particles along the lines of magnetic and current; finally, incompressible MHD are usually coupled with thermal equations because of the temperature difference between the conductive current in the momentum equation and the inability to ignore the buoyancy effect. In this way, velocity, pressure, magnetic induction, and temperature in multiple physical fields are coupled. Thermally coupled incompressible MHD model has been widely used in industries and engineering such as magnetic propulsion devices, nuclear reactor technology, semiconductor manufacturing, metal hardening, casting, melting and crystal growth; see [1,3,[6][7][8][9].
In this paper, we consider the following unsteady thermally coupling incompressible MHD equations [1,2,4,5]: where Ω is a bounded polygonal domain in R d (d = 2, 3) and ∂Ω is a polygon boundary, therefore, our main solution region is as shown in Figure 1, T is the final time. The functions of υ denotes the velocity field, p denotes the pressure, b denotes the magnetic field and ϑ denotes the temperature. f is a forcing term for the magnetic induction, g is the known applied current with ∇ · g = 0, Ψ is a given heat source. υ 0 , b 0 and ϑ 0 are the initial velocity, the initial magnetic and the initial temperature, respectively. The initial magnetic induction b 0 satisfies ∇ · b 0 = 0. n is the normal vector outside the unit of ∂Ω. For considering the heat equation, as early as 1994, the existence and uniqueness of the solution of the stationary thermally coupled MHD equations were studied in [1,10,11]. In 1999, the influence of magnetic Prandt number on fluid in a three-dimensional nonlinear convection model under strong vertical magnetic field was studied in [5]. In 2010, the steady-state thermally coupled MHD equation under two gravity models was studied, and the existence and uniqueness results of corresponding weak solutions through data under certain preconditions were studied in [12]. In 2011, the use of a stable finite element method to approximate the thermally coupled MHD problem was proposed, and a numerical formula to solve this equation was proposed in [2]. In 2018, the MHD equation from the heat equation at each time step by using a partitioning method was decoupled, refer to [3]. Recently, the convergence analysis of the thermally coupled MHD Crank-Nicolson extrapolation full-discrete scheme was proposed in [4]. A linearized projection scheme for non-stationary incompressible coupled the MHD with heat equations is introduced; meanwhile, a linearized fully discrete scheme is proposed in [13]. Although considerable work has been conducted to develop efficient schemes for both steady and unsteady thermally coupled MHD model. Little attention has been paid to the fully decoupled scheme of unsteady thermally coupled MHD model.
In order to design an efficient and stable numerical scheme for the thermally coupled MHD equations, the main difficulties we will face are: (i) strong nonlinear terms in the MHD equations and heat equation; (ii) velocity and temperature are coupled due to buoyancy effects; due to Lorentz force and Ohm's law, the velocity field is coupled with the magnetic field; (iii) due to the incompressible condition of velocity; namely, ∇ · υ = 0, the velocity and pressure are strongly coupled together, forming a saddle point problem; (iv) due to the artificial Neumann boundary condition of pressure, the numerical boundary layer will be generated, so that the L 2 -norm of pressure and the H 1 -norm of velocity cannot reach the optimal order; (v) due to the incompressible condition of the magnetic field; namely, ∇ · b = 0, this cause the model to produce singular solutions. Therefore, it is necessary to design a linear, fully decoupled and stable format for the model. To deal with (iii) difficulty encountered in the model, a common strategy to decouple the computation of the pressure from the velocity is a projection-type schemes as in the case for Navier-Stokes equations, refer to [14,15].
Projection methods can be viewed as fractional/splitting step methods, where convection diffusion and incompressibility are dealt with in two steps, refer to [14][15][16][17][18][19]. For velocity correction projection methods, sticky item is displayed processed or ignored. The pressure is made explicit in the first sub-step and is corrected in the second one by projecting the provisional velocity onto the space of incompressible vector fields. The velocity obtained in the convection-diffusion sub-step is projected in order to satisfy the weak incompressibility condition. It's well known that SVC projection methods which the projection step precedes the viscous step, one could also refer to these methods as "projection-diffusion" methods as in [20]. This method creates artificial Neumann boundary conditions for pressure which cuts down the accuracy of the pressure approximation. The RVC scheme is proposed in [18,21,22], that leads to improved pressure approximation. More importantly and appealing, using projection methods only solves a sequence of decoupled elliptic equations for the velocity and the pressure at each time step; meanwhile, it is very efficient for large scale numerical simulations. In order to decouple the velocity and pressure, we choose the SVC and RVC in the velocity correction projection methods. The RVC scheme not only solves the problem of velocity and pressure coupling in the model, but also overcomes the limitation of artificial Neumann boundary condition for the pressure, thus improving the error accuracy of pressure.
To overcome the above difficulties, we use first-order scheme marching in time of RVC projection method to solve time dependent thermally coupled MHD problem. We deal with difficulties (i) and (ii) by using implicit-explicit format processing for nonlinear terms refer to [2,4,23], so we can decouple velocity, magnetic, temperature from the model. In addition, the velocity correction projection formats are used to solve difficulties (iii)-(iv). In order to overcome the difficulty (v), H(curl, Ω) is selected for the magnetic field and H(curl)conforming Nédélect element approximation is used for the magnetic field selection, thus the singularity of solutions in non-convex regions is avoided. Therefore, the difficulties (i)-(v) are successfully solved. After the above processing, the complex nonlinear saddlepoint system is transformed into a series of simple linear elliptic problems, which greatly improves the computational efficiency. For space approximation, we use uniform finite elements: P b 1 to discrete υ and ϑ, the b is implemented by H(curl)-conforming Nédélec element, the continuous P 1 finite element for discretizing p. Besides, we provide the energy stability for the velocity correction projection methods. Last, numerical experiments verify the stability and convergence of the RVC algorithm.
The rest of this paper is organized as follow. In Section 2, we propose some notations for the magneto-thermal coupling model. In Section 3, we put forward a linear full decoupling velocity correction algorithms for system (1). Meanwhile, we give the corresponding stability of the proposed algorithms and derive a detailed proofs. In Section 4, to further verify the stability and effectiveness of the considered model and RVC algorithm we conduct corresponding numerical experiments, and fully demonstrate the advantages of the RVC algorithm. Finally, we summarize the conclusions of the article and make an outlook for future research work in Section 5.

Functional Setting for the Magneto-Thermal Coupling Model
To begin with, the following spaces are defined by: In references [24,25], we know the following Poincaré inequalities and embedded inequality for C is a generic coefficient, where the space L 2 (Ω) is equipped with the inner product (·, ·) and L 2 -norm · .

Linear Full Decoupling Algorithms and Their Stabilities
In the section, firstly, we present linear fully decoupling SVC scheme and stability for thermally coupled incompressible MHD equations. However, SVC scheme can cause obviously an artificial Neumann boundary condition for the pressure, which cuts down the accuracy of the pressure approximation. Therefore, in order to deal with the limitation of artificial Neumann boundary conditions for the pressure, we propose an RVC scheme for this model and give its stability.

Remark 1.
In fact, the solution for the magnetic introduction still satisfies the weakly divergence free property. Since ∇ϕ ∈ W for all ϕ ∈ H 1 0 (Ω), by choosing C = ∇ϕ, we can obtain Thus, there is no need to add a Lagrange multiplier in the magnetic equation as in [27].

Remark 2.
The SVC scheme's boundary condition for the pressure filed is artificial Neumann boundary [15,20]. From the Equation (6)

Remark 3.
In a bid to prove the SVC scheme's stability, we need to define Theorem 1 (SVC scheme's stability). For all 0 ≤ k ≤ N, for any (υ, p, b, ϑ) that satisfy the Algorithm 1 is stable in the sense that

Remark 4.
The RVC scheme's boundary condition for the pressure filed is consistent boundary [21]. From the Equation (18) (17), pressure field boundary condition is obtained

Remark 5.
To decouple the magnetic field b and the velocity field υ, the velocity υ k can be computed directly from the Equation (15).
we get a linear equation as follows thus, we can calculate b k+1 , by the aid of b k+1 × n| ∂Ω = 0. Remark 6. Sinceυ k is known, we can directly calculate the temperature in (16) and take it into (17) to apply directly. We use a special technique to obtain p k+1 from the third equation of (17) and Equation (18). Using (17) taking divergence for the above equation with ∇ · (∇ × υ) = 0 and ∇ · υ k+1 = 0, we can get it is worth noting that we can calculateυ k−1 with the initial value and backward Euler scheme.

Remark 7.
Finally, we updateυ k+1 from the original equation with the boundary condition then, all the unknowns variables υ, b, ϑ, p are fully calculating.

Remark 8.
In a bid to prove the stability, we use the Gauge-Uzawa format. Then, we recommend a Gauge variable ξ k and an instrumental variable α k , namely Theorem 2. (RVC scheme's stability) For all 0 ≤ k ≤ N, for any (υ, p, b, ϑ) that satisfy the Algorithm 2 is stable in the sense that

Algorithm 2 RVC Scheme
Step 1. Find υ k+1 ∈ X, b k+1 ∈ W such that: Step 2. Solve ϑ k+1 ∈ Y 0 such that: Step 3. Solve υ k+1 ∈ X, p k+1 ∈ Q such that: Step 4. Update υ k+1 ∈ X by working outυ k+1 ∈ X: Proof. Taking inner product from both sides for the first equation of (15) the second equation of (15), taking inner product with 2∆tυ k , we find For the first equation of (16), taking inner product with 2∆tϑ k+1 , using Lemma 2 and the property (υ · ∇s, s) = 0, we achieve For the first equation of (17), taking inner product with 2∆tυ k+1 , we achieve Taking inner product with itself for the first equation of (14), we get with the definition of ξ k , we obtain Finally, summing up (19)- (22), (24) from k = 0 to N − 1 and getting rid of some positive terms, there holds thus, the proof is ended.

Numerical Experiments
In this section, we provide some numerical examples to validate the efficiency and accuracy of the RVC scheme for the time-dependent thermally coupled MHD problem. The first is to prove the optimal convergence performance of the scheme. The second is a thermal driven cavity flow problem without an exact solution. The last is sinusoidal hot cylinder problem. We use uniform finite elements P b 1 to discrete υ and ϑ, the b is implemented by H(curl)-conforming Nédélec element, the continuous P 1 finite element for discretizing p.

Exact Solution with a Smooth Solution
In the first experiment, we take into account exact solution problem to test the convergence results. A smooth solution is presented in Ω = [0, 1] × [0, 1], we assume the following functions , b = (a sin(πx 1 ) cos(πx 2 ) cos(t), −a sin(πx 2 ) cos(πx 1 ) cos(t)), We set as a = R e = R m = S = κ = 1, β = (0, 1), the time size ∆t = h 2 and the decreasing mesh sizes 1 50 , 1 60 , 1 70 , 1 80 . The numerical error results of each variable is shown in Table 1. We can get that each variable reaches the optimal error; namely, the error order of the 2 ) which is consistent with reference [15]. Since one of the advantages of the RVC scheme overcomes the difficulty caused by the artificial pressure Neumann boundary condition. Then we mainly changed the equation coefficients, namely R e , R m , κ and a to verify that the change of the coefficients will reduce the L 2 convergence order of the υ, and the convergence order of other variables remains optimal, see Tables 2-5.
In addition, we show the pressure error field at T = 1, the time size ∆t = h 2 , h = 1 64 for a typical time step using the SVC and the RVC schemes refer to [15]. Figure 2a produces a numerical boundary layer, there is no numerical boundary layer in Figure 2b, but we observe large spikes at the four corners of the domain. This test suggests that the divergence correction of the rotational form successfully cured the numerical boundary layer problem. However, the large errors at the four corners degrade the global convergence rate of the pressure approximation. Finally, we also do the basic agreement between the exact solutions and numerical solutions of the variables at different grid in Figure 3. From this, we find that the coarseness of the grid has little effect on the error.
In Figures 5-8, we compare the results of the magnetic, temperature, velocity and pressure fields for T = 0.1, 0.5, 1 when the coefficients are the same. With the change of time, the magnetic field streamlines gradually bend, and the vortex of the velocity streamlines gradually move to the right. There will be three more vortex in the square cavity except the lower left corner, and the temperature and pressure will also produce obvious changes at T = 1. In Figures 9-12, we compare the effect of R e with β = (0, 1), S = R m = κ = 1. With the continuous increase of R e , the magnetic field streamlines have obvious bend, at the same time, the velocity vortex moves to the left and a small vortex is generated at the corner of the area when R e = 1000. The isotherm is bent because the high temperature liquid flows to the low temperature and the low temperature liquid flows to the high temperature. In Figure 12, the pressure shows a gradual decrease.

Sinusoidal Hot Cylinder Problem
Last, we test an important practical problem, called sinusoidal hot cylinder problem, which has been investigated in [28]. We define the test domain Ω = [0, 1] × [0, 1], the source terms f = g = 0, Ψ = 0, the initial values υ 0 = b 0 = 0, ϑ 0 = 0, the time sizes ∆t = 0.000625, the mesh size 1 40 in Figure 4b-d and the boundary conditions are given as follows: where b D = (0, 1). Firstly, we set the radius of inner cylinder (r in = 0.1), R e = S = R m = κ = 1 for different β = (0, 10), (0, 100), (0, 1000). In Figures 13-16, we find that the magnetic field lines are slightly curved. The high temperature circle gradually shifts upward in the center of the Figure 14. Two main vortices are captured in velocity filed, along with β increases the velocity streamlines follow the shape of a cylinder, the two main vortices, it is obvious that the vortex in the center has changed from slender to stubby. In Figure 16, the pressure difference is decreasing as the fluid flows.     Then, we set β = (0, 1000) and R e = S = R m = κ = 1 for different radius of inner cylinder (r in = 0.1, 0.2, 0.3) in Figures 17-20. As the radius of the inner cylinder gradually increases, the streamlines of the magnetic field gradually flatten from bending, the fluid does not have enough space for rotation, so the conduction mode dominates. Besides, four vortices are captured. The isotherm follows the shape of the cylinder, and the upward movement of the two main vortices gradually disappears. The pressure changes obviously with the increasing of the radius of the cylinder.
The above phenomena found in Figures 13-20 are very consistent with [28][29][30]. Therefore, Our method can simulate this model well.

Conclusions
In this paper, we accomplished the following goals: (i) give a linear full decoupling velocity correction schemes for thermally coupled incompressible MHD system; (ii) the stability of the proposed algorithms are derived; (iii) the numerical results verified the reliability and accuracy of the proposed RVC algorithm. The above numerical experiments show that the RVC method is a powerful tool for solving the problem and can deal with complex problem. In the later work, we will go on study about full decoupling numerical method for settling unsteady thermally coupled incompressible MHD equations.  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.