A Mixed Finite Element Method for Stationary Magneto-Heat Coupling System with Variable Coefficients

In this article, a mixed finite element method for thermally coupled, stationary incompressible MHD problems with physical parameters dependent on temperature in the Lipschitz domain is considered. Due to the variable coefficients of the MHD model, the nonlinearity of the system is increased. A stationary discrete scheme based on the coefficients dependent temperature is proposed, in which the magnetic equation is approximated by Nédélec edge elements, and the thermal and Navier–Stokes equations are approximated by the mixed finite elements. We rigorously establish the optimal error estimates for velocity, pressure, temperature, magnetic induction and Lagrange multiplier with the hypothesis of a low regularity for the exact solution. Finally, a numerical experiment is provided to illustrate the performance and convergence rates of our numerical scheme.


Introduction
In this paper, we develop a steady-state magnetohydrodynamics (MHD) system coupling thermal problem with the physical parameters dependent on temperature in R 3 [1][2][3], as follows: − div [κ(θ)∇θ] + u · ∇θ = ψ in Ω, where n is the outer unit normal of ∂Ω, and Ω is a bounded domain with a Lipschitz boundary ∂Ω; (u, p, B, r, θ) denote the velocity field, pressure, magnetic induction, scalar fuction and temperature; (ν, σ, µ, κ, β) denote the kinematic viscosity, electric conductivity, coupling number, thermal conductivity and thermal expansion coefficient; (ψ, f , g) denote a given heat source, a forcing term for magnetic induction and the known applied current with div g = 0. In fact, the scalar function r is a virtual function. The purpose of adding r is related to the constraint div B = 0 [4,5]. Many numerical methods for incompressible MHD problems have been widely studied. According to our investigation, a considerable number of scholars in the literature use continuous Lagrange finite elements to approximate velocity and magnetic induction unknowns, see, e.g., [6][7][8][9][10]. However, it is well known that using continuous elements to approximate magnetic induction may lead to an inaccurate approximation when the regularity of the magnetic unknown is lower than H 1 (Ω), cf. [11,12], which may often be encountered in non-convex polyhedral or with a non-C 1,1 boundary. A novel method to overcome these shortcomings was put forward in [13] by using Nédélec finite elements to approximate the magnetic unknown B. This seems attractive and has been used in [14][15][16][17] and the references therein. In addition, the authors had proposed different formulas to keep the magnetic solutions divergence-free in the numerical schemes in the papers [18,19].
Because of the movement of the fluid with viscosity, viscosity will produce heat. Under the action of the external magnetic field, the incompressible MHD problem is usually coupled with the thermal system through the famous Boussinesq approximation. For example, Meir et al. [20,21] studied the mixed finite element method for the thermally coupled MHD equation by adopting continuous elements to approximate the unknowns of magnetic, fluid and thermal problems, which is pioneering work. Coupling fluid systems or electromagnetic models with coefficients dependent on temperature face the mathematical challenge of solving strongly nonlinear partial differential equations, as studied by some physicists and mathematicians, see, e.g., [22][23][24][25][26][27][28][29] and the references therein. In addition, in many practical applications of the MHD system, the change in temperature will lead to a change in the coefficients in the fluid field and electromagnetic field. For the coupling MHD problem whose coefficients depend on temperature, it is important to study its reliable finite element method.
In this paper, we aim to give a rigorous well-posed analysis of the solution to the continuous problem and error estimates for the MHD system with variable coefficients by the finite element method. The fluid field was approximated by Taylor-Hood-type finite elements, and the thermal system was approximated by Lagrange finite elements. Considering the highly nonlinearity brought by the physical parameters dependent on temperature and the Lorentz terms in the magnetic equation, as well as a possible nonconvex domain or a non-C 1,1 boundary, we choose H(curl)-conforming Nédélec edge elements to approximate the magnetic equation to capture the physical solutions. The optimal error estimates of velocity, pressure, temperature, magnetic induction and Lagrange multiplier are established. As far as we know, there still lacks rigorous analysis in the literature for the error estimate of the stationary MHD thermally coupled model with temperature-dependent coefficients.
The remaining parts of this paper are organized as follows. In Section 2, we introduce some notations and basic finite element estimation used in the discussion and show the uniqueness of the solution to the continuous system. In Section 3, we propose a discrete finite element method for the variable coefficient system consisting of Equations (1)- (6). In Section 4, we give the convergence of all variables under the slightly smooth regularity assumption. In Section 5, we verify the effectiveness of the proposed numerical method through a numerical experiment.

Notations for the Variable Coefficients Model
Firstly, we introduce some symbols that will be used throughout this article. For all m ∈ N + , 1 ≤ p ≤ ∞, let W m,p (Ω) denote the standard Sobolev space, and when p = 2, it can be written as H m (Ω). The notation (·, ·) is expressed as an inner product, namely (φ, ψ) = Ω φψ dx, and the norm in L 2 (Ω) defined by · 0 . Vector-valued quantities will be denoted in boldface notations, such as u = (u 1 , u 2 , u 3 ) and L 2 (Ω) = (L 2 (Ω)) 3 . We use C and c to denote generic positive constants independent of the mesh size h, and it may adopt different values in different places.
To simplify, we define the following Sobolev spaces The norm of the following types still need to be defined: We then set H(div; Ω) = {b ∈ L 2 (Ω), div b ∈ L 2 (Ω)}, and H(Ω) = W 0 ∩ H(div; Ω), which is equipped with the following norm: To facilitate our analysis, the following embedding results (see, e.g., Proposition 3.7 of [30] or [31]) are introduced here, which are also valid for the Lipschitz polyhedron domain.
In order to better demonstrate the stability of energy, the following trilinear terms are denoted Next, the definition of the weak solution to the magneto-heat coupling system with variable coefficients (Equations (1)-(6)) is given.

Remark 2.
The existence of a solution to the system in Equations (7)- (9) can be referred to the work of [34] (Section 4), albeit under additional smoothness assumptions on the domain. Here, we restrict ourselves to proving the following (more straightforward) uniqueness result.

Uniqueness of Continuous Problems
Throughout the paper, we set that σ(θ), κ(θ), β(θ) and ν(θ) are Lipschitz continuous and satisfy 0 Before proving the uniqueness of continuous problems, we first introduce some basic knowledge.

Lemma 2.
For any u, B and θ that satisfy Equations (7)-(9), the following estimates hold For the convenience of the subsequent analysis, setting The norms are defined as: It is well known that (cf. [31,35]) both H 1 0 (Ω) × Q and W × H 1 0 (Ω) satisfy the corresponding inf-sup conditions, namely, and where the generic constant C only depends on Ω.
Before we begin the proof, the following assumptions should be given: where C lip is the Lipschitz constant. Now, we are going to prove that the solution of the continuous problem is unique.

Finite Element Analysis for Magneto-Heat Coupling System
In this section, we introduce a mixed finite element approximation of the MHD system coupled thermal problem in Equations (7)- (9). The approximation is based on the Nédélec first family of elements for the discretization of the magnetic induction.
Throughout, the domain Ω is partitioned into a finite number of open non-overlapping subdomains with regular and quasi-uniform meshes T h of mesh-size h that partition Ω into tetrahedra K. Each tetrahedron K is supposed to be the image of a reference tetrahedronK under an affine map F K . P k (K) and P k (K) represent the space of polynomials of the total degree at most k ≥ 0 and homogeneous polynomials k on K, respectively.
Given the generalized Taylor-Hood element (X k h , Q k−1 h ) with k ≥ 2, where X k h is the k order vectorial Lagrange finite element subspace of X, and Q k−1 h is the k − 1 order scalar Lagrange finite element subspace of Q. For k = 1, the velocity and pressure pair can be approximated by the well-known stable mini-elements, cf. [31,35,36]. Furthermore, Y k h is the k order scalar Lagrange finite element subspace of Y, refer to [31,35].
For D k (K), it denotes the polynomials q in P k (K) that satisfy q(x) · x = 0 on K. Define the following space where 1 ≤ k. Using Nédélec H(curl)-conforming finite element space (see [33,37]) and we can define the following weakly divergent space In addition, the following discrete Poincaré-Friedrichs inequality is established with a constant C p > 0 independent of the mesh size h. The link between the spaces W k 0h and W (Ω) is accomplished by the Hodge mapping Z : H(curl; Ω) → W (Ω)-refer to [12]-where W (Ω) = {C ∈ H(Ω), div C = 0 in Ω} such that curl Z (C) = curl C ∀ C ∈ W.
Furthermore, there exists l = l(Ω) > 0, In addition, the discrete kernel space of the divergence operator is given by On a quasi-uniform mesh, there holds (see Theorem 3.2.6 of [38]) where C inv > 0 is a generic constant independent of the mesh size h, ı and m are two real numbers with 0 ≤ ı ≤ m ≤ 1, and p and q are two integers with 1 ≤ p ≤ q ≤ ∞.
The following discrete inf-sup conditions (see Chapter 2 of [36] or [12]) are founded where β * is a generic positive constant depending on the domain Ω.
Remark 3. The temperature-dependent coefficients will greatly enhance the nonlinearity of the problem and make the analysis more complicated. The existence of a solution to Equations (25)- (28) can be displayed by Brouwer's fixed point theorem-for details, refer to Section 4.3 of [39].

Remark 4.
For all s h ∈ S k h , by selecting C = ∇s h in Equation (27), it can be derived −(∇r h , ∇s h ) = (g, ∇s h ). Thus, for a solenoidal source term g, it is natural to deduce r h = 0.
In order to estimate the error in the next section, the stability of the numerical scheme in Equations (25)-(28) should be given here.
Proof. The proof is parallel to that of Theorem 1.

Convergence Analysis of the Magneto-Heat Coupling Problem
In this section, the convergence of the MHD system coupled the heat equation with variable coefficients is considered by employing the finite element method. We strictly establish optimal error estimates of velocity, pressure, temperature, magnetic induction and Lagrange multiplier under the assumption that the exact solution has low regularity.
Here, it is necessary to make the following regularity assumptions for the weak solution of Equations (7)-(9), which will facilitate the error estimate of the discrete solution. Assumption 1. Assum that the solution (u, p, B, θ, r) satisfies the following regularity: where the exponent s > 1/2 depends on Ω.
For the convenience of the subsequent analysis, we will assume there exists a constant C f depending on f , g, ψ and Ω such that u 1+s,2 + p s,2 + curl B s,2 + θ 1+s,2 + ∇r s,2 ≤ C f .
Let = min{k, s}, for k ≥ 1 and s > 1/2, with the help of [31,40], then we have the following approximation properties: where k and s are the order index of the finite element spaces and the regularity of the exact solution, respectively.
Let (e u , e p , e θ , e B , e r ) = ( . A combination of Equations (7)- (9) and Equations (25)-(28) yields the following truncation error equations: With the preparations of the above work, we now begin to study the optimal error estimation of each variable.
is satisfied. Then, the weak formulation made up of Equations (7)-(9) and the discretization scheme in Equations (25)-(28) has a unique solution (u, p, θ, Proof. As a first step, the test functions of the momentum and magnetic equations are constrained in the discrete kernel space. Let (v, C) ∈ X k 0h × W k 0h . Using the orthogonality property, with Equation (32), we have According to Lemma 1 and the property of Hodge mapping (Equation (21)), by setting 1/(3 + δ 1 ) + 1/(6 − δ 2 ) = 1/2, it is easy to see that and similarly As a consequence of the previous calculation, we can estimate Equation (36) as follows Since h . The right-hand side of Equation (36) has the following estimate Using Equation (33), we can arrive at The left-hand side of Equation (39) has the following estimate The right-hand side of Equation (39) has the following estimate Using Equation (34), we derive The left-hand side of Equation (42) has the following estimate Since C − B h belongs to the kernel space W k 0h , we have (∇e r , C − B h ) = (∇(r − j), C − B h ) for any j ∈ S k h . The right-hand side of Equation (42) has the following estimate Combining Equations (37), (40), (43), (38), (41) and (44), together with (32)- (34), it can be checked that By the triangle inequality, there holds for all (v, C) belongs to the discrete kernel space (X k 0h , W k 0h ), ϕ ∈ Y k h , q ∈ Q k−1 h and j ∈ S k h . In the next step, let (v, for any (q, j) ∈ Q k−1 h × S k h . From the inf-sup condition and the continuity of b(·, ·) and a(·, ·), there exists a solution to this problem that satisfies See [36] for more details. Then, (w + v, d + B) ∈ X k 0h × W k 0h can be inserted into Equation (45). With the triangle inequality, we deduce Applying Equations (46) and (47) shows that the proof is completed.
In addition, by using the inf-sup condition (Equation (23)), we obtain the following estimates for pressure and the Lagrange multiplier. Theorem 4. Suppose that Assumption 1 and Theorem 3 are satisfied. The continuous system of Equations (7)-(9) and the numerical scheme in Equations (25)- (28) has the unique solution (p, r) and (p h , r h ), respectively, which satisfies Proof. From Equations (32) and (34), we know Combining Assumption 1 and Theorem 3, we can conclude Theorem 4.
With the help of (31), the following theorem draws the conclusion of this paper.

Numerical Experiment
In this section, we consider a numerical experiment to test the convergence rate of the numerical scheme proposed in Section 3. The parallel code is developed based on the finite element package Parallel Hierarchical Grids (PHG), cf. [41,42]. The computations are carried out on the LSSC-IV Cluster of the State Key Laboratory of Scientific and Engineering Computing, Chinese Academy of Sciences. The domain is Ω = (0, 1) 3 , and the finite element mesh is obtained by a uniform tetrahedral partition. Let T i , i = 0, 1, 2, 3 be four successively refined meshes listed in Table 1.
From Tables 2 and 3, we find that the convergence rates for u h , p h , B h and θ h are given by Here, the velocity u and the temperature θ are discretized by using the continuous P 2 finite elements, the pressure p is discretized by using the continuous P 1 finite elements and the magnetic induction B is discretized by using the first-order edge elements method. This means that optimal convergence rates are obtained for all variables.