A decoupled Crank-Nicolson time-stepping scheme for thermally coupled magnetohydrodynamic system

Article History: Received 20 March 2016 Accepted 11 August 2017 Available 23 October 2017 Thermally coupled magneto-hydrodynamics (MHD) studies the dynamics of electro-magnetically and thermally driven flows, involving MHD equations coupled with heat equation. We introduce a partitioned method that allows one to decouple the MHD equations from the heat equation at each time step and solve them separately. The extrapolated Crank-Nicolson time-stepping scheme is used for time discretization while mixed finite element method is used for spatial discretization. We derive optimal order error estimates in suitable norms without assuming any stability condition or restrictions on the time step size. We prove the unconditional stability of the scheme. Numerical experiments are used to illustrate the theoretical results.


Introduction
Thermally coupled magneto-hydrodynamics has many applications including in electromagnetic pumping design [35], electromagnetic filtration [4], contact-less electromagnetic stirring [32] and damping convective flow in metal-like melt [34].Magnetohydrodynamics in general has broad applications including fusion [19], underwater propulsion [18], nuclear reactors [13], metallurgy [1,2,11,31] and astrophysics [30].In all of these applications, qualitative and quantitative understanding of the dynamics is important to achieve optimal operating conditions.This has led to considerable research efforts over the past three decades into the development of theoretical, see e.g [16,24,26,27,29] and efficient and accurate computational techniques, see e.g.[8,9,20,21] for MHD equations.Majority of the numerical analysis work done on the equations has been for steady state equations.In [17,23,25,33], time stepping schemes for unsteady MHD equations have been analyzed.However, these work consider MHD equations where thermal effects are negligible.Thermally coupled MHD equations model a complex flow phenomena which is in general three dimensional, highly nonlinear and represents multi-physics.
In this work, we propose and analyze a decoupled time stepping scheme for the thermally coupled MHD equations.It uses a semi-implicit Crank-Nicolson scheme, which combines an implicit treatment of the second derivative terms, a semi-implicit second order extrapolation of the nonlinear convective terms and an explicit treatment of the temperature coupling term in the Navier-Stokes equations.The proposed scheme solves the MHD equations and the heat equation separately in each time step (without iteration) allowing the possibility of optimizing the subproblem's respective physics.We show unconditional stability of the scheme and provide a complete error analysis for fully discrete scheme using finite element spatial discretization.

*Corresponding Author
The remaining of the paper is organized as follows: The continuum problem and some preliminaries are presented in Section 2. In Section 3, we present the decoupled time-stepping scheme and analyze its stability, accuracy and convergence.Finally, we present a numerical example that illustrates our theoretical results.

Continuum problem and preliminaries
To begin with,we present some notations and basic results that will be used throughout the article.

Continuum problem
The non-dimensional Boussinesq equations describing thermally coupled MHD equations are (see for e.g.[15]) in (0, T ], where T denotes time and Ω ⊂ R I d (d = 2, 3) a bounded region with Lipschitz-continuous boundary Γ.Moreover the different fields appearing in the equations are u(x, t) the fluid velocity, B(x, t) the magnetic field, θ the temperature, p(x, t) the pressure, f the source and i 3 the unit basis vector.The non-dimensional numbers that appear in the MHD equations are S := P r B P r θ H 2 , the Hartman number H, the Rayleigh number Ra, the thermal Prandtl number P r θ and the magnetic Prandtl number P r B .The MHD system we consider is supplemented with the initial conditions u(x, 0) = u 0 (x) , θ(x, 0) = θ 0 (x) and along with the boundary conditions

Function spaces
For a Banach space X, we denote by L p (0, T ; X) the time-space function space endowed with the norm w L p (0,T ;X) := T 0 w p X dt We will often use the abbreviated notation L p (X) := L p (0, T ; X) for convenience.The symbol C([0, T ]; X) denotes the set of continuous functions u : [0, T ] → X endowed with the norm u C(0,T ;X) := max 0≤t≤T u(t) X .For any integer k ≥ 1, let W k,p (Ω) be the Sobolev space of functions in L p (Ω) with derivatives up-to the k th order endowed with the norm φ m,p := where We denote by H k (Ω) the space W k,2 (Ω), when p = 2, and drop the subscripts p(= 2) in referring to the norm in H k (Ω).Moreover, we will use the following simplified norm notations: We introduce the time discrete space l p (Z) associated with L p (0, T ; Z); l p (Z) is the space of Z-valued sequences w := {w n ; n = 1, . . ., N } with norm • l p (Z) defined by For later purposes, we recall the inequality the Poincaré inequality the Gagliardo-Nirenberg interpolation inequality [3] We define the explicitly skew-symmetrized trilinear forms and the trilinear form Notice that the trilinear form d(•, •, •) is skewsymmetric with respect to the first and last arguments, i.e., d(B, C, v) = −d(v, C, B).
We end this section with a result regarding the existence and uniqueness of solutions to the initialboundary value problem (1)-( 3) whose proof can be furnished by using Galerkin approximations, a-priori estimates and compactness methods.
Proposition 1. Assume that the given functions f , g, k, q, q, u 0 and B 0 satisfy . In two-spatial dimension (d = 2), these solutions are unique.

Properties of finite element spaces and projections
In order to keep the exposition simple, we restrict our attention to convex polyhedral domains.Let T h be a family of subdivisions (e.g.triangulation) of Ω ⊂ R d satisfying Ω = ∪ K∈T h K so that diameter(K) ≤ h and any two closed elements K 1 and K 2 ∈ T h are either disjoint or share exactly one face, side or vertex.Suppose further that T h is a shape regular and quasi-uniform triangulation.That is, there exists a constant C > 0 such that the ratio between the diameter h K of an element K ∈ T h and the diameter of the largest ball contained in K is bounded uniformly by C, and h K is comparable with the mesh size h = max K∈T h h K for all K ∈ T h .For example, T h consists of triangles for d = 2 or tetrahedra for d = 3 that are nondegenerate as h → 0. We choose families of finite dimensional spaces , parameterized by a parameter h such that 0 < h < 1.Let g h , q h and q h be approximations of g, q and q, respectively, such that there exists We also define the discretely divergence free space is given by We make the following assumptions on the finite dimensional subspaces X h , Y h , Z h and Q h : Assumption A1.We have the approximation properties: there exists an integer k and a constant C, independent of h, v, B, θ and r, such that inf For every r h ∈ Q h , there exists a nonzero function v h ∈ X h and β > 0 such that with an inf-sup constant β > 0 that is independent of the mesh size h.Assumption A3.For any integers l and m (0 ≤ l ≤ m ≤ 1) and any real numbers p and q (1 ≤ p ≤ q ≤ ∞) it holds that There are many conforming finite element spaces satisfying the assumptions (A1)-(A3).One may choose, for example, the Taylor-Hood element pair for the velocity and pressure (i.e, piecewise quadratic polynomial for velocity and piecewise linear polynomial for pressure), and piecewise quadratic polynomials for the magnetic field and temperature.Then, hypothesis (A1)-(A3) hold with k = 2.We define Stokes, Maxwell and Ritz projections as follows: Given (u, p) ∈ H 1 (Ω)×L 2 0 (Ω), θ ∈ H 1 (Ω) and B ∈ H 1 (Ω), we define the Stokes projection (P s h u, P s h p) ∈ X h,g h × Q h as the solution of the problem the Maxwell projection P m h B ∈ Y h,q h as the solution of the problem and the Ritz projection P r h θ ∈ Z h, q h as the solution of the problem We have the following convergence and boundedness results for these projections.
Lemma 1. Suppose that assumptions (A1)-(A2) hold with a positive integer k, and that (u, p) ∈ and B ∈ H k+1 (Ω).Then, for any h ∈ (0, h 0 ] the Stokes projection (P s h u, P s h p) of (u, p) satisfies Proof.The proof of ( 8)- (10) follows by the regularity properties of the Stokes, Maxwell and Ritz projections and by duality argument.In order to prove (11)-( 13), we first notice that Gagliardo-Nirenberg's inequality yields 2 .Therefore the approximation properties ( 8)- (1) together with Agmon's inequality yield the desired result.
Let ∆t denote the step size for t so that t n = n∆t, n = 0, 1, 2, . . ., N .For notational convenience, we denote Moreover, let P s h u be the Stokes projection of u, P m h B the Maxwell projection of B and P r h θ the Ritz projection of θ.If assumptions (A1)-(A2) hold with a positive integer k, then Proof.The proof of (i)-(iii) follows by Taylor expansion with integral remainder whereas the proof of (iv)-(vi) follows as a consequence of Lemma 1.
We will need the following well known discrete Grönwall lemma.Lemma 3. (Discrete Grönwall lemma) Let d, ∆t, {a n } n≥0 , {b n } n≥0 , {c n } n≥0 , and {d n } n≥0 be nonnegative numbers such that A proof of this result can be found, for e.g, in [12].

Decoupled Crank-Nicolson time-stepping scheme
We discretize the system (1) by Crank-Nicholson scheme in time and Galerkin finite element in space.The time discretization combines an implicit treatment of the second derivative terms, a semi-implicit second-order extrapolation for the nonlinear convective terms and explicit treatment of the temperature coupling term in the Navier-Stokes equations.

Stability analysis
In this section, we demonstrate the unconditional energy stability of the decoupled scheme proposed in Section 2. We first recall a few basic facts and some notation that are needed below.Let us define the discrete trace spaces of X h , Y h and Z h by [10,28] .Similarly, we can define discrete extension operators E h and

Then there exists a discrete extension operator
In order to prove, we first define suitable boundary extensions.
We make the following assumptions about the ex- Assumption A4.

The extension operators satisfy
Suppose assumption (A4) holds and let {( < ∞ and ) and using the skew-symmetry of )) )) Let us next bound each term on the right-hand side of (15) 1 except the last two.The first five terms can be estimated using Cauchy/Duality and Young's inequalities to obtain We estimate A n 6 and A n 7 using Hölder's, Gagliardo-Nirenberg and Young's inequalities as follows ) 4   1 and ) 4   1 Collecting these estimates in (15) 1 , we obtain We employ similar arguments to bound the terms on the right-hand-side of (15) 2 and (15) 3 to obtain and Finally we estimate the last terms in ( 16)- (18) using assumption (A4) and Young's inequality to obtain where ǫ is a suitably chosen positive constant.Employing these estimates in ( 16)-( 18), we obtain 1 ) , Now summing each of the inequalities in (20) from n = 2 to m, using the skew symmetry of d(•, •, •) and the telescoping property, we obtain that for some constant M > 0 by the assumptions.The required stability bound follows by setting ( ) and applying triangle inequality.

Error analysis
In this section we discuss the accuracy and convergence of the decoupled Crank-Nicolson scheme.In the subsequent analysis, we will assume the boundary data is independent of time for simplicity.
Theorem 2. Suppose that the assumption (A1)-(A3) hold with a positive number h 0 and a positive integer k, that the solution (u, B, p, θ) of ( 1)-( 3) ) and that the initial conditions (u i h , B i h , θ i h ) , i = 0, 1 satisfy Then, for any h ∈ (0, h 0 ] the approximate solutions (u h , B h , θ h ) of ( 14) satisfy the following error estimates for some constant C independent of the mesh size h and time step ∆t.

Proof. Let (P
)) and Using the definition of Stokes, Maxwell and Ritz projections, we obtain the basic error equations of the method for all We next split the nonlinear terms < ℵ n h , v h >, < ℵ n h , φ h > and < ℵ n h , ψ h > on the right-hand side of ( 22) into several terms as follows: )) )) , ψ h ) =: ), e ), v h ) .
Therefore, we have Estimating ℵ n 1 − ℵ n 4 similarly, we obtain Employing ( 24)-( 29) into ( 23) and using Young's inequality, we obtain ), ) )) + (e ) , where We next add the three equations in (30) and use the identity ) where From the assumptions on the solution (u, p, B, θ) it holds that ∆t Therefore summing (31) from n = 1 to m and the discrete Grönwall inequality (Lemma 3), we have that The required error estimate now follows from (33) and triangle inequality.
Theorem 3.Under the assumptions in Theorem 2, the approximate pressure p h of ( 14) satisfies some constant c independent of mesh size h and time step ∆t.
Proof.From ( 22) 1 and the inf-sup condition it holds that We start estimating ) X h * below.First, by Hölder's and Gagliardo-Nirenberg inequalities, we obtain Before estimating the other two terms, notice that by the inverse estimate (Assumption (A3)) and (33), we obtain Similarly, we can show Therefore, by Hölder's, Gagliardo-Nirenberg inequalities and ( 35)-(36), we obtain Estimating other terms in (34)  ) + e n+1/2 1h The required error estimate now follows from last inequality by using Theorem 2 and triangle inequality.
The error estimate for the pressure in the previous theorem can be improved under stronger regularity properties of the solution.To this end, we next derive optimal order error estimates for the time derivatives of velocity, magnetic field and temperature.
In addition, assume the initial conditions Then for any h ∈ (0, h 0 ] the approximate velocity u n h , magnetic field B n h and temperature θ n h satisfy for some constant c independent of the mesh size h and time step ∆t Moreover, we have for some constant c independent of the mesh size h and time step ∆t. Proof.Putting v h = D(e n 1h ), φ h = D(e n 3h ), ψ h = D(e n 4h ) into (22) and splitting the nonlinear terms as in the proof of Theorem 2, we obtain 2 )] Let us start estimating < ℵ n i , D(e n 1h ) > for i = 1, . . ., 14. First using Hölder's inequality and Gagliardo-Nirenberg inequality, we obtain From the inverse inequality (Assumption (A3)) and Gagliardo-Nirenberg inequality, it follows that Using (39), we estimate < ℵ Alternatively, we can estimate < ℵ n 6 , D(e n 1h ) > as follows ), ), ) 1 e n 1h 1 e n−1 1h 1 . (41) Combining ( 40) and (41), we have where Estimating other terms as before, we obtain .
Let us next start estimating ℵ 1 − ℵ 6 .First, we rewrite them using integration by parts formula and then we estimate them using Hölder's inequality and Gagliardo-Nirenberg inequality .
(47) Summing (44) from n = 1 to m and the assumptions about initial conditions (u i h , B i h , θ i h ), i = 0, 1, we obtain Therefore using (49) in (37), we obtain the required estimate.

Numerical results
In this section, we present a numerical example to illustrate the theoretical results of the previous section.We set Ω := (0, 1) × (0, 1) and choose the standard piecewise quadratic finite space for approximating the magnetic field and temperature.We also choose the Taylor-Hood element pair, i.e., continuous piecewise-quadratic and continuous piecewise linear finite element space for the fluid velocity and pressure approximations, respectively.Uniform triangular meshes are created by first dividing the rectangular domain Ω into identical small squares and then dividing each square into two triangles.We set the exact solutions to u = ((y + y 2 )e −t , (x + x 2 )e −t ) B = ((sin(y) + y)e −t , (sin(x) + x 2 )e −t ) p = (x + y)e −t θ = (1 + xy)e −t .
The right-hand side data in the MHD system, initial conditions and boundary conditions are then chosen correspondingly.For simplicity, we set the parameters P r θ , S, P r B , Ra equal to 1.0.In order to determine the order of convergence α with respect to the time step ∆t, we fix the spatial spacing h and use the following approximation

Table 1 .
set of values of α are listed in Table4.1 with a fixed spacing h = 1/32 and varying time step ∆t = 1/20, 1/40, 1/80, 1/160, 1/320, which clearly suggest the concerned orders of convergence in time are all O(∆t 2 ) for the decoupled scheme.Thus, the numerical experiments clearly suggest that the orders of convergence in time in error estimates in Theorem 2 for the L 2 − norm of u, B and θ are optimal.Convergence order of O(∆t α ) of the partitioned scheme at time t N = 1.0, with the fixed spacing h = α ≈ log 2 v h,∆t (x, t N ) − v h, ∆t 2 (x, t N ) v h, ∆t 2 (x, t N ) − v h, ∆t 4 (x, t N ) .(50)A