A FULLY DISCRETE FINITE ELEMENT METHOD FOR A CONSTRAINED TRANSPORT MODEL OF THE INCOMPRESSIBLE MHD EQUATIONS

. In this paper, we propose and analyze a fully discrete finite element method for a constrained transport (CT) model of the incompressible magnetohydrodynamic (MHD) equations. The spatial discretization is based on mixed finite elements, where the hydrodynamic unknowns are approximated by stable finite element pairs, the magnetic field and magnetic vector potential are discretized by 𝐻 ( curl )-conforming edge element. The time marching is combining a backward Euler scheme and some subtle implicit-explicit treatments for nonlinear and coupling terms. With these treatments, the fully discrete scheme is linear in the implementation and the computation of the magnetic vector potential is decoupled from the whole coupled system. The most attractive feature of this scheme that it can yield the exactly divergence-free magnetic field and current density on the discrete level. The unique solvability and unconditional stability of the scheme are also proved rigorously. By utilizing the energy argument, error estimates for the velocity, magnetic field and magnetic vector potential are further demonstrated under the low regularity hypothesis for the exact solutions. Numerical results are provided to verify the theoretical analysis and to show the effectiveness of the proposed scheme


Introduction
The incompressible MHD system models the dynamic behaviors of an electrically conducting fluid under the influence of an external magnetic field.It couples the incompressible Navier-Stokes equations and the Maxwell's equations via the Lorentz force and the Ohm's law.There are lots of related applications in different fields, such as astrophysics, engineering related to liquid metal, controlled thermonuclear fusion.Theoretical analysis and mathematical modeling of the MHD equations can be found in [9,14] and the references therein.
Let Ω be a simply connected bounded domain with Lipschitz-continuous boundary Γ := Ω in R 3 .In this paper, we consider numerical approximation of the following incompressible MHD equations: where  > 0 is the final time,  = (0,  ],  is the density of fluid,  is the velocity of fluid,  is the pressure,  is the magnetic field,  is the magnetic induction,  is the current density, and  is the electric field.The physical parameters are the dynamic viscosity , the magnetic permeability  and the electric conductivity . In this paper, we consider the following initial and boundary conditions, (, 0) =  0 (), (, 0) =  0 () in Ω, where the initial values satisfy ∇ •  0 = ∇ •  0 = 0 and  is the unit outer normal vector on Γ.
Recently, divergence-free finite element methods have attracted more and more interest in numerical solving the MHD equations.As for the MHD model in (1), there are three divergence-free constraints to be satisfied.The first constraint is div = 0, that is, the mass conservation of fluid.From [12,25], when the numerical methods lack of point-wise impressibility, they may introduce a pressure-dependent consistency error which can potentially pollute the computed velocity.Under certain extreme situations, they can lead to spurious results and compromise the fundamental invariant property of incompressible flows.Cockburn et al. [5] proposed and analyzed a class of discontinuous Galerkin methods for the incompressible Navier-Stokes equations yielding exactly divergence-free solutions, in which the velocity is approximated by face finite element spaces.Based on this work, Greif et al. [16] proposed a mixed DG finite element method for the stationary MHD model.
The next constraint is the divergence-free condition of .It is well-known that it is a precise physical law in electromagnetism which implies that there is no monopole of magnet.From [4,8,37], the violation of the divergence-free condition on the discrete level will introduce a strong non-physical force and small perturbations to this condition can lead to large errors in numerical simulation.In the literature, the divergence-free constraint condition for  is relaxed by adding a grad-div augmented term −∇∇ •  [17,19,42] or a Lagrange multiplier  ∈  1 [13,34,36,41] to numerical formulations.By doing so, the discrete magnetic induction  ℎ is only weakly divergence-free in the best case.In [23], Hu et al. proposed a magnetic-electric based MHD model where the electric field  is kept and regarded as an independent variable.By doing so, div ℎ = 0 is ensured exactly and directly by using mixed finite element methods together with techniques from finite element exterior calculus.Later, Hiptmair et al. [22] proposed a fully divergence-free method for the incompressible MHD system by introducing a new vector potential formulation.For the MHD vector potential model, the constraint condition for  is kept exactly by the definition,  ℎ := curl  ℎ .We also refer to the review paper [37] and the references therein to other approaches.
The last constraint is div = 0, which is obtained by taking divergence of the ampere equation  = curl .Physically, this condition means that the charge is conservative.In addition, it also plays an important role in keeping the calculation accuracy for the numerical simulation [26,31,32].For the -based MHD model, the discrete current density is exactly divergence-free naturally as a consequence of  ℎ := curl  ℎ .For the inductionless MHD model, the current density  is a primitive variable and the divergence-free constraint on it can be satisfied either by postprocessing [32,33] or mixed finite element methods in (div, Ω) ×  2 (Ω) [27,28,[38][39][40].We refer to [26,29] for further arguments for the importance of the three divergence-free conditions.
Based on the discussion carried out in the previous paragraphs, it is a challenging but interesting work to develop stable numerical schemes which satisfies the three divergence-free conditions simultaneously.Recently, Li et al. [29] utilized the potential-based method and constrained transport method to propose a stable finite element method for stationary incompressible MHD equations, in which the velocity, the current density, and the magnetic induction are exactly divergence-free on the discrete level.Compared with traditional potentialbased and CT methods, the new method treats magnetic field  and magnetic vector potential  as individual variables by means of edge finite element discretization.The discrete current density  ℎ := curl  ℎ and the discrete magnetic induction  ℎ := curl  ℎ are divergence-free naturally.To fulfill the divergence-free constraint on , they approximated the velocity by combing divergence-conforming finite element spaces with DG approach.
Here and afterwards, we refer to this system simply as CT-MHD model.The motivation of this work is to extend the divergence-free finite element method to the non-stationary setting.For the sake of simplicity, we are mainly concerned with the divergence-free constraints on  and  .
The purpose of this paper is to present and analyze a fully discrete finite element method for the nonstationary CT-MHD system.The finite element approximation is based on mixed finite elements, where fluid stable elements are used for the Navier-Stokes equations, (curl)-conforming edge elements are used for the magnetic equation and magnetic vector potential equations.For the temporal discretization, we use a semiimplicit Euler scheme.At each time step, the fully discrete scheme is linear in the implementation and the magnetic vector potential is solved separately from the whole system.As a result, the proposed scheme ensures the exactly divergence-free approximation of the magnetic induction and current density.Moreover, we prove that the scheme is uni-solvent and unconditionally energy stable.Under the low regularity hypothesis for the exact solutions, we utilize the energy argument to establish the error estimates for the velocity, magnetic field and magnetic vector potential.A series of numerical tests are carried out to confirm our theoretical results and verify the performance of the scheme.
Here we would like to specify the main differences between this paper and the most relevant [22,29].On the one hand, compared to the work in [22], they only used the magnetic vector potential as the variable, then the magnetic field and current density are defined by  := curl  and  := −  +  × .Thus, the scheme therein only can yield the exactly divergence-free magnetic field on the discrete level.While our scheme can ensure the exactly divergence-free approximation of the magnetic induction and current density.In addition, there is a little different in guaranteeing the uniqueness of the magnetic vector potential.They use the temporal gauge while we adopt the Coulomb's gauge.What's more, they only provide a plain convergence for the finite element solutions and no error estimates are given.On the other hand, compared to the work in [29], they focused stable finite element discretization for stationary incompressible MHD equations and corresponding fast solvers.Although some convergence results of the finite element solution are established, no error analysis are presented.It is much more complicated and difficult to carry out the error analysis as we have to deal with the issues due to the nonlinear coupling of multiple variables.Therefore, more delicate analyses are needed for the error estimates.
The paper is organized as follows.In Section 2, we introduce some notations and present the energy estimate fo the MHD equations.In Section 3, we propose a fully discrete finite element method for the MHD equations.In Section 4, we state the main results of this work.In Section 5, we deduce some detailed proofs in Section 4. In Section 6, we present some numerical experiments to confirm the stability and accuracy of the scheme.In Section 7, we conclude with a few remarks.

Preliminaries
We start by introducing some notations and Sobolev spaces.Throughout the paper, the vector-valued functions and vector-valued function spaces in R 3 are denoted in boldface.As usual, the inner product and norm in  2 (Ω) are denoted by (•, •) and ‖•‖, respectively.Let  , (Ω) stand for the standard Sobolev spaces equipped with the standard Sobolev norms ‖•‖ , .For  = 2, we write   (Ω) for  ,2 (Ω) and its corresponding norm is ‖•‖  .For a given Sobolev space , we write   (0,  ; ) for the Bochner space.We will use the following notations for some spaces, In addition, the standard norm in (curl, Ω) is denoted by ‖•‖ curl .Here and what follows, we use  to denote generic positive constants independent of the discretization parameters, which may take different values at different places.
Then we present the CT-MHD model in a dimensionless form.Since div = 0 in Ω and  × = 0 on Γ, we invoke with the vector potential theorem in [15] to introduce a magnetic vector potential  with the Coulomb's gauge such that Putting ( 2) into (1) and eliminating  and , we can obtain the following CT-MHD model The system are complemented with the following initial and boundary conditions, where the initial values satisfy Let ,  0 ,  0 ,  0 = / 0 be characteristic quantities of length, time, magnetic field, and fluid velocity, respectively.We introduce the dimensionless variables as follows Then, the MHD system (4) can be written in a dimensionless form, where   =  0 / is the Reynolds number,   =  0 is the magnetic Reynolds number and )︀ is the coupling number between the fluid and the magnetic field.Note that we have introduced a Lagrange multiplier  with homogeneous Neumann boundary condition to deal with divergence-free constraint of div = 0. Indeed, the equation (5e) is equivalent to the following equation This is because if we apply divergence operator to (5e), we get ∆ = 0. From the boundary condition for , we deduce  = 0 in  1 /R.In this paper, we adopt (5e) since this formulation allows us easily deal with the divergence-free condition for .Similar ideas can be found in [29,36].Now we derive a weak formulation of (5).For convenience, we introduce some notations for function spaces  :=  1 0 (Ω),  :=  2 0 (Ω),  :=  0 (curl, Ω),  := (curl, Ω),  :=  1 (Ω)/R.
On the one hand, we note that the energy law describes the variation of the total energy caused by energy conversion.The total energy consists of the fluid kinetic energy 1  2 ‖‖ 2 and the magnetic field energy  2 ‖‖ 2 .
The dissipation stems from the friction losses  −1  ‖∇‖ 2 and the Ohmic losses  −1  ‖curl ‖ 2 .On the other hand, we see that there is no term related to the magnetic vector potential  in the energy law.This observation indicates that  is only a intermediate variable in this system.Even so, we still have some stability results for it.
Remark 2.2.Note that the boundary we considered for the magnetic induction and electric field is  ×  = 0 and  •  = 0.As to another frequently used set of boundary conditions for the magnetic induction and electric field,  •  = 0 and  ×  = 0, we need make some modifications on the boundary conditions for the magnetic field and magnetic vector potential.Invoking with the vector potential theorem in [15], we deduce the following boundary conditions for  and , Furthermore, our forthcoming analysis can be easily extended to this type boundary conditions.

A semi-implicit mixed FEM
This section is devoted to giving the fully discrete finite element scheme for the MHD equations.For the space discretization, let  ℎ be a quasi-uniform and shape-regular tetrahedral mesh of Ω.As usual, we introduce the local mesh size ℎ  = diam() and the global mesh size ℎ := max ∈ ℎ ℎ  .For any integer  ≥ 0, let   () be the space of polynomials of degree  on element  and define   () =   () 3 .To approximate the velocity and pressure (, ), we use the conforming finite element pair ( ℎ ×  ℎ ) ⊂ ( × ), which satisfies the discrete inf-sup condition, inf where   is a constant independent of mesh size ℎ.To approximate the magnetic field , we employ the conforming finite element space  ℎ ⊂  .To approximate the magnetic vector potential and Lagrange multiplier (, ), we adopt the conforming finite element pair where   is a constant independent of mesh size ℎ.We assume that the finite element spaces possess the following approximation properties ulteriorly.There exists an integer  ≥ 1 such that the following standard approximation properties hold for  ≥ 1, inf inf inf inf Many existing stable finite element pairs can be available in present setting [3,15,30].In the numerical experiments of Section 6, we choose Mini-element to discrete the velocity and pressure, the lowest order Nédélec edge element for the magnetic field and magnetic vector potential, and linear nodal element for the Lagrange multiplier.To be specific, where  1, () =  1 () span{ 1  2  3  4 },   ,  = 1, 2, 3, 4 are the barycentric coordinate functions on .
For convenience, we define the discrete kernel space of the divergence operator on  ℎ by Similarly, we introduce the discretely solenoidal function space on  ℎ by From the de-Rham sequence, we know that there exists a Lagrange finite element space  ℎ ⊂  1 0 (Ω) such that ∇ ℎ ⊂  ℎ .For the choice of (17), Then the discretely solenoidal function space on  ℎ is defined by Further, there holds the following Poincaré's inequality [1,18,21,34,36], In order to deal with the convection term in (5a), we define the following trilinear form For the time discretization, let {  =  :  = 0, 1, • • • ,  },  =  /, be an equidistant partition of the time interval [0,  ].Given a function (, ), the full discrete approximation to (,   ) will be denoted by   ℎ .For  ≥ 1 and any function , define     =   − −1

𝜏
. With the discrete spaces and above notions, a fully discrete finite element scheme for problem (5) reads as follows.For any  ≥ 1, we find At the initial time step,  0 ℎ =  ℎ  0 and  0 ℎ = ℱ ℎ  0 , where  ℎ and ℱ ℎ are projection operators defined in Section 5. Finally, we give several remarks concerning the proposed scheme.Remark 3.1.From (19e), we know that the discrete magnetic vector potential is divergence-free weakly.In fact, the discrete solution for the magnetic field also satisfies the weakly divergence free property.For any  ℎ ∈  ℎ , by taking  ℎ = ∇ ℎ in (19c), we obtain (    ℎ , ∇ ℎ ) = 0, namely, . Thus, we do not need to add a Lagrange multiplier in the magnetic equation as in [29,36].Remark 3.2.After obtaining the magnetic field   ℎ and the magnetic vector potential   ℎ , we can define the current density as   ℎ := curl   ℎ and the magnetic induction as   ℎ := curl   ℎ .Thus the discrete current density and magnetic induction are divergence-free exactly.
Remark 3.3.Compared with the semi-implicit scheme for the -based model in [13,19], the extra cost is the computation of a saddle problem for the magnetic vector potential and Lagrange multiplier at each time step.Moreover, we have decoupling them from other variables by lagging the magnetic field in (19d), which makes the implementation simpler and the scheme more efficient.In the future, we will consider using the regularization strategy proposed in [11] with  =  to remove the Lagrangian multiplier  for the magnetic vector potential  while keeping the numerical accuracy.
Remark 3.4.As for the detailed implementation process of (19), we can first solve the discrete problem (19d) and (19e) for (, ) and then solve the discrete problem (19a)-(19c) for (, , ).Thus, there exits a similar fully discrete scheme, whose implementation process is really just the opposite.A brief discussion about this scheme is laid down in Appendix A.

Main results
In this section, we state the main results of this paper.First, we establish the unique solvability of the proposed scheme (19) is the following theorem.
Then, we present the stability of the discrete problem (19) in the following theorem.
For the error estimates, we assume that the exact solution of MHD system (5) uniquely exists and the unknowns have following regularity property: where  > 1 2 .Under this assumption, our main error estimate result can be summarized as follows.Theorem 4.3.Let (, , , , ) be the exact solution of (5) with the regularity (27) holds.Let (  ℎ ,   ℎ ,   ℎ ,   ℎ ,   ℎ ) be the numerical solution of the discrete system (19).Then we have the following error estimates for all with  = min{, } and  is independent of the discrete parameters  and ℎ.If we further assume that   ∈  ∞ (0,  ;  2 (Ω)), there holds that The detailed proof of this theorem will be given in Section 5.
Remark 4.4.In Theorem 4.3, the regularity assumption for the exact solutions is necessary to deliver the error estimates.Compared with the existing works [19,42], the regularity for the exact solutions is relative weak and may be a reasonable hypothesis in such a setting of the general polyhedral domain with Lipschitz boundary.
The interested readers can refer to [14] for regularity of classical MHD model and [6,7] for Maxwell operator.

Error analysis
In this section, we carry out a rigorous error analysis for the proposed scheme.

Auxiliary estimates
In this subsection, we gather the necessary tools for the error estimates.First we present the classical and discrete Sobolev inequalities needed for the error estimates in the next subsection.
For  ∈   (Ω) with  > 1 2 , we have Next we define some projections of the unknowns and gather their approximation properties.For the fluid pair (, ), we define the Stokes projection ( ℎ , For the magnetic field , we utilize the Fortin operator ℱ ℎ :  →  ℎ , see [2,43].Given  ∈  , we define For the magnetic vector potential and Lagrange multiplier (, ), we define the Maxwell projection From [10,15,35,36,43], we have the following approximation property result for the projections.
As a consequence of the above result, we have that the initial errors satisfy Finally, we will frequently use the following discrete version of the Gronwall lemma [24].
Lemma 5.3.Let   ,   ,   and   be non-negative numbers with  ≥ 0 such that where  and  are two positive constants.Then

Error estimates
In this section, we derive error estimates for the fully discrete scheme of the MHD equations.In order to facilitate the subsequent error estimates, we define the following notations here and hereafter where the errors are decomposed into approximation errors    and discretization errors    .The splittings are performed with the Stokes projection, the standard  2 -projection and the Fortin operator, First we notice that the exact solutions of the system (5) satisfy the following variational equations at   for all ( ℎ ,  ℎ ,  ℎ ,  ℎ ,  ℎ ) where the truncation errors    ,    and    are defined by In the following lemma, we estimate the residual terms    ,    and    .Lemma 5.4.Suppose that (27) is satisfied, then we have Proof.From the definition of    , using Hölder's inequality, Sobolev inequalities in Lemma 5.1 and ( 27), we arrive at Thus we have Similarly, which means that For the term    , we estimate it as which implies that This completes the proof.
Subtracting the numerical system (19) from the above system (35) and using the definitions of projections ( 31)- (33), we obtain the following error equations for all ( ℎ , where the abbreviated terms are gathered as After the above preparation, we are now in a position to prove the error estimates. Proof of Theorem 4.3.Taking ( ℎ ,  ℎ ,  ℎ ) = (   ,    ,    ) in (38a)-(38c) and adding the three equations together, we obtain We next bound the terms on the right hand side of the above identity one by one.For the first two terms on the right hand side of (39), using Cauchy-Schwarz inequality, Young inequality and (18), we have Similarly, for the third and forth terms on the right hand side of (39), we get For the fifth term, we rearrange it as follows, The terms in the last step can be bounded using Hölder's inequality, Sobolev inequalities in Lemma 5.1 and the approximation properties of the projections in Lemma 5.2 as, Thus, we conclude that For the last two terms on the right hand side of (39), we rearrange them as follows: Next we will estimate  1 +  3 and  2 +  4 separately.For the terms  1 +  3 , we have where we have used the fact ( × , ) + ( × , ) = 0 to cancel the last two terms out in last equality.Using Hölder's inequality and Young's inequality, we get For the terms  2 +  4 , we insert the following identity into them, Applying Hölder's inequality, Young's inequality and the stability of the projections in Lemma 5.2, we have To deal with the term ‖curl    ‖ on the right hand side of (41), we take ( ℎ ,  ℎ ) = (   ,    ) in (38d) and (38e) to get Thus, we have Hence, we derive that Combining all the above estimates, we arrive at Before proceeding with the analysis, we notice that the following estimates hold )︁ .
Multiplying 2 on the estimate (45), summing over  = 1, • • • , , using the fact  0  =  0  = 0, the regularity assumption (27) and the approximation properties of the projections in Lemma 5.2, we have where Ĉ is given by )︁ By the discrete Gronwall inequality, we have We complete the proof of the error estimates for  and  in (28) by applying the triangle inequality, the approximation properties of the projections in Lemma 5.2 together with above estimates.With the above estimates for  and , we invoke with (42) and the approximation properties of the projections in Lemma 5.2 to get ‖curl Using (37) with the triangle inequality, we obtain the error estimate (29) for .With a slightly stronger regularity assumption with   ∈  ∞ (0,  ;  2 (Ω)), we have With this result, we can regain the full order of  in (30) by using (49) and triangle inequality.This completes the proof.
Remark 5.5.For the error estimate for  in the scheme (19), the truncation term    enters into the error equations, which further leads to the estimate (49).Hence, with the regularity (27), we only get the suboptimal convergence order in time for  with half an order lost.With a slightly stronger regularity assumption with   ∈  ∞ (0,  ;  2 (Ω)), we improve the estimate of    and regain the full order of  .But this issue does not arise during the error estimate for the scheme (A.1), we refer to Appendix A for more details.Here we use the notation   (•) to indicate that the constant   depend on the "•",  = 1, 2. Clearly, one can see that the constant  on the right-hand side of (48) heavily depends on the the Reynolds number, the magnetic Reynolds number and the coupling number.Therefore, the error estimates obtained will deteriorate with these physical parameters increase.

Numerical experiments
In this section, we give some numerical experiments to validate the accuracy and stability of the scheme.All simulations are run using the finite element software FreeFEM [20].
Since the exact solutions are linear or constant in spatial variables, the errors only come from the discretization of the time variable.We fix a tetrahedron mesh with ℎ = 1/4 and the terminal time  = 1.The numerical results with different time steps are listed in Table 1.We observe that the expected first-order convergence is obtained for the scheme, which is consistent with the error estimates in Theorem 4.3.Example 6.2 (Convergence rate for full discretization).This example is to test the convergence rates for the fully discretization.The computational domain is taken as Ω = (0, 1) 3 , the physical parameters are set by   =   =  = 1 and the terminal time is given by  = 1.The right-hand sides, initial and boundary conditions are chosen so that the exact solutions are given by  = (cos() sin(), sin() exp(−), sin() cos()),  = sin( +  + ) cos(),  = 0,  = (sin() cos(), cos() sin(), cos() exp(−)),  = (cos() cos(), sin() exp(−), sin() sin()).The initial mesh size and time step are set by ℎ = 1/2 and  = 1/5.We fist refine the time step and the mesh size at the same time such that  = ℎ to test the convergence of  1 -error for , (curl)-errors for  and , the  2 -errors for ,  and .Next, to test the convergence of  2 -error for , we refine the time step and the mesh size simultaneously such that  = ℎ 2 .The corresponding results are displayed in Table 2.We can see that the discrete solutions have the asymptotic behaviors, In other words, the expected convergence rates of the proposed numerical scheme can be observed, which agrees with Theorem 4.3.The unconditional and optimal error estimates will be our further work.
Example 6.3 (Unconditional energy stability).This example is to test the unconditional energy stability of the proposed scheme.We take the domain as Ω = (0, 1) 3 and set the initial conditions for  and  to be With the given data, we test the energy stability of our scheme under different physical parameters and time steps on a fixed mesh with ℎ = 1/16.Denote the discrete energy at   by E 1 shows the time evolution of the discrete energy with the finial time  = 5.It can be observed that all energy curves decrease monotonically, which means that our scheme is unconditionally energy stable.where ℎ is the mesh size.Let   = (0, 0, −), the boundary conditions are set by To compare our results with those in [29], the physical parameters are set by   = 100,   = 10 and  = 1.In the computation, we set the mesh size ℎ = 1/32, the time step  = 0.01, and the finial time  = 5.In Figure 2, we display the streamlines of velocity on the three cross-sections  = 0.5,  = 0.5 and  = 0.5, respectively.They show clearly the vortex structure of the fluid.Note that our result on the cross-section  = 0.5 is in agreement well with the previous result in [29], which studied the steady CT-MHD model.Figure 3 shows the distributions of | ℎ | and Figure 4 shows the magnetic vector potential  ℎ on the three cross-sections  = 0.5,  = 0.5 and  = 0.5, respectively.From these figures, we can know how the fluid influences the electromagnetic fields.

Conclusions
In this paper, we study a fully discrete finite element scheme for the CT-MHD model.The scheme is designed by combining an Euler semi-implicit scheme for the temporal discretization and mixed finite elements for the spatial discretization.A novel feature is that it preserves the divergence-free property exactly for the magnetic induction and current density on the discrete level.Error estimates for the velocity, magnetic field and magnetic vector potential are further established under the low regularity hypothesis for the exact solutions.Some numerical experiments are presented to verify the theoretical analysis and show the performance of the proposed scheme.
Note that we use a semi-implicit backward Euler scheme in time, which is first-order accurate in the temporal direction.Apparently, one replace it by higher-order time integrator with some subtle implicit-explicit treatments for nonlinear and coupling terms, such Crank-Nicholson scheme and second-order backward difference formula   (BDF2).For instance, a semi-implicit scheme based on BDF2 is to find where the second-order difference operator  2  and extrapolated operator ˆare defined by , ζ := 2 −1 −  −2 .
For  = 1, we can use the semi-implicit Euler scheme (19) to compute . Provided enough temporal regularity, one can easily obtain the error estimate by a similar analysis.We believe that their analysis shall not present difficulties beyond those already encountered in this work.In the future, robust and fast solvers for this scheme will be considered.

Appendix A.
In this section, we give a similar fully discrete scheme to (19) and present the results of error analysis in a nutshell.Different from the treatment skills in (19), here we lag the magnetic vector potential in (19a)-(19c) to decouple it from other variables.Hence, a new fully discrete finite element scheme for problem (5)  )︀ .
Next, we state the unique solvability and stability of the scheme (A.1) is the following theorem, which mimics Theorems 4.1 and 4.2.
In the proof of the error estimate for  in the scheme (A.1), the truncation term    vanishes in the error equations due to the treatment of  and  in (A.1d) at the same time level.As a result, we can derive the optimal convergence rate in time straightway.But it's a pity that we have introduce a slightly stronger regularity assumption with curl   ∈  2 (0,  ;  2 (Ω)) by comparing (27) with (A.7).

Example 6 . 4 (
Driven cavity flow).In this example, we simulate the benchmark problem of driven cavity flow.The problem models the flow of a conducting fluid driven by the movement of the lid at the top of the

Table 1 .
Errors and convergence rates for time discretization.Remark 5.6.Invoking with (45)-(47) and Lemma 5.3, we derive carefully to get a more elaborate expression for the constant on the right-hand side of (48),  ,   ,    2 )︀ .

Table 2 .
Errors and convergence rates for the fully discrete scheme.