CONVERGENCE ANALYSIS OF A FULLY DISCRETE FINITE ELEMENT METHOD FOR THERMALLY COUPLED INCOMPRESSIBLE MHD PROBLEMS WITH TEMPERATURE-DEPENDENT COEFFICIENTS

. In this paper, we study a fully discrete finite element scheme of thermally coupled in-compressible magnetohydrodynamic with temperature-dependent coefficients in Lipschitz domain. The variable coefficients in the MHD system and possible nonconvex domain may cause nonsmooth solutions. We propose a fully discrete Euler semi-implicit scheme with the magnetic equation approximated by N´ed´elec edge elements to capture the physical solutions. The fully discrete scheme only needs to solve one linear system at each time step and is unconditionally stable. Utilizing the stability of the numerical scheme and the compactness method, the existence of weak solution to the thermally coupled MHD model in three dimensions is established. Furthermore, the uniqueness of weak solution and the convergence of the proposed numerical method are also rigorously derived. Under the hypothesis of a low regularity for the exact solution, we rigorously establish the error estimates for the velocity, temperature and magnetic induction unconditionally in the sense that the time step is independent of the spacial mesh size.


Introduction
Magnetohydrodynamic (MHD) is the theory of macroscopic interaction of conductive fluid and electromagnetic induction.It consists of a viscous, incompressible fluid which has the property of electric current conduction and interacting with electromagnetic inductions.There are lots of applications in astronomy and geophysics as well as engineering problems, such as metallurgical engineering, electromagnetic pumping, stirring of liquid metals, liquid metal cooling of nuclear reactors, refer to [16,21,36,42], and the references therein.However, in many cases, the effect of temperature can not be ignored.Especially, the change of temperature will cause the change of fluid coefficients [52,54] as well as the magnetic field coefficients [10,34].
In this work, we consider the following transient incompressible Navier-Stokes equations and Maxwell's equations coupled to the heat equation with temperature-dependent coefficients in R 3 [3,10,16,42] as follows, div  = 0 in   , (1.4) where   = Ω × (0,  ),  > 0 is a given finite final time, Ω is a bounded, simply-connected and Lipschitz polyhedral domain. denotes the velocity field,  the pressure,  the magnetic induction,  the temperature,  the kinematic viscosity,  the electric conductivity,  the magnetic permeability,  the thermal conductivity,  a given heat source,  the thermal expansion coefficient,  a forcing term for the magnetic induction,  the known applied current with div  = 0.The system is considered in conjunction with the following initial values and boundary conditions, (, 0) =  0 , (, 0) =  0 , (, 0) =  0 ∀ ∈ Ω, (1.6) where   = Ω × (0,  ),  is the outer unit normal of Ω and the initial magnetic induction  0 satisfies div  0 = 0. Our results in this paper are also valid for another frequently used set of boundary conditions of the magnetic induction  for (1.1)- (1.5) given by  •  = 0,  × curl  = 0 on   .
In the last several decades, various finite element methods for MHD problems regardless of heat effects have been extensively developed in the literature.Let us review the references and try to summarize them, which is unlikely to be complete and accurate, of course.We can mainly classify them into two formulations based on the chosen discrete finite element spaces.The first formulation was proposed by Gunzburger et al. [26] and they made use of standard Lagrange  1 finite element spaces to approximate both the hydrodynamic unknowns and the magnetic induction.It is therefore easy to implement and has been extensively employed in computational MHD, see e.g.[20,25,35,55] for stationary models and [4,17,27,58] for time-dependent models.However, it is known that the nodal finite element method discretizations of the magnetic operator cannot be correctly approximated when the magnetic induction components may have regularity below  1 (Ω), cf.[14,30], which may be frequently encountered in non-convex polyhedral or with a non  1,1 boundary.A possible way to overcome these difficulties was proposed in [48] by virtue of Nédélec finite elements for the magnetic induction , which leads to another natural formulation and is valid for non-smooth magnetic solution.This variational formulation seems to be attractive and has been employed in [19,24,45,50] and the references therein.We also mention the recent papers [31,32], where different formulations were proposed to maintain the divergence free magnetic solutions in the numerical schemes.
When buoyancy effects cannot be neglected in the momentum equation due to temperature differences in the conductive flow, the incompressible MHD is usually coupled to the heat equation by the well-known Boussinesq approximation.For instance, as a pioneer work, Lagrange continuous finite element methods for the stationary heat coupled MHD equations with constant coefficients had been studied in [39,40].On the other hand, in many practical applications of MHD problems, the change of temperature will affect the coefficients in the fluid field as well as the electromagnetic filed.It is of great significance to study reliable finite element methods for the coupled MHD system with the coefficients depend on the temperature.The coupled fluid systems or electromagnetic models with temperature-dependent coefficients are faced with nonlinear PDEs with great mathematical challenges, which has attracted attention by both physicists and mathematicians, see e.g.[10,33,43,54] and the references therein.As far as we know, the first work to study error estimates of finite element methods for the thermally coupled MHD equations with temperature-dependent coefficients is given in [47], where a fully discrete Crank-Nicolson scheme is proposed and investigated.More recently, the fully discrete Euler semi-implicit scheme has been studied in [46].The proposed schemes in the above two papers are based on the magnetic induction approximated by Lagrange  1 finite element method and all the error estimates are conducted under sufficiently smooth assumption on the exact solutions.However, in view of the highly nonlinearity brought by the temperature-dependent coefficients and the Lorentz terms in the magnetic equation, as well as a possible non-convex domain or a non  1,1 boundary, we can not expect to a smooth solution for the magnetic induction belong to  1 (Ω) in such situations (see [12,13]).Thus (curl)-conforming Nédélec edge element is a natural choice to approximate the magnetic equation in order to capture the physical solutions.Furthermore, Nédélec finite elements seems to be a better choice since it can treat the boundary conditions of magnetic induction easier than Lagrange finite element discretization.
In this work, we will give a rigorous convergence analysis and error estimates of a fully discrete finite element method for the MHD system described by (1.1)-(1.7)based on the magnetic induction approximated by (curl)-conforming Nédélec edge element.The time discretization is based on a backward Euler semi-implicit scheme and the stable Taylor-Hood type finite elements are used to approximate the fluid field and Lagrange finite element to approximate the heat equation.In the first half of the paper, we show that the numerical solution converges to a weak solution of the continuous system without any further regularity assumption as both meshwidth and timestep tend to zero.In fact, the first convergence result of finite element discretization for MHD is attributed to [45], where some weak and strong convergence of the subsequence of discrete velocity and magnetic field are proved by the standard compact argument.In this paper, we will extend the results to MHD models with temperature-dependent coefficients.Strong convergence of all the discrete fields (velocity, magnetic induction and heat) is proved rigorously.So the result of this paper is also an improvement of [31], which has not proved any strong convergence of the discrete fields.Furthermore, we show the uniqueness of weak solution for incompressible MHD models with temperature-dependent coefficients provided it satisfies a smoother condition, which seems to be new in the literature.Then we can prove that the whole sequence of the discrete solution converges (strongly) to the unique weak solution.In the second half of the paper, under a weak regularity hypothesis on the exact solution, a rigorous error estimate of the velocity, temperature and magnetic induction are established unconditionally in the sense that the timestep is independent on the spacial meshwidth.In this paper, we will confront with two major difficulties in the analysis.The first difficulty arises from the nonlinearity in the model caused by variable temperature-dependent coefficients, and the other is the very low regularity of the exact solution caused by the nonlinearity structure or non-convex domain.In order to deal with these difficult problems, some technical tools need to be developed in this paper.
A brief overview of this paper is provided as follows.In the next section, we describe the notations and some preliminary knowledge to be used throughout the paper.In Section 3, we propose a fully finite element discrete method with a backward Euler semi-implicit scheme for the system (1.1)-(1.7),and some basic lemmas and theorems are recalled.In Section 4, the well-posedness and stability of the fully discrete scheme are presented.We show that the fully discrete solution converges to a weak solution of the continuous problem as ∆ and ℎ tend to zero.The uniqueness of the continuous problem is also established under a slight smoother assumption on the weak solution.In Section 5, we prove error estimates for all variables under a weak regularity assumption.In Section 6, we present two numerical examples to illustrate our theoretical results.Finally, we close the paper with some concluding remarks in Section 7.

Functional setting for the magneto-heat coupling model
For mathematical setting of problem (1.1)-(1.5)with the initial values and boundary conditions (1.6)-(1.7),we first introduce some notations that will be used throughout the paper.For all  ∈ N + , 1 ≤  ≤ ∞, let  , (Ω) denote the standard Sobolev space and it is written as   (Ω) when  = 2.The norm in  , (Ω) is defined by ‖•‖ , such that where For the function spaces   (0,  ; ), 1 ≤  ≤ ∞, the norms are denoted as where  is a real Banach space with the norm ‖•‖  .The inner product will be denoted by (•, •), that is (, ) = ∫︀ Ω  d, the norm in  2 (Ω) defined by ‖•‖ 0 .⟨•, •⟩ is for the dual product between a Banach space and the dual space.Vector-valued quantities will be denoted in boldface notations, such as  = ( 1 ,  2 ,  3 ) and We use  and , with or without subscripts, bars, tildes or hats, to denote generic positive constants independent of the discretization parameters, which may take different values at different places.
We introduce the following classical Sobolev spaces: Here and what in follows, we define the following norm We also need to define the following Sobolev spaces and ℋ(Ω) =  0 ∩ (div; Ω), which is equipped with the following norm We recall the following embedding result (see e.g., Prop.3.7 of [2] or [23]) which is valid for a Lipschitz polyhedron.
It is well known that the following Poincaré type and embedding inequalities are valid in bounded polyhedral domains (see Chapter 3 of [41] for more details), For every  ∈  1 2 (Ω), let  0 () be an affine space of  defined by where   ∈  is an extension of .Similarly, let  =  0 , for every  ∈  1 2 (Ω), let  () be an affine space of  defined by where   ∈  is an extension of .
We denote by C  (Ω) the space of functions  times continuously differentiable in Ω.The space C  (︀ Ω )︀ consists of functions in C  (Ω) bounded, and uniformly continuous in Ω with derivatives up to the  th order.The space C ,1 (︀ Ω )︀ consists of functions in C  (︀ Ω )︀ that are Lipschitz continuous in Ω with derivatives up to the th order, refer to [23].

A fully discrete finite element method based on Euler scheme
In this section, we introduce a mixed finite element method which describes a spatial discretization of the problem (2.5)-(2.7)based on the backward Euler scheme.
Let Ω be a polyhedral domain, and the domain is partitioned into a mesh  ℎ , where ℎ is the diameter of the element.Each tetrahedron  is supposed to be the image of a reference tetrahedron K under an affine map   .The family of meshes { ℎ } ℎ>0 is assumed to be regular and quasi-uniform.Let   () be the space of polynomials of total degree at most  ≥ 0 on  and ̃︀   () the space of homogeneous polynomials  on .We first introduce the generalized Taylor-Hood element where   ℎ is the  order vectorial Lagrange finite element subspace of ,  −1 ℎ is the  − 1 order scalar Lagrange finite element subspace of .And   ℎ is the  order scalar Lagrange finite element subspace of  , refer to [6,23].For the case  = 1, we use the well-known stable mini-elements to approximate velocity and pressure, cf.[6,8,23].
The space   () denotes the polynomials  in ̃︀   () that satisfy () •  = 0 on .For 1 ≤ , we define the space To approximate the magnetic induction, we use Nédélec (curl)-conforming finite element space (see [41,44]), which is defined by Note that the above definition is the first family (curl)-conforming discrete space.We can also employ the second family finite element spaces with   () is chosen by   (). Setting , we introduce the discretely solenoidal function space Furthermore, the space   0ℎ is known to satisfy the following discrete Poincaré-Friedrichs inequality (see [48] or Thm.4.7 of [30]), with a constant  * > 0 independent of the mesh-size ℎ.
The link between the spaces   0ℎ and (Ω) is accomplished by the Hodge mapping  : (curl; Ω) → (Ω), refer to [30], where Furthermore, the Hodge mapping satisfies the following approximation property, there exists  = (Ω) > 0, We use  −1 ℎ to approximate the pressure , and use a finite element affine space where Π ℎ =  3 ℎ , to find the velocity .Further, the discrete kernel space of the divergence operator can be defined by We recall the following inverse estimate from Theorem 3.2.6 of [11].On a quasi-uniform mesh there holds where   > 0 is a generic constant independent of the mesh-size ℎ,  and  are two real numbers with 0 ≤  ≤  ≤ 1,  and  are two integers with 1 ≤  ≤  ≤ ∞.
From the Fortin criterion, the following discrete inf-sup condition (see Chap. 2 of [8] or [30]) is established, where  * is a generic positive constant depending on the domain Ω.
Remark 3.1.The coefficients depend on the temperature will increase the nonlinearity of the model and make the problem more intricate.The existence of solution to problem (3.5)-(3.8)will be proved in next section.Furthermore, we will prove that the discrete solution converges to a weak solution of the continuous model as ∆ and ℎ tend to zero.
Remark 3.2.When Ω is a non-convex polyhedron, the difficulty comes from the fact that the magnetic induction is in general not in  1 (Ω) (see Rem. 3.3.1 of [21] for an explanation).The approximation of the magnetic induction with classical  1 -conforming Lagrange finite elements can not capture the singularities and may converge to a wrong solution, refer to [14,15,21], etc.This is the reason why we choose Nédélec edge element to approximate the magnetic induction in (3.7).
We now show that the solution of (3.5)-(3.8)enjoys a stability property, which holds regardless of the sizes of ℎ and ∆, and will play an important role in proving the convergence and well-posedness of the fully discrete solution.
The following estimates for the trilinear term in Navier-Stokes equations can be referred to Lemma 4.5 of [53].

Well-posedness and convergence of the fully discrete solution
In this section, we will prove the well-posedness of the numerical solution to the problem (3.5)-(3.8)by using the Lax-Milgram theorem (see [22] for more details).Utilizing the stability of the numerical scheme and the compactness method, the existence of weak solution to the thermally coupled MHD model in three dimensions is established.Furthermore, the uniqueness of weak solution and the convergence of the proposed numerical method are also derived.
We first prove the well-posedness result for the discrete solution in the following theorem.Proof.We divide this proof into two steps: Let us define It is easy to verify that (  ℎ ,  ℎ ) satisfies the ellipticity and boundedness.An application of the Lax-Milgram theorem shows that problem (3.8) attains a unique solution   ℎ .
Remark 4.3.A priori stability estimate of the above time derivatives for discrete finite element solution plays a key role in the subsequent strong convergence of the scheme.Such type of estimate for MHD model was first developed in [45].The proof therein seems to be not complete and there are some minor gaps in the bounds on the important Lorentz force terms (e.g., the estimates in line 19 and line 20 of page 1073 are not valid and similar problem arises in line 2-3 of page 1074).Here we give a rigorous proof with a different index .
We also need to recall the Aubin-Lions' compactness result for Bochner spaces (refer to Lem. 2.8 of [21]).
Lemma 4.5.Let  be a Banach space,  0 and  1 be two reflexive Banach spaces.Assume  0  with compact injection,  ⊂  1 with continuous injection.Then the space Next, we present some basic convergence results for the fully discrete solution in the following theorem.
Theorem 4.6.There exist functions  ∈  ∞ (︀ 0,  ;  2 (Ω) We say that the above three subsequences (still denoted by the same notations) enjoy the same accumulation function (, , ).In fact, applying the interpolation inequality, Hölder inequality and Lemma 3.5, we have with  = 6− 5 , and we continue to obtain Considering that the viscosity coefficients (•), (•), (•) and (•) are assumed to be in C 0,1 (︀ Ω × R; R )︀ , here we may define | • | C 0,1 ( Ω×R;R) by where  can be taken as , ,  and .Now, we will prove that the accumulation function (, , ) is indeed a weak solution to (2.5)-(2.7),which provides a numerical version of the existence analysis of the thermally coupled incompressible MHD problems with temperature-dependent coefficients.
We will also show the continuous system (2.5)-(2.7)has a unique solution.To this end, let us recall the Gronwall lemma in differential form.where () and  0 () are non-negative integrable functions, then there holds We will prove the uniqueness of the continuous system (2.5)-(2.7),provided that the exact solution is under a slight smooth assumption.More precisely, we need to make a smoother assumption on the weak solution for the magneto-thermal coupling model with temperature-dependent coefficients.

Error estimates for the magneto-heat coupling model with temperature-dependent coefficients
In this section, we mainly consider the error estimates of the fully discrete finite element method for the MHD system coupled the thermal equation with temperature-dependent coefficients.Under the hypothesis of a low regularity for the exact solution, we rigorously establish the error estimates for the velocity, temperature and magnetic induction unconditionally in the sense that the time step is restricted but is independent of the spacial mesh size.We also prove a sub-optimal error estimate for the pressure as a supplementary result, which is consistent with the pioneering work [54].
We first recall a discrete version of the Gronwall inequality in a slightly more general form than usually found in the literature, and this detailed proof can be found in [29].
Lemma 5.1.Let  * , ∆,   ,   ,   and   be non-negative numbers with  ≥ 0 such that Suppose that ∆  < 1, for all , and set   = (1 − ∆  ) −1 .Then Before proceeding, we need to make a regularity assumption for the weak solution of (2.5)-(2.7),which will be helpful for the error analysis of the numerical solution.
Remark 5.3.The above regularity assumption for the solution (, , , ) is even weaker than most of the hypotheses in the literature, see, e.g.[27,29,38,49,54].We hope it is a reasonable assumption and may be valid for a general polyhedral with Lipschitz boundary.In fact, it is enough to facilitate the subsequent error analysis of the magneto-thermal coupling model with temperature-dependent coefficients.
By virtue of the properties of projection, we present that the solution to (2.5)-(2.7)has the following estimates, which hold regardless of the sizes ℎ and ∆.Moreover, the stability property is beneficial to the subsequent error analysis in full discretization.
Lemma 5.4.Suppose Assumption 5.2 holds.Let (, , ) be the unique solution of (2.5)-(2.7),then the following estimates are established, where   is a generic constant depending on the regularity of the domain Ω.
Proof.Let  =  ℎ  ∈   ℎ , with  is the standard Lagrange nodal interpolant.By using the finite element approximations, including (5.3), inverse inequality (3.3), Sobolev's embedding theorem and Assumption 5.2, we directly see We can verify that other terms are bounded in a similar way, this bring the proof to an end.
With the above preparations, we can now establish the following error estimates for the velocity, temperature and magnetic induction.
As far as error estimate for the pressure, we will prove the following sub-optimal result by following a similar argument developed in [54].
Proof.From (5.5), we know that Making use of the inf-sup condition (3.4), we derive where    1ℎ =   1ℎ +   1ℎ +   ℎ +   1ℎ .We will estimate the terms in the last line of (5.27) one by one.According to the estimate (5.25), there holds Concerning the term ⟨  1ℎ ,  ℎ ⟩, by (3.15), similar to the previous theorem, we can deduce A combination of the inverse inequality, Poincaré type inequality and (5.25) yields Hence it holds that In a similar way, we can prove )︁ .
To bound the last term in the last line, we continue to derive By virtue of (5.)︁ .

Numerical experiments
In this section, we consider two numerical experiments to verify the convergence properties of the fully finite element discretization for the MHD coupled thermal equation with temperature-dependent coefficients (1.1)-(1.7).The parallel code is developed based on the finite element package-Parallel Hierarchical Grids (PHG), cf.[56,57].The computations were (partly) done on the high performance computers of State Key Laboratory of Scientific and Engineering Computing, Chinese Academy of Sciences.
The two examples are used to verify the optimal error estimates of the fully discrete finite element method proposed in Section 3. In all examples, the domain under consideration is Ω = (0, 1) 3 and the finite element mesh is obtained by a uniform tetrahedral partition.We employ the continuous  2 finite element for discretizing the velocity  and the temperature , the continuous  1 element for discretizing the pressure , and the second order edge element method for discretizing magnetic induction .Since the exact solution is linear in spatial variables, the major error comes from the discretization of the time variable.We fix a tetrahedral mesh with ℎ 0 = 0.433, the terminal time  = 1, and test the convergence rate at each time step.We list the discretization error for all unknowns at the last moment  = 1 in Table 1, from which it shows perfectly that the temporal convergence rate of the Euler discrete scheme is first-order.
Setting ∆ 0 = ℎ 2 0 , Table 2 shows that the convergence rate for the backward Euler scheme at the terminal time.The initial conditions, boundary conditions and source terms are determined by the analytical solution.Both the timestep and the meshwidth are refined at the same time such that ∆ =  (︀ ℎ 2 )︀ .The corresponding convergent results are displayed in Table 2 and an  (︀ ℎ 2 )︀ convergence of the proposed numerical scheme can be observed, which agrees with the theoretical results developed in this paper.

Summary
We have studied a fully discrete finite element scheme for the 3D thermally coupled incompressible MHD problems with variable coefficients problems.The proposed scheme has the nice features that it only needs to solve one linear system at each time step and the magnetic induction is approximated by (curl)-conforming Nédélec edge element, which make it quite attractive to solve these highly nonlinear MHD models with possibly non-smooth magnetic induction solution.We prove that the fully discrete solution converges to a weak solution of the continuous problem as both meshwidth and timestep size tend to zero, and it is unique under a further smooth assumption.Thus we have given a numerical verification of the existence of weak solution to this model, which is still missing in the literature.Under a quite low regularity assumption for the exact solution, we rigorously establish the error estimates for the velocity, temperature and magnetic induction unconditionally in the sense that the time step is restricted but is independent of the spacial mesh size.Whether the plain convergence or error estimate, the technique and results of this paper have some improvements on that of relevant papers.

Table 1 .
6, the following error bounds are denoted by The convergence rate of Euler scheme at terminal time  = 1 (Example 6.1).

Table 2 .
The convergence rate of Euler scheme at terminal time  = 1 (Example 6.2).