Loading [MathJax]/jax/output/SVG/jax.js
Research article

A modified characteristics projection finite element method for unsteady incompressible Magnetohydrodynamics equations

  • This paper provides a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations. In this method, modified characteristics finite element method and the projection method will be combined for solving the unsteady incompressible MHD equations. Both the stability and the optimal error estimates both in L2 and H1 norms for the modified characteristics projection finite element method will be shown. In order to demonstrate the effectiveness of our method, we will present some numerical results at the end.

    Citation: Shujie Jing, Jixiang Guan, Zhiyong Si. A modified characteristics projection finite element method for unsteady incompressible Magnetohydrodynamics equations[J]. AIMS Mathematics, 2020, 5(4): 3922-3951. doi: 10.3934/math.2020254

    Related Papers:

    [1] Lanyin Sun, Kunkun Pang . Numerical solution of unsteady elastic equations with C-Bézier basis functions. AIMS Mathematics, 2024, 9(1): 702-722. doi: 10.3934/math.2024036
    [2] Hyam Abboud, Clara Al Kosseifi, Jean-Paul Chehab . Stabilized bi-grid projection methods in finite elements for the 2D incompressible Navier-Stokes equations. AIMS Mathematics, 2018, 3(4): 485-513. doi: 10.3934/Math.2018.4.485
    [3] Zhengang Zhao, Yunying Zheng, Xianglin Zeng . Finite element approximation of fractional hyperbolic integro-differential equation. AIMS Mathematics, 2022, 7(8): 15348-15369. doi: 10.3934/math.2022841
    [4] Jinghong Liu . $ W^{1, \infty} $-seminorm superconvergence of the block finite element method for the five-dimensional Poisson equation. AIMS Mathematics, 2023, 8(12): 31092-31103. doi: 10.3934/math.20231591
    [5] Chunjuan Hou, Zuliang Lu, Xuejiao Chen, Xiankui Wu, Fei Cai . Superconvergence for optimal control problems governed by semilinear parabolic equations. AIMS Mathematics, 2022, 7(5): 9405-9423. doi: 10.3934/math.2022522
    [6] Jie Liu, Zhaojie Zhou . Finite element approximation of time fractional optimal control problem with integral state constraint. AIMS Mathematics, 2021, 6(1): 979-997. doi: 10.3934/math.2021059
    [7] Mingliang Liao, Danxia Wang, Chenhui Zhang, Hongen Jia . The error analysis for the Cahn-Hilliard phase field model of two-phase incompressible flows with variable density. AIMS Mathematics, 2023, 8(12): 31158-31185. doi: 10.3934/math.20231595
    [8] Haifeng Zhang, Danxia Wang, Zhili Wang, Hongen Jia . A decoupled finite element method for a modified Cahn-Hilliard-Hele-Shaw system. AIMS Mathematics, 2021, 6(8): 8681-8704. doi: 10.3934/math.2021505
    [9] Zuliang Lu, Xiankui Wu, Fei Huang, Fei Cai, Chunjuan Hou, Yin Yang . Convergence and quasi-optimality based on an adaptive finite element method for the bilinear optimal control problem. AIMS Mathematics, 2021, 6(9): 9510-9535. doi: 10.3934/math.2021553
    [10] Lingling Sun, Hai Bi, Yidu Yang . A posteriori error estimates of mixed discontinuous Galerkin method for a class of Stokes eigenvalue problems. AIMS Mathematics, 2023, 8(9): 21270-21297. doi: 10.3934/math.20231084
  • This paper provides a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations. In this method, modified characteristics finite element method and the projection method will be combined for solving the unsteady incompressible MHD equations. Both the stability and the optimal error estimates both in L2 and H1 norms for the modified characteristics projection finite element method will be shown. In order to demonstrate the effectiveness of our method, we will present some numerical results at the end.


    Magnetohydrodynamics(MHD) is a subject that studies the motion of conductive fluids in magnetic fields. It is mainly used in astrophysics, controlled thermonuclear reactions and industry. The system of MHD equations is a coupled equation system of Navier-Stokes equations of fluid mechanics and Maxwell's equations of electrodynamics (see[1,2]). In this paper, we consider the unsteady incompressible MHD equations as follows:

    {ut+(u)uνu+μB×curlB+p=f,xΩ,u=0,xΩ,Bt+Rm1curlcurlBcurl(u×B)=g,xΩ,B=0,xΩ, (1.1)

    where u is velocity, ν kinematic viscosity, μ magnetic permeability, p hydrodynamic pressure, B magnetic field, f body force, and g applied current. The coefficients are the magnetic Reynolds number Rm. As in [3], we assume that Ω is a convex domain in Rd, d=2,3 and t[0,T]. Here, ut=ut, Bt=Bt. The term u×B explains the effect of hydrodynamic flow on current. The magnetic Reynolds number Rm is moderate or low (from 102 to 1) (see[4]). For the purposes of simplicity, we choose Rm=1 in our numerical analysis.

    System (1.1) is considered in conjunction with the following initial boundary conditions

    u(x,0)=u0(x),B(x,0)=B0(x),xΩ,u=0,Bn=0,n×curlB=0,xΩ,t>0,

    where n denotes the outer unit normal of Ω.

    Over that last few decades, there has been a lot of works devoted on the numerical solution of MHD flows. In [5], some mathematical equations related to the MHD equations have been studied. In [6], a weighted regularization approach was applied to incompressible MHD problems. Reference [7,8] decoupled the linear MHD problem involving conducting and insulating regions. A mixed finite element method for stationary incompressible MHD equations was given in [9]. In order to avoid in-sup condition, Gerbeau [10] presented a stabilized finite element method and Salah et al. [11] proposed a Galerkin least-square method. In [12], a two-level finite element method with a stabilizing sub-grid for the incompressible MHD equations has been shown.

    Recently, a convergence analysis of three finite element iterative methods for 2D/3D stationary incompressible MHD equations has been given by Dong and He (see[13]). In order to continue the in-depth explore of three-dimensional incompressible MHD system, an unconditional convergence of the Euler semi-implicit scheme for the three-dimensional incompressible MHD equations was studied by He [14]. In [15], a two-level Newton iteration method for 2D/3D steady incompressible MHD equation was proposed by Dong and He. In addition, Yuksel and Ingram have discussed the full discretization of Grank-Nicolson scheme at small magnetic Reynolds numbers (see[16]). The defect correction finite element method for the MHD equations was shown by Si et. al. [17,18,19]. In [20], a decoupling penalty finite element method for the stationary incompressible MHD equation was presented. To further study this method, we can refer to [21,22,23].

    In this paper, a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations will be given. In this method, modified characteristics finite element method and the projection method will be combined for solving the incompressible MHD equations. Both the stability and the optimal error estimates both in L2 and H1 norms for the modified characteristics projection finite element method will be derived. In order to demonstrate the effectiveness of our method, we will present some numerical results at the end of the manuscript. The contents of this paper are divided into sections as follows. To obtain the unconditional stability and optimal error estimates, in section 2, we introduce some notations and present a modified characteristics projection algorithm for the MHD equations. Then, in section 3, we give stability and error analysis and show that this method has optimal convergence order by mathematical induction. In order to demonstrate the effectiveness of our method, some numerical results are presented in the last section.

    In this section, we introduce some notations. We employ the standard Sobolev spaces Hk(Ω)=Wk,2, for nonnegative k with norm vk=(k|β|=0Dβv20)1/2. For vector-value function, we use the Sobolev spaces Hk(Ω)=(Hk(Ω))d with norm vk=(di=1vi2k)1/2(see [24,25]). For simplicity, we set:

    X=H10(Ω)2,X0={vX:v=0},M=L20(Ω)={qL2(Ω)d:Ωq(x)dx=0},H={vL2(Ω)d:v=0},W={ψ(H1(Ω))d:ψn|Ω=0},W0={ψW:ψ=0},

    which are equipped with the norms

    vX=v0,qM=q0,ψW=ψ1.

    In this paper, we will use the following equality

    curl(ψ×B)=BψψB,ψ,BW0.

    Then, we recall the following Poincaré-Friedrichs inequality in W: there holds

    c0Ccurlc0,cW,

    with a constant C>0 solely depending on the domain Ω; (see[26]). Here, we define the following trilinear form a1(,,) as follows, for all uV, w,vX,

    a1(u,v,w)=(uw+12(u)w,v)=12(uw,v)12(uv,w)

    and the properties of the trilinear form a1(,,) [27] are helpful in our analysis

    a1(u,v,w)=a1(u,w,v),  uX0,v,wX,|a1(u,v,w)|N2u0(v0wL+vL6wL3),  uL2(Ω)2, vX,wL(Ω)2X,|a1(u,v,w)|N2(uLv0+uL3vL6)w0,  uL(Ω)2X, vX,wL2(Ω)2,v0γ0v0,vL6Cv0,wL4C0w0,  wX,

    where N>0 is a constant, γ0 (only dependent on Ω) is a positive constant, and C0 (only dependent on Ω) is an embedding constant of H1(Ω)L4(Ω) ( denotes the continuous embedding).

    Furthermore, we define the Stokes operator A:D(A)˜H such that A=P, where D(A)=H2(Ω)dX0 and P is an L2-projection from L2(Ω)d to ˜H={vL2(Ω)d;v=0,vn|Ω=0}.

    Next, we recall a discrete version of Gronwall's inequality established in [14].

    Lemma 2.1 Let C0,an,bn and dn, for integers n0, be the non-negative numbers such that

    am+τmn=1bnτm1n=1dnan+C0   m1. (2.1)

    Then

    am+τmn=1bnC0exp(τm1n=1dn)   m1. (2.2)

    Throughout this paper, we make the following assumptions on the prescribed date for the problem (1.1), which satisfies the regularity of the data needed for our main results.

    Assumption (A1)(see[27,28]): Assume that the boundary of Ω is smooth so that the unique solution (v,q) of the steady Stokes problem

    Δv+q=f,v=0inΩ,v|Ω=0,

    for prescribed fL2(Ω)3 satisfies

    v2+q1cf0,

    and Maxwell's equations

    curlcurlB=g,B=0inΩ,n×curlB=0,BnonΩ,

    for the prescribed gL2(Ω)3 admits a unique solution BW0 which satisfies

    B2cg0.

    Assumption (A2): The initial data u0,B0,f satisfy

    Au00+B02+sup0tTf0+ft0C. (2.3)

    Assumption (A3): We assume that the MHD problem admits a unique local strong solution (u,p,B) on (0,T) such that

    sup0tTAut0+Bt2+p1+ut0+Bt0+T0(ut21+Bt21+pt21+utt20+Btt20)dtC. (2.4)

    Let 0=t0<t1<<tN=T be a uniform partition of the time interval [0,T] with tn=nτ, and N being a positive integer. Set un=(,tn), For any sequence of functions {fn}Nn=0, we define

    Dτfn=fnfn1τ.

    At the initial time step, we start with U0=˜U0=u0,B0=˜B0=b0 and an arbitrary P0MH1(Ω). Then, we can give

    P0=f0u0u0+νΔu0μB0×curlB0.

    Then, we give the algorithm as follows

    Algorithm 2.1 (Time-discrete modified characteristics projection finite element method)

    Step 1. Solve Bn+1, for instance

    Bn+1ˆBnτ+curlcurlBn+1(Bn)Un=gn+1, (2.5)
    Bn+1n=0,n×curlBn+1=0, on Ω, (2.6)

    where

    ˆln={ln(ˆx),   ˆx=xUnτΩ,0,   othrewise.

    Step 2. Find ˜Un+1, such that

    ˜Un+1ˆUnτνΔ˜Un+1+Pn+μBn×curlBn+1=fn+1, (2.7)
    ˜Un+1=0, on Ω. (2.8)

    Step 3. Calculate (Un+1,Pn+1), for example

    Un+1˜Un+1τ+(Pn+1Pn)=0, (2.9)
    Un+1=0, (2.10)
    Un+1=0 on Ω,. (2.11)

    The corresponding weak form is

    Step 1. Solve Bn+1, for instance

    (Bn+1ˆBnτ,ψ)+(curlBn+1,curlψ)((Bn)Un,ψ)=(gn+1,ψ),   ψW. (2.12)

    Step 2. Find ˜Un+1, such that

    (˜Un+1ˆUnτ,v)+ν(˜Un+1,v)+(Pn,v)+μ(Bnv,Bn+1)μ(vBn,Bn)=(fn+1,v),   vX, (2.13)

    Step 3. Calculate (Un+1,Pn+1), for example

    (Un+1˜Un+1τ,v)+((Pn+1Pn),v)=0,   vX, (2.14)
    (Un+1,ϕ)=0,   ϕM. (2.15)

    Remark 2.1 As we all know (Bn×curlBn+1,v)=(Bnv,Bn+1)(vBn,Bn+1). In order to improve our algorithm, we change it as (Bnv,Bn+1)(vBn,Bn).

    We denote Th be a regular and quasi-uniform partition of the domain Ω into the triangles for d=2 or tetrahedra for d=3 with diameters by a real positive parameter h(h0). The finite element pair (Xh,Mh,Wh) is constructed based on Th.

    Let Wh0n=Xh×Wh. There exits a constant β0>0 (only dependent on Ω) such that

    sup0(v,ψ)Wh0nd(v,q)v0β0q0,  qMh. (2.16)

    There exits a mapping rhL(H2(Ω)2V,Xh) satisfying

    ((vrhv),q)=0,(vrhv)0Chv2,  vH2(Ω)2V,  qMh,

    and a L2-projection operator ρh:MMh satisfying

    qρhq0Chq1,  qH1(Ω)M,

    and a mapping RhL(H2(Ω)2Vn,Wh) satisfying

    (×RhΦ,×Ψ)+(RhΦ,Ψ)=(×Φ,×Ψ)+(Φ,Ψ)=(×Φ,×Ψ),  ΨWh,ΦRhΦ0+hΦRhΦ1Ch2Φ2,  H2(Ω)Vh.

    Here, we define the discrete Stokes operator Ah=Phh and h (see [15] and the references therein) defined by

    (huh,vh)=(uh,vh),uh,vhVh.

    Meanwhile, we define the discrete operator A1hBh=R1h(h××Bh+hBh)Wh (see [15]) as follows

    (A1hBh,Ψ)=(A121hBh,A121hBh)=(×Bh,×Ψ)+(Bh,Ψ),Bh,ΨWh.

    We define the Stokes projection (Rh(u,p),Qh(u,P)):(X,M)(Xh,Mh) by

    ν(Rh(u,p)u,vh)(Qh(u,p)p,vh)=0, (2.17)
    ((Rh(u,p)u),ϕh)=0. (2.18)

    By classical FEM theory ([25,29,30]), we have the following results.

    Lemma 2.2 Assume that uH10(Ω)dHr+1(Ω)d and pL20(Ω)Hr(Ω). Then,

    Rh(u,p)u0+h((Rh(u,p)u)0+Qh(u,p)p0)Chr+1(ur+1+pr), (2.19)

    and

    Rh(u,p)LC(u2+p1). (2.20)

    Lemma 2.3 If (u,p)W2,k(Ω)d×W1,k(Ω) for k>d,

    Rh(u,p)LC(uW1,+pL). (2.21)

    Algorithm 2.2 (The fully discrete modified characteristics projection finite element method)

    Step 1. Solve Bn+1h, for instance

    (Bn+1hˆBnhτ,ψh)+(curlBn+1h,curlψh)((Bnh)unh,ψh)=(gn+1,ψh),   ψhWh. (2.22)

    where

    ˆlnh={lnh(ˆx),   ˆx=xunhτΩ,0,   othrewise.

    Step 2. Find ˜Un+1, such that

    (˜un+1hˆunhτ,vh)+ν(˜un+1h,vh)+(pnh,vh)+μ(Bnhvh,Bn+1h)μ(vhBnh,Bnh)=(fn+1,vh),   vhXh. (2.23)

    Step 3. Calculate (un+1h,pn+1h), for example

    (un+1h˜un+1hτ,vh)+((pn+1hpnh),vh)=0,   vhXh, (2.24)
    (un+1h,ϕh)=0,   ϕhMh. (2.25)

    Here, we make use of the following notation:

    en=unUn,˜e=un˜Un,n=0,1,...,N,εn=bnBn,ˆεn=bnˆBn,n=0,1,...,N,ξn=pnPn,n=0,1,...,N.

    Lemma 3.1 (see[31]). Assume that g1,g2,ρ are three functions defined in Ω and ρ|Ω=0. If

    τ(g1W1,+g2W1,)12,

    then

    ρ(xg1(x)τ)ρ(xg2(x)τ)LqCτρW1,q2g1g2Lq1,ρ(xg1(x)τ)ρ(xg2(x)τ)1Cτρ0g1g2W1,,

    where 1/q1+1/q2=1/q,1<q<.

    The following theorem gives insight on the error estimate between the time-discrete solution and the solution of the time-dependent MHD system (1.1).

    Theorem 3.1 Suppose that assumption A2–A3 are satisfied, there exists a positive C>0 such that

    max0nm(em+120+μεm+120)+mn=0(en+1en20+μεn+1εn20),+τmn=0(ν˜en+120+μcurlεn+120)Cτ2,max0nm(νem+121+μcurlεm+120)+τmn=1(νen+1en21+μcurlεn+1curlεn20)+τmn=1(Dτen+120+μDτεn+120)Cτ2,τmn=0(en+122+˜en+122+εn+122)Cτ2,τmn=0(DτUn+122+Dτ˜Un+122+DτPn+121+DτBn+122+Un+12W2,d+Bn+12W2,d)+max0nm(Un+12+˜Un+12+Pn+11+Bn+12)C.

    Proof. Firstly, we prove the following inequality

    em2+τ3/4UmW2,d1, (3.1)
    εm2+τ3/4BmW2,d1, (3.2)

    by mathematical induction for m=0,1,...,N.

    Since U0=u0,B0=b0 the above inequality hold for m=0.

    We assume that (3.1) and (3.2) hold for mn for some integer n0. By Sobolev embedding theory, we get

    UmLumL+Cem2C, (3.3)
    BmLbmL+Cεm2C, (3.4)

    and

    τUmW1,1/4, (3.5)
    τBmW1,1/4, (3.6)

    when ττ1 for some positive constant τ1.

    To prove (3.1) and (3.2) for m=n+1, we rewrite (1.1) by

    Dτun+1νΔun+1+pn+1=fn+1unˉunτ+Rn+1tr1μbn+1×curlbn+1, (3.7)
    un+1=0, (3.8)
    Dτbn+1+curlcurlbn+1=gn+1bnˉbnτ+Rn+1tr2+bn+1un+1. (3.9)

    where un+1=U(tn+1),bn+1=B(tn+1), and

    ˉln={ln(ˉx),ˉx=xunτΩ,0,otherwise.

    and

    Rn+1tr1=un+1ˉunτun+1t(un+1)un+1,Rn+1tr2=bn+1ˉbnτbn+1t(un+1)bn+1,

    define the truncation error. We can refer to ([31]), there holds that

    τNm=1Rmtr120Cτ2,τNm=1Rmtr220Cτ2. (3.10)

    The corresponding weak form is

    (Dτun+1,v)+ν(un+1,v)(pn+1,v)=(fn+1,v)(unˉunτ,v)+(Rn+1tr1,v)μ(bn+1×curlbn+1,v), (3.11)
    (un+1,ϕ)=0, (3.12)
    (Dτbn+1,ψ)+(curlbn+1,curlψ)=(gn+1,ψ)(bnˉbnτ,ψ)+(Rn+1tr2,ψ)+(bn+1(un+1un),ψ)+(bn+1un,ψ). (3.13)

    Combining (2.13) and (2.14), we obtain

    (Un+1ˆUnτ,v)+ν(˜Un+1,v)(Pn+1,v)+μ(Bnv,Bn+1)μ(vBn,Bn)=(fn+1,v). (3.14)

    Then, we rewrite (3.14) and (2.12) as

    (DτUn+1,v)+ν(˜Un+1,v)(Pn+1,v)=(fn+1,v)(UnˆUnτ,v)μ(Bnv,Bn+1)+μ(vBn,Bn), (3.15)

    and

    (DτBn+1,ψ)+(curlBn+1,curlψ)=(gn+1,ψ)+(BnUn,ψ)(BnˆBnτ,ψ). (3.16)

    Subtracting (3.15) and (3.16) from (3.11) and (3.13), respectively, leads to

    (Dτen+1,v)+ν(˜en+1,v)(ξn+1,v)=(unˉun(UnˆUn)τ,v)+(Rn+1tr1,v)μ((bn+1bn)vv(bn+1bn),bn+1)+μ(vbn,bn+1bn)+μ(vεn,bn)+μ(vBn,εn)μ(εnv,bn+1)μ(Bnv,εn+1), (3.17)

    and

    (Dτεn+1,ψ)+(curlεn+1,curlψ)=(bnˉbn(BnˆBn)τ,ψ)+(Rn+1tr2,ψ)+(bn+1(un+1un)+(bn+1bn)un+εnun+Bnen,ψ). (3.18)

    Testing (2.9) by τen+1, since en+1=0, we arrive at

    (en+1,˜en+1)=(en+1,en+1). (3.19)

    Testing (2.9) by τen, we get

    (en,˜en+1)=(en,en+1). (3.20)

    Subtracting (3.20) from (3.19), we can obtain

    (en+1en,˜en+1)=(en+1en,en+1). (3.21)

    Let v=2τ˜en+1 in (3.17), we obtain

    en+120en20+en+1en20+2ντ˜en+1202τ(ξn+1,˜en+1)=2((en(ˉunˆUn),˜en+1)+2τ(Rn+1tr1,˜en+1)2μτ((bn+1bn)˜en+1˜en+1(bn+1bn),bn)+2μτ(˜en+1bn,bn+1bn)+2μτ(˜en+1εn,bn)+2μτ(˜en+1Bn,εn)2μτ(εn˜en+1,bn+1)2μτ(Bn˜en+1,εn+1). (3.22)

    Let ψ=2μτεn+1 in (3.18), we get

    μεn+120μεn20+μεn+1εn20+2μτcurlεn+120=2μ(εn(ˉbnˆBn),εn+1)+2μτ(Rn+1tr2,εn+1)+2μτ(bn+1(un+1un)+(bn+1bn)un+εnun+Bnen,εn+1). (3.23)

    Taking sum of (3.22) and (3.23) yields

    en+120en20+μεn+120μεn20+en+1en20+μεn+1εn20+2ντ˜en+120+2μτcurlεn+1202τ(ξn+1,˜en+1)=2(enˆen,˜en+1)2(ˆunˉun,˜en+1)+2τ(Rn+1tr1,˜en+1)2μτ((bn+1bn)˜en+1,bn)+2μτ(˜en+1(bn+1bn),bn)+2μτ(˜en+1bn,bn+1bn)+2μτ(˜en+1εn,bn)+2μτ(˜en+1Bn,εn)2μτ(εn˜en+1,bn+1)2μτ(Bn˜en+1,εn+1)2μ(εnˆεn,εn+1)2μ(ˆbnˉbn,εn+1)+2μτ(Rn+1tr2,εn+1)+2μτ(bn+1(un+1un),εn+1)+2μτ((bn+1bn)un,εn+1)+2μτ(εnunεn+1)+2μτ(Bnen,εn+1). (3.24)

    From (2.9), we have

    ˜en+1=en+1+τ(ξn+1ξn)τ(pn+1pn).

    Then, we can obtain

    2τ(ξn+1,˜en+1)=2τ(ξn+1,en+1)+2τ2(ξn+1,(ξn+1ξn))4τ2(ξn+1,(pn+1pn))=τ2(ξn+120ξn20+ξn+1ξn20)2τ2(ξn+1,(pn+1pn))τ2(ξn+120ξn20+(ξn+1ξn)20)Cτ2ξn+10(pn+1pn)0=τ2(ξn+120ξn20+(ξn+1ξn)20)τ(Cτξn+10(pn+1pn)0)τ2(ξn+120ξn20+(ξn+1ξn)20)τ(Cτ2ξn+120+14(pn+1pn)20)=τ2(ξn+120ξn20+(ξn+1ξn)20)Cτ3ξn+120+14(pn+1pn)20).

    On the other hand, using Lemma 3.1, we deduce that

    2|(enˆen,˜en+1)|2enˆen1˜en+10Cτen0UnW1,˜en+10ντ12˜en+120+Cτen20,2|(ˆunˉun,˜en+1)|CˆunˉunL6/5˜en+1L6CτunL3en0˜en+10ντ12˜en+120+Cτen20,2τ|(Rn+1tr1,˜en+1)|Rn+1tr10˜en+10ντ12˜en+120+CτRn+1tr120,2μτ|((bn+1bn)˜en+1,bn+1)|Cτbn+1bn0˜en+10bn+12ντ12˜en+120+Cτ3bn+1t20bn+122,2μτ|(˜en+1(bn+1bn),bn+1)|Cτ˜en+1L6bn+1L3bn+1bn0ντ12˜en+120+Cτ3bn+1t20bn+122,2μτ|(˜en+1bn,bn+1bn)|Cτ˜en+1L6bn+1L3bn+1bn0ντ12˜en+120+Cτ3bn+1t20bn+122,2μτ|(˜en+1εn,bn)|Cτ˜en+10bn2εn0ντ12˜en+120+Cτεn20bn22,2μτ|(˜en+1Bn,εn)|Cτ˜en+10BnLεn0ντ12˜en+120+Cτεn20Bn2L,2μτ|(εn˜en+1,bn+1)|Cτεn0˜en+10bn+1Lντ12˜en+120+Cτεn20bn+122,2μτ|(Bn˜en+1,εn+1)|CτBnL˜en+10εn+10ντ12˜en+120+Cτεn20,2μτ|(εnˆεn,εn+1)|Cτεn0BnW1,εn+10μτ8curlεn+120+Cτεn20,2μτ|(ˆbnˉbn,εn+1)|CτˆbnˉbnL6/5εn+1L6CτbnL3εn0εn+10μτ8curlεn+120+Cτεn20,2μτ|(Rn+1tr2,εn+1)|CτRn+1tr20εn+10μτ8curlεn+120+CτRn+1tr220,2μτ|(bn+1(un+1un),εn+1)|Cτbn+12un+1un0curlεn+10μτ8curlεn+120+Cτun+1un20,2μτ|((bn+1bn)un,εn+1)|Cτbn+1bn0unLεn+10μτ8curlεn+120+Cτ3bn+1t20un2L,2μτ|(εnun,εn+1)|Cτεn0unL3curlεn+10μτ8curlεn+120+Cτεn20un2L3,2μτ|(Bnen,εn+1)|CτBnLen0curlεn+10μτ8curlεn+120+Cτen20.

    Substituting these above inequality into (3.24), we obtain

    en+120en20+μεn+120μεn20+τ2ξn+120τ2ξn20+en+1en20+μεn+1εn20+τ2(ξn+1ξn)20+ντ˜en+120+μτcurlεn+120Cτ(en20+εn20+εn+120)+Cτ(Rn+1tr120+Rn+1tr220)+Cτ3(ξn+120+ξn20)+Cτ3. (3.25)

    Taking sum of the (3.25) for n from 0 to mN and using discrete Gronwall's inequality Lemma 2.1, we get

    max0nm(em+120+μεm+120+τ2ξm+120)+mn=0(en+1en20+μεn+1εn20)+τmn=0(ν˜en+120+μcurlεn+120)Cτ2. (3.26)

    Owing to (3.21), we have

    ((en+1en),˜en+1)=((en+1en),en+1). (3.27)

    To acquire the H1 estimates, we take v=2τDτen+1 in (3.17) to get

    ν(en+121en21+en+1en21)+2τDτen+120=2((en(ˉunˆUn),Dτen+1)+2τ(Rn+1tr1,Dτen+1)2μτ((bn+1bn)Dτen+1Dτen+1(bn+1bn),bn+1)+2μτ(Dτen+1bn,bn+1bn)+2μτ(Dτen+1εn,bn)+2μτ(Dτen+1Bn,εn)2μτ(εnDτen+1,bn+1)2μτ(BnDτen+1,εn+1). (3.28)

    Then, we take ψ=2μτDτεn+1 in (3.18) to get

    μ(curlεn+120curlεn20+curlεn+1curlεn20)+2μτDτεn+120=2μ(εn(ˉbnˆBn),Dτεn+1)+2μτ(Rn+1tr2,Dτεn+1)+2μτ(bn+1(un+1un)+(bn+1bn)un+εnun+Bnen,Dτεn+1). (3.29)

    Combining (3.28) and (3.29), we obtain

    ν(en+121en21+en+1en21)+μ(curlεn+120curlεn20+curlεn+1curlεn20)+2τDτen+120+2ντDτεn+120=2(enˆen,Dτen+1)2(ˆunˉun,Dτen+1)+2τ(Rn+1tr1,Dτen+1)2μτ((bn+1bn)Dτen+1,bn+1)+2μτ(Dτen+1(bn+1bn),bn+1)+2μτ(Dτen+1bn,bn+1bn)+2μτ(Dτen+1εn,bn)+2μτ(Dτen+1Bn,εn)2μτ(εnDτen+1,bn+1)2μτ(BnDτen+1,εn+1)2μ(εnˆεn,εn+1)2μ(ˆbnˉbn,Dτεn+1)+2μτ(Rn+1tr2,Dτεn+1)+2μτ(bn+1(un+1un),Dτεn+1)+2μτ((bn+1bn)un,Dτεn+1)+2μτ(εnun,Dτεn+1)+2μτ(Bnen,Dτεn+1). (3.30)

    Similarly, by Lemma 3.1, we have

    2|(enˆen,Dτen+1)|2enW1,2UnLDτen+10τ16Dτen+120+Cτen21,2|(ˆunˉun,Dτen+1)|CunW1,4enL4Dτen+10CunW1,4en1Dτen+10τ16Dτen+120+Cτen21,2τ|(Rn+1tr1,Dτen+1)|Rn+1tr10Dτen+10τ16Dτen+120+CτRn+1tr120,2μτ|((bn+1bn)Dτen+1,bn+1)|Cτ(bn+1bn)0Dτen+10bn+12τ16Dτen+120+Cτ3bn+1t20bn+122,2μτ|(Dτen+1(bn+1bn),bn+1)|CτDτen+10(bn+1bn)0bn+12τ16Dτen+120+Cτ3bn+1t20bn+122,2μτ|(Dτen+1bn,bn+1bn)|CτDτen+10bnLbn+1bn0τ16Dτen+120+Cτ3bn+1t20bn2L,2μτ|(Dτen+1εn,bn)|CτDτen+10εn0bnLτ16Dτen+120+Cτcurlεn20bn2L,2μτ|(Dτen+1Bn,εn)|CτDτen+10BnL3εnL6τ16Dτen+120+Cτεn20Bn2L3,τ16Dτen+120+Cτcurlεn20Bn2L3,2μτ|(εnDτen+1,bn+1)|Cτεn0Dτen+10bn+12τ16Dτen+120+Cτcurlεn20bn+122,2μτ|(BnDτen+1,εn+1)|CτBnLεn+10Dτen+10τ16Dτen+120+Cτcurlεn+120Bn2L,2μτ|(εnˆεn,Dτεn+1)|CτεnW1,2BnLDτεn+10μτ8Dτεn+120+Cτcurlεn20,2μτ|(ˆbnˉbn,Dτεn+1)|CτbnL4εnL4Dτεn+10μτ8Dτεn+120+Cτcurlεn20,2μτ|(Rn+1tr2,Dτεn+1)|CτRn+1tr20Dτεn+10μτ8Dτεn+120+CτRn+1tr220,2μτ|(bn+1(un+1un),Dτεn+1)|Cτbn+12(un+1un)0Dτεn+10μτ8Dτεn+120+Cτ(un+1un)20,2μτ|((bn+1bn)un,Dτεn+1)|Cτbn+1bn0unLDτεn+10μτ8Dτεn+120+Cτ3bn+1t20un2L,2μτ|(εnun,Dτεn+1)|Cτεn0unLDτεn+10μτ8Dτεn+120+Cτcurlεn20un2L,2μτ|(Bnen,Dτεn+1)|CτBnLen0Dτεn+10μτ8Dτεn+120+Cτen20Bn2L.

    Substituting these above inequality into (3.30), we obtain

    ν(en+121en21+en+1en21)+μ(curlεn+120curlεn20+curlεn+1curlεn20)+τDτen+120+μτDτεn+120Cτ(en21+curlεn+120+curlεn20)+Cτ(Rn+1tr120+Rn+1tr220)+Cτ3. (3.31)

    Taking sum of the (3.31) for n from 0 to mN and using discrete Gronwall's inequality Lemma 2.1, we arrive at

    max0nm(νem+121+μcurlεm+120)+τmn=1(νen+1en21+μcurlεn+1curlεn20)+τmn=1(Dτen+120+μDτεn+120)Cτ2. (3.32)

    Furthermore, the standard result for (3.17) and (3.18) with p=2, respectively, leads to

    ˜en+12CDτen+10+Cτenˆen0+CτˆUnˉUn0+CRn+1tr10+Cbn+1bn0bn+1L+Cbn+1bn0bn+1L+CbnLbn+1bn0+Cεn0bnL+CBnW1,εn0+Cεn0bn+1L+CBnLcurlεn+10+ξn+11CDτen+10+CRn+1tr10+ξn+11+Cτ, (3.33)
    εn+12CDτεn+10+Cτεnˆεn0+Cτˆbnˉbn0+CRn+1tr20+Cbn+1Lun+1un0+Cbn+1bn0unL+Cεn0unL+CBnLen0CDτen+10+CRn+1tr0+Cτ. (3.34)

    Thanks to (3.19), we can obtain

    (Δen+1,Δen+1)=(Δ˜en+1,Δen+1),en+12˜en+12.

    Then, we have

    τmn=0(en+122+˜en+122+εn+122)Cτ2. (3.35)

    From (3.35), we can arrive at

    max0nmUn+12max0nm(un+12+en+12)C, (3.36)
    max0nm˜Un+12max0nm(un+12+˜en+12)C, (3.37)
    max0nmPn+11max0nm(pn+11+ξn+11)C, (3.38)
    max0nmBn+12max0nm(bn+12+εn+12)C, (3.39)
    τmn=0DτUn+1222τmn=0(Dτun+122+Dτen+122)C, (3.40)
    τmn=0Dτ˜Un+1222τmn=0(Dτun+122+Dτ˜en+122)C, (3.41)
    τmn=0DτPn+1212τmn=0(Dτpn+121+Dτξn+121)C, (3.42)
    τmn=0DτBn+1222τmn=0(Dτbn+122+Dτεn+122)C. (3.43)

    From (2.7) and (2.5), we have

    DτUn+1νΔ˜Un+1+Pn+1=fn+1UnˆUnτμBn×curlBn+1, (3.44)

    and

    DτBn+1+curlcurlBn+1=gn+1+(Bn)UnBn+1ˆBnτ. (3.45)

    By (2.9), we have

    ΔUn+1=×טUn+1. (3.46)

    Considering the identity 2˜Un+1=˜Un+1×טUn+1 and (3.46), we can obtain

    Δ˜Un+1=˜Un+1ΔUn+1. (3.47)

    Now, the standard result for the Stokes system (3.17) and (3.18) with p=d>2, respectively, leads to

    Un+1W2,d+Pn+1˜Un+1W1,dCDτUn+1Ld+CτˆUnUnLd+Cfn+1Ld+CBn×curlBn+1LdC(DτUn+1Ld+UnLdUnL+fn+1Ld+BnLdcurlBn+1Ld), (3.48)
    Bn+1W2,dCDτBn+1Ld+Cgn+1Ld+CτˆBnBnLd+CBnUnLdC(DτBn+1Ld+gn+1Ld+BnLdBnL+BnLdUnLd). (3.49)

    By (3.48) and (3.49), we get

    τmn=0(Un+12W2,d+Bn+12W2,d)C. (3.50)

    Thus, we can obtain

    en+12+τ3/4Un+1W2,d1, (3.51)
    εn+12+τ3/4Bn+1W2,d1. (3.52)

    The proof is complete.

    For the sake of simplicity, we denote Rnh=Rh(Un,Pn),Qnh=Qh(Un,Pn),R0h=R0hBn,n=1,2,...,N, and enh=unhRnh,˜enh=˜unhRnh,εnh=BnhRn0h,n=1,2,...,N, where R0h:L2(Ω)3Wh is the L2-orthogonal projection, where

    R0hww0+h(R0hww)0Chr+1wr+1, (3.53)
    R0hwLCw2, (3.54)
    R0hwLCwW1,. (3.55)

    Theorem 3.2 Suppose that assumption A2–A3 are satisfied, there exists a positive C>0 such that

    um+1h20+μBm+1h20+mn=0(un+1hunh20+μBn+1hBnh20)+τmn=0(ν˜un+1h20+μcurlBn+1h20)C,um+1h21+μcurlBm+1h20+mn=0(un+1hunh21+μcurlBn+1hcurlBnh20)+τmn=0(νAun+1h20+μcurlBn+1h20)C,em+1h20+μεm+1h20+mn=0(en+1henh20+μεn+1hεnh20)+τmn=0(ν˜en+1h20+μcurlεn+1h20)Ch4,un+1h2L+μBn+1h2L+τmn=0(un+1h2W1,+μBn+1h2W1,)C.

    Proof. We prove following estimates

    emh20+μεmh20+τnm=0(νemh20+μcurlεmh20)Ch4, (3.56)

    by mathematical induction for m=0,1,...,N.

    When m=0,u0h=u0,B0h=B0. It is obvious (3.56) holds at the initial time step. We assume that (3.56) holds for 0mn for some integer n0. By (2.20), (2.21), (3.54), (3.55), Theorem 3.1 and inverse inequality, we infer

    τumhW1,(emhW1,+RmhW1,)C(hd2emh0+UmW1,+PmL)14, (3.57)
    τBmhW1,(εmhW1,+Rm0hW1,)C(hd2εmh0+BmW1,)14, (3.58)

    and

    τumhL(emhL+RmhL)C(hd2emh0+Um2+Pm1)C, (3.59)
    τBmhL(εmhL+Rm0hL)C(hd2εmh0+Bm2)C. (3.60)

    Combining (2.23) and (2.24), we obtain

    (un+1hˆunhτ,vh)+ν(˜un+1h,vh)(pn+1h,vh)+μ(Bnhvh,Bn+1h)μ(vhBnh,Bnh)=(fn+1,vh). (3.61)

    When m=n+1, letting vh=2τ˜un+1h in (3.61), we can obtain

    2(un+1hunh+unhˆunh,˜un+1h)+2ντ(˜un+1h,˜un+1h)2τ(pn+1h,˜un+1h)+2μτ(Bnh˜un+1h,Bn+1h)2μτ(˜un+1hBnh,Bnh)=2τ(fn+1,˜un+1h). (3.62)

    Letting vh=τun+1h in (2.24), since un+1h=0, it follows that

    (un+1h,˜un+1h)=(un+1h,un+1h). (3.63)

    Taking vh=τunh in (2.24), we get

    (unh,˜un+1h)=(unh,un+1h). (3.64)

    Subtracting (3.64) from (3.63), we can obtain

    (un+1hunh,˜un+1h)=(un+1hunh,un+1h). (3.65)

    Using 2(ab,a)=a20b20+ab20, leads to

    un+1h20unh20+un+1hunh20+2ντ˜un+1h202τ(pn+1h,˜un+1h)=2τ(fn+1,˜un+1h)+2(ˆunhunh,˜un+1h)2μτ(Bnh˜un+1h,Bn+1h). (3.66)

    Taking vh=pn+1h in (2.24), we deduce

    2τ(pn+1h,˜un+1)=2τ2pn+1h202τ2pnh20+2τ2pn+1hpnh20. (3.67)

    When m=n+1, letting ψh=2μτBn+1h in (2.22), we can obtain

    2μ(Bn+1hBnh+BnhˆBnh,Bn+1h)+2μτ(curlBn+1h,curlBn+1h)2μτ(Bnhunh,Bn+1h)=2μτ(gn+1,Bn+1h). (3.68)

    Using 2(ab,a)=a20b20+ab20, we have

    μBn+1h20μBnh20+μBn+1hBnh20+2μτcurlBn+1h20=2μτ(gn+1,Bn+1h)+2μτ(Bnhunh,Bn+1h)+2μ(ˆBnhBnh,Bn+1h). (3.69)

    Taking sum of (3.66) and (3.69) yields

    un+1h20unh20+μBn+1h20μBnh20+2τ2pn+1h202τ2pnh20+un+1hunh20+μBn+1hBnh20+2τ2pn+1hpnh20+2ντ˜un+1h20+2μτcurlBn+1h20=2τ(fn+1,˜un+1h)+2μτ(gn+1,Bn+1h)+2(ˆunhunh,˜un+1h)2μτ(Bn˜un+1h,Bn+1h)+2μ(ˆBnhBnh,Bn+1h)+2μτ(Bnhunh,Bn+1h). (3.70)

    Then, we can obtain

    2τ|(fn+1,˜un+1h)|Cτfn+10˜un+1h0ντ4˜un+1h20+Cτfn+120,2μτ|(gn+1,Bn+1h)|Cτgn+10Bn+1h0μτ4curlBn+1h20+Cτgn+120,2|(ˆunhunh,˜un+1h)|Cτˆunhunh6/5˜un+1hL6CτunhL3unh0˜un+1h0ντ4˜un+1h20+Cτunh20,2μτ|(Bn˜un+1h,Bn+1h)|CτBnL˜un+1h0Bn+1h0ντ4˜un+1h20+CτBn+1h20,2μ|(ˆBnhBnh,Bn+1h)|CτˆBnhBnh1Bn+1h1CτBnh0BnhW1,curlBn+1h0μτ4curlBn+1h20+CτBnh20Bnh2W1,,2μτ|(Bnunh,Bn+1h)|CτBnLunh0curlBn+1h0μτ4curlBn+1h20+Cτunh20.

    Substituting these above estimates into (3.70), we obtain

    un+1h20unh20+μBn+1h20μBnh20+τ2pn+1h20τ2pnh20+un+1hunh20+μBn+1hBnh20+2τ2pn+1hpnh20+ντ˜un+1h20+μτcurlBn+1h20Cτfn+120+Cτgn+120+Cτunh20+CτBnh20. (3.71)

    Taking sum of the (3.71) for n from 0 to mN and using discrete Gronwall's inequality Lemma 2.1, we get

    um+1h20+μBm+1h20+τ2pm+1h20+mn=0(un+1hunh20+μBn+1hBnh20)+τmn=0(ν˜un+1h20+μcurlBn+1h20)C.

    Letting vh=2τAun+1h in (3.61), we obtain

    2(un+1hunh+unhˆunh,Aun+1h)+2ντ(A˜un+1h,Aun+1h)+2μτ(BnhAun+1h,Bn+1h)2μτ(Aun+1hBnh,Bnh)=2τ(fn+1,Aun+1h). (3.72)

    From (3.63), we arrive at

    (A˜un+1h,Aun+1h)=(Aun+1h,Aun+1h),(A˜un+1h,Aun+1h)=Aun+1h20.

    Using 2(ab,a)=a20b20+ab20, leads to

    un+1h21unh21+un+1hunh21+2ντAun+1h20=2τ(fn+1,Aun+1h)+2(ˆunhunh,Aun+1h)2μτ(BnhAun+1h,Bn+1h). (3.73)

    Letting ψh=2μτcurlcurlBn+1h in (2.22), we can obtain

    μcurlBn+1h20μcurlBnh20+μcurlBn+1hcurlBnh20+2μτcurlcurlBn+1h20=2μτ(gn+1,curlcurlBn+1h)+2μτ(Bnhunh,curlcurlBn+1h)+2μ(ˆBnhBnh,curlcurlBn+1h). (3.74)

    Taking sum of (3.73) and (3.74) yields

    un+1h21unh21+μcurlBn+1h20μcurlBnh20+un+1hunh21+μcurlBn+1hcurlBnh20+2ντAun+1h20+2μτcurlcurlBn+1h20=2τ(fn+1,Aun+1h)+2(ˆunhunh,Aun+1h)+2τ(gn+1,curlcurlBn+1h)2μτ(BnhAun+1h,Bn+1h)+2μτ(Bnhunh,curlcurlBn+1h)+2μ(ˆBnhBnh,curlcurlBn+1h). (3.75)

    Then, we get

    2τ|(fn+1,Aun+1h)|Cτfn+10Aun+1h0ντ6Aun+1h20+Cτfn+120,2μτ|(gn+1,curlcurlBn+1h)|Cτgn+10curlcurlBn+1h0μτ4curlcurlBn+1h20+Cτgn+120,2|(ˆunhunh,Aun+1h)|CτunhW1,4unhL4Aun+1h0ντ6Aun+1h20+Cτunh21,2μτ|(BnhAun+1h,Bn+1h)|CτBnhW1,Aun+1h0Bn+1h0ντ6Aun+1h20+CτcurlBn+1h20,2μτ|(Bnhunh,curlcurlBn+1h)|CτBnhLunh0curlcurlBn+1h0μτ4curlcurlBn+1h20+Cτunh20,2μ|(ˆBnhBnh,curlcurlBn+1h)|CτBnhW1,4BnhL4curlcurlBn+1h0μτ4curlcurlBn+1h20+CτBnh2W1,4Bnh2L4μτ4curlcurlBn+1h20+CτcurlBnh20.

    Substituting these above estimates into (3.75), we obtain

    un+1h21unh21+μcurlBn+1h20μcurlBnh20+un+1hunh21+μcurlBn+1hcurlBnh20+ντAun+1h20+μτcurlcurlBn+1h20Cτfn+120+Cτgn+120+Cτunh21+Cτ(curlBn+1h20+curlBnh20). (3.76)

    Taking sum of the (3.76) for n from 0 to mN and using discrete Gronwall's inequality Lemma 2.1, we get

    um+1h21+μcurlBm+1h20+mn=0(un+1hunh21+μcurlBn+1hcurlBnh20)+τmn=0(νAun+1h20+μcurlcurlBn+1h20)C.

    Then, we rewrite (3.61) and (2.22) as

    (Dτun+1h,vh)+ν(˜un+1h,vh)(Pn+1h,vh)=(fn+1,vh)(unhˆunhτ,vh)μ(Bnhvh,Bn+1h)+μ(vhBnh,Bnh), (3.77)

    and

    (DτBn+1h,ψh)+(curlBn+1h,curlψh)=(gn+1,ψh)+(Bnhunh,ψh)(BnhˆBnhτ,ψh). (3.78)

    Subtracting (3.15) and (3.16) from (3.77) and (3.78), respectively, leads to

    (Dτen+1h,vh)+ν(˜en+1h,vh)(pn+1hQn+1h,vh)=(Dτ(Un+1Rn+1h),vh)+(ˆunhunhτ,vh)+(UnˆUnτ,vh)μ(εnhvh,Bn+1h)μ((Rn0hBn)vh,Bn+1h)μ(Bnvh,εn+1h)μ(Bnvh,Rn+10hBn+1)+μ(vhεnh,Bnh)+μ(vhBn,εnh)μ(vhBn,Rn0hBn)+μ(vh(Rn0hBn),Bnh), (3.79)

    and

    (Dτεn+1h,ψh)+(curlεn+1h,curlψnh)=(Dτ(Bn+1Rn+10h),ψh)+(BnhˆBnhτ,ψh)+(ˆBnBnτ,ψh)+(εnhunh,ψh)+((Rn0hBn)unh,ψh)+(Bn(RnhUn),ψh)+(Bnenh,ψh). (3.80)

    Taking vh=2τen+1h in (2.24), we can obtain

    (en+1h,˜en+1h)=(en+1h,en+1h),(Dτen+1h,˜en+1h)=(Dτen+1h,en+1h),(en+1h,˜en+1h)=(en+1h,en+1h)=en+1h20.

    Let θn=UnRnh,ηn=BnRn0h. By taking vh=2τen+1h in (3.79), we deduce

    en+1h20enh20+en+1henh20+2ντen+1h20=2τ(DτUn+1Rh(DτUn+1,DτPn+1),en+1h)+2(ˆenhenh,en+1h)+2(θnˆθn,en+1h)2μτ(εnhen+1h,Bn+1h)2μτ(ηnen+1h,Bn+1h)2μτ(Bnen+1h,εn+1h)2μτ(Bnen+1h,ηn+1)+2μτ(en+1hεnh,Bnh)+2μτ(en+1hBn,εnh)+2μτ(en+1hBn,ηn)+2μτ(en+1hηn,Bnh). (3.81)

    By taking ψh=2μτεn+1h in (3.80), we have

    μεn+1h20μεnh20+μεn+1hεnh20+2μτcurlεn+1h20=2τ(DτBn+1R0hDτBn+1),εn+1h)+2μ(εnhˆεnh,εn+1h)+2μ(ηnˆηn,εn+1h)+2μτ(εnhunh,εn+1h)2μτ(ηnunh,εn+1h)+2μτ(Bn(RnhUn),εn+1h)+2μτ(Bnenh,εn+1h). (3.82)

    Combining (3.81) and (3.82), we obtain

    en+1h20enh20+μεn+1h20μεnh20+en+1henh20+Sεn+1hεnh20+2ντen+1h20+2μτcurlεn+1h20=2τ(DτUn+1Rh(DτUn+1,DτPn+1),en+1h)+2(ˆenhenh,en+1h)+2(θnˆθn,en+1h)2μτ(εnhen+1h,Bn+1h)2μτ(ηnen+1h,Bn+1h)2μτ(Bnen+1h,Bn+1h)2μτ(Bnen+1h,ηn+1)+2μτ(en+1hεnh,Bnh)+2μτ(en+1hBn,εnh)+2μτ(en+1hBn,ηn)+2μτ(en+1hηn,Bnh)+2τ(DτBn+1R0hDτBn+1),εn+1h)+2μ(εnhˆεnh,εn+1h)+2μ(ηnˆηn,εn+1h)+2μτ(εnhunh,εn+1h)2μτ(ηnunh,εn+1h)+2μτ(Bn(RnhUn),εn+1h)2μτ(Bnenh,Bn+1h). (3.83)

    Due to (2.19) and Theorem 3.1, we can get θn0Ch2(Un2+Pn1)Ch2, ηn0Ch2Bn2. By Lemma 3.1, there holds that

    2τ|(DτUn+1Rh(DτUn+1,DτPn+1),en+1h)|Cτh2en+1h0(DτUn+12+DτPn+11)ντ20en+1h20+Cτh4(DτUn+122+DτPn+121),2|(ˆenhenh,en+1h)|Cenhˆenh1en+1h0Cτenh0unhW1,en+1h0ντ20en+1h20+Cτenh20,2|(θnˆθn,en+1h)|Cθnˆθn1en+1h0Cτθn0unhW1,en+1h0ντ20en+1h20+Cτθn20,2μτ|(εnhen+1h,Bn+1h)|Cτεnh0en+1h0Bn+1hLντ20en+1h20+Cτεnh20,2μτ|(ηnen+1h,Bn+1h)|Cτηn0en+1h0Bn+1hLντ20en+1h20+Cτηn20,2μτ|(Bnen+1h,εn+1h)|CτBnLen+1h0εn+1h0ντ20en+1h20+Cτεn+1h20,2μτ|(Bnen+1h,ηn+1)|CτBn2en+1h0ηn+10ντ20en+1h20+Cτηn+120,2μτ|(en+1hεnh,Bnh)|Cτen+1h0εnh0BnhLντ20en+1h20+Cτεnh20,2μτ|(en+1hBn,εn+1h)|Cτen+1h0Bn2εnh0ντ20en+1h20+Cτεnh20,2μτ|(en+1hBn,ηn)|Cτen+1h0Bn2ηn0ντ20en+1h20+Cτηn20,2μτ|(en+1hηn,Bnh)|Cτen+1h0ηn0BnhLντ20en+1h20+Cτηn20,2τ|(DτBn+1R0hDτBn+1,εn+1h)|Cτh2εn+1h0DτBn+12ηn0μτ16curlεn+1h20+Cτh4DτBn+122,2μ|(εnhˆεnh,εn+1h)|Cεnhˆεnh0εn+1h0Cτcurlεnh0BnhLεn+1h0μτ16curlεn+1h20+Cτεn+1h20,2μ|(ηnˆηn,εn+1h)|Cηnˆηn1εn+1h0Cτηn0BnhW1,curlεn+1h0μτ16curlεn+1h20+Cτηn20,2μτ|(εnhunh,εn+1h)|Cτεnh0unhL3εn+1hL6μτ16curlεn+1h20+Cτεn+1h20,2μτ|(ηnunh,εn+1h)|Cτηn0unhL3εn+1hL6μτ16curlεn+1h20+Cτηn20,2μτ|(Bn(RnhUn),εn+1h)|CτBn2θn0εn+1h0μτ16curlεn+1h20+Cτθn20,2μτ|(Bnenh,εn+1h)|CτBnLenh0curlεn+1h0μτ16curlεn+1h20+Cτenh20.

    Substituting these above estimates into (3.83), we obtain

    en+1h20enh20+μεn+1h20μεnh20+en+1henh20+μεn+1hεnh20+ντen+1h20+μτcurlεn+1h20Cτ(enh20+εnh20+εn+1h20)+Cτh4. (3.84)

    Taking the sum of the (3.76) for n from 0 to mN and using the discrete Gronwall's inequality Lemma 2.1, we get

    em+1h20+μεm+1h20+mn=0(en+1henh20+μεn+1hεnh20)+τmn=0(νen+1h20+μcurlεn+1h20)Ch4.

    Then, we can deduce

    max0nmun+1hLmax0nm(Rn+1hL+en+1hL)Cmax0nmUn+12+Cmax0nmPn+11+Ch12max0nmen+1h0C,max0nmBn+1hLmax0nm(Rn+10hL+εn+1hL)Cmax0nmBn+12+Ch12max0nmεn+1h0C,τmn=0un+1h2W1,2τmn=0(Rn+1h2W1,+en+1h2W1,)Cτmn=0(Un+12W1,+Pn+12L)+Cτh2mn=0en+1h20C,τmn=0Bn+1h2W1,2τmn=0(Rn+10h2W1,+εn+1h2W1,)Cτmn=0Bn+12W1,+Cτh2mn=0εn+1h20C.

    Thus, all the results in Theorem 3.2 have been proved.

    In order to show the effect of our method, we present some numerical results. Here, we consider the

    u1=(y+y4)cos(t),u2=(x+x2)cos(t),p=(2.0x1.0)(2.0y1.0)cos(t),B1=(sin(y)+y)cos(t),B2=(sin(x)+x2)cos(t).

    The boundary conditions and the forcing terms are given by the exact solution. The finite element spaces are chosen as P1b-P1-P1b finite element spaces. Here, we choose τ=h2 and h=1/n, n=8,16,24,32,40,48. Different Reynolds number Re and magnetic Reynolds number Rm are chosen to show the effect of our method, the numerical results were shown in Tables 16. Here, we use the software package FreeFEM++ [32] for our program.

    Table 1.  The numerical results for Re=1.0 and Rm=1.0 for different h at T=1.0.
    1/h u1hu0 (u1hu1)0 p1hp0 B1hB0 (B1hB)0 E1hE0 |ΩBdx| |Ωudx|
    8 0.00395939 0.104518 0.0633311 0.00146071 0.0443105 0.0402426 4.95209e-17 1.99024e-14
    16 0.00109704 0.0527177 0.0212408 0.000384435 0.0221793 0.0208612 1.41995e-16 4.15508e-15
    24 0.000473611 0.0347495 0.0118957 0.000175071 0.0148256 0.0141051 5.03032e-17 5.95751e-15
    32 0.000256132 0.0258999 0.010131 9.82405e-05 0.0111161 0.010656 3.35182e-16 7.38832e-15
    40 0.000160786 0.0206654 0.00831346 6.2778e-05 0.00889124 0.00856254 7.97451e-16 5.55468e-15
    48 0.000110796 0.0171977 0.00648101 4.35568e-05 0.00740848 0.00715659 1.00015e-16 4.15077e-15

     | Show Table
    DownLoad: CSV
    Table 2.  Convergence order for Re=1.0 and Rm=1.0 for different h at T=1.0.
    1/h uL2 uH1 PL2 BL2 BH1
    16 1.85166 0.987389 1.57608 1.92586 0.998435
    24 2.07166 1.02792 1.42982 1.93995 0.99343
    32 2.13671 1.0217 0.5582 2.00837 1.00097
    40 2.08665 1.01181 0.88607 2.00685 1.00084
    48 2.04246 1.00748 1.36571 2.00491 1.00066

     | Show Table
    DownLoad: CSV
    Table 3.  The numerical results for Re=100.0 and Rm=100.0 for different h at T=1.0.
    1/h u1hu0 (u1hu1)0 p1hp0 B1hB0 (B1hB)0 E1hE0 |ΩBdx| |Ωudx|
    8 0.0172667 0.279161 0.0354359 0.0127351 0.19731 0.00779509 3.41253e-17 3.17954e-14
    16 0.00233046 0.0697115 0.01409 0.00240034 0.0479792 0.00159673 9.63754e-17 1.3337e-14
    24 0.000831767 0.0398145 0.00892872 0.000977778 0.0269431 0.000711847 8.03834e-18 8.42217e-15
    32 0.000422379 0.0280313 0.00649219 0.000527044 0.0184095 0.000385869 2.36579e-16 5.31439e-15
    40 0.000256947 0.0217553 0.00511173 0.000330314 0.013926 0.000239563 7.44033e-16 2.85026e-15
    48 0.000173343 0.0178297 0.00421705 0.000226765 0.0112143 0.000162916 7.14376e-17 2.41418e-15

     | Show Table
    DownLoad: CSV
    Table 4.  Convergence order for Re=100.0 and Rm=100.0 for different h at T=1.0.
    1/h uL2 uH1 PL2 BL2 BH1
    16 2.88931 2.00163 1.33054 2.4075 2.03998
    24 2.54095 1.38147 1.12512 2.21495 1.42315
    32 2.35555 1.21978 1.10773 2.1482 1.32391
    40 2.22742 1.13589 1.07134 2.09391 1.2508
    48 2.1588 1.09144 1.05529 2.063 1.1878

     | Show Table
    DownLoad: CSV
    Table 5.  The numerical results for Re=1000.0 and Rm=1000.0 for different h at T=1.0.
    1/h u1hu0 (u1hu1)0 p1hp0 B1hB0 (B1hB)0 E1hE0 |ΩBdx| |Ωudx|
    8 0.0397337 1.68458 0.0330041 0.029124 1.17999 0.0213862 1.41814e-16 1.99053e-14
    16 0.00470743 0.368592 0.0137702 0.00403407 0.219747 0.00228304 7.44847e-17 1.26213e-14
    24 0.00158913 0.177215 0.00901776 0.00143967 0.106852 0.000843484 8.59434e-17 8.328e-15
    32 0.000750714 0.102493 0.00644056 0.000756826 0.0713168 0.000456685 2.75943e-16 5.37049e-15
    40 0.000430867 0.0675245 0.00509083 0.000460115 0.0510193 0.000286991 8.76007e-16 2.73806e-15
    48 0.000276193 0.0481754 0.00421193 0.000306103 0.038122 0.000195792 1.25266e-16 2.4751e-15

     | Show Table
    DownLoad: CSV
    Table 6.  Convergence order for Re=1000.0 and Rm=1000.0 for different h at T=1.0.
    1/h uL2 uH1 PL2 BL2 BH1
    16 3.07735 2.19229 1.26109 2.8519 2.42487
    24 2.6783 1.80614 1.04402 2.54119 1.77829
    32 2.60675 1.90339 1.16997 2.23523 1.40541
    40 2.48819 1.8701 1.05392 2.23021 1.50095
    48 2.43909 1.85191 1.03949 2.23536 1.59834

     | Show Table
    DownLoad: CSV

    Tables 1 and 2 give the numerical results for Re=1.0 and Rm=1.0. In order to show the effect of our method for high Reynolds number, we give the numerical results for Re=100.0, 1000.0 and Rm=100.0, 1000.0. Tables 3 and 4 give the numerical results for Re=100.0 and Rm=100.0. Tables 5 and 6 give the numerical results for Re=1000.0 and Rm=1000.0. It shows that our method is effect for high Reynolds numbers. It shows that the errors are goes small as the space step goes small and the convergence orders are optimal. We can see that |ΩBdx| and |Ωudx| are small, which means that our method can conserve the Gauss's law very well.

    In this paper, we have given a modified characteristics projection finite element method for the unsteady incompressible magnetohydrodynamics(MHD) equations. In this method, the modified characteristics finite element method and the projection method were combined for solving the incompressible MHD equations. Both the stability and the optimal error estimates in L2 and H1 norms for the modified characteristics projection finite element method have been given. In order to demonstrate the effectiveness of our method, some numerical results were given at the end of the manuscript.

    The authors would like to thank the anonymous referees for their valuable suggestions and comments, which helped to improve the quality of the paper. This work is supported in part by National Natural Science Foundation of China (No. 11971152), China Postdoctoral Science Foundation (No. 2018M630907) and the Key scientific research projects Henan Colleges and Universities (No. 19B110007).

    The authors declare there is no conflict of interest in this paper.



    [1] R. Moreau, Magnetohydrodynamics, Kluwer Academic Publishers, 1990.
    [2] W. Hughes, F. Young, The Electromagnetodynamics of Fluids, Wiley: New York, 1966.
    [3] M. D. Gunzburger, A. J. Meir, J. S. Peterson, On the existence, uniqueness, and finite element approximation of solutions of the equations of stationary, incompressible magnetohydrodynamics, Math. Comput., 56 (199), 523-563.
    [4] J. Gerbeau, C. Bris, T. Lelièvre, Mathematical methods for the Manetohydrodynamics of liquid metals, Oxford University Press, 2006.
    [5] M. Sermange, R. Temam, Some mathematical questions related to the magnetohydrodynamic equations, Comput. Compacts, 1 (1983), 212. doi: 10.1016/0167-7136(83)90286-X
    [6] U. Hasler, A. Schneebeli, D. Schötzau, Mixed finite element approximation of incompressible MHD problems based on weighted regularization, Appl. Numer. Math., 51 (2004), 19-45. doi: 10.1016/j.apnum.2004.02.005
    [7] J. L. Guermond, P. D. Minev, Mixed Finite Element approximation of an MHD problem involving conducting and insulating regions: The 2D case, ESAIM: Math. Modell. Numer. Analy., 36 (2002), 517-536. doi: 10.1051/m2an:2002024
    [8] J. L. Guermond, P. D. Minev, Mixed finite element approximation of an MHD problem involving conducting and insulating regions: The 3D case, Numer. Methods Partial Differential Equations, 19 (2002), 709-731. doi: 10.1002/num.10067
    [9] D. Schötzau, Mixed finite element methods for stationary incompressible magneto-hydrodynamics, Numer. Math., 96 (2004), 771-800. doi: 10.1007/s00211-003-0487-4
    [10] J. F. Gerbeau, A stabilized finite element method for the incompressible magnetohydrodynamic equations, Numerische Mathematik, 87 (2000), 83-111. doi: 10.1007/s002110000193
    [11] N. B. Salah, A. Soulaimani, W. G. Habashi, A finite element method for magnetohydrodynamics, Comput. Methods Appl. Mech. Eng., 190 (2001), 5867-5892. doi: 10.1016/S0045-7825(01)00196-7
    [12] A. I. Nesliturk, M. Tezer-Sezgin, Two-evel finite element method with a stabilizing subgrid for the incompressible MHD equations, Int. J. Numer. Methods Fluids, 62 (2010), 188-210.
    [13] X. J. Dong, Y. N. He, Y. Zhang, Convergence analysis of three finite element iterative methods for the 2D/3D stationary incompressible magnetohydrodynamics, Comput. Method. Appl. Math., 276 (2014), 287-311.
    [14] Y. N. He, Unconditional convergence of the Euler semi-implicit scheme for the three-dimensional incompressible MHD equations, IMA. J. Numer. Anal., 35 (2015), 767-801. doi: 10.1093/imanum/dru015
    [15] X. J. Dong, Y. N. He, Two-level newton iterative method for the 2D/3D stationary incompressible magnetohydrodynamics, J. Sci. Comput., 63 (2015), 426-451. doi: 10.1007/s10915-014-9900-7
    [16] R. Ingram, Numerical analysis of a finite element Crank-Nicolson discretization for MHD flows at small magnetic Reynolds numbers, Int. J. Numer. Analy. Model., 10 (2012), 74-98.
    [17] Y. F. Lei, Y. Yang, Z. Y. Si, Error estimate of fully discrete dc-fem for unsteady incompressible magnetohydrodynamics equations, Appl. Anal., 97 (2018), 2355-2376. doi: 10.1080/00036811.2017.1366990
    [18] Z. Y. Si, S. J. Jing, Y. X. Wang, Defect correction finite element method for the stationary incompressible magnetohydrodynamics equation, Appl. Math. Comput., 285 (2016), 184-194.
    [19] Z. Y. Si, C. Liu, Y. X. Wang, A semi-discrete defect correction finite element method for unsteady incompressible magnetohydrodynamics equations, Math. Method Appl. Sci., 40 (2017), 4179-4196. doi: 10.1002/mma.4296
    [20] J. E. Deng, Z. Y. Si, A decoupling penalty finite element method for the stationary incompressible magnetohydrodynamics equation, Int. J. Heat. Mass. Tran., 128 (2019), 601-612. doi: 10.1016/j.ijheatmasstransfer.2018.08.096
    [21] W. Layton, H. Tran, C. Trenchea, Numerical analysis of two parititioned methods for uncoupling evolutionary MHD flows, Numer. Methods Partial Differ. Equ., 30 (2014), 1083-1102. doi: 10.1002/num.21857
    [22] G. D. Zhang, Y. N. He, D. Yang, Analysis of coupling iterations based on the finite element method for stationary magnetohydrodynamics on a general domain, Comput. Math. Appl., 68 (2014), 770-788. doi: 10.1016/j.camwa.2014.07.025
    [23] G. D. Zhang, Y. N. He, Decoupled schemes for unsteady MHD equations. II: Finite element spatial discretization and numerical implementation, Comput. Math. Appl., 69 (2015), 1390-1406. doi: 10.1016/j.camwa.2015.03.019
    [24] R. A. Adams, Sobolev Space, In: Pure and Applied Mathematics, vol. 65, Academic Press, New York, 1975.
    [25] V. Girault, P. Raviart, Finite element methods for Navier-Stokes equations, 1986.
    [26] J. Heywood, G. John, The Navier-Stokes equations: On the existence, regularity and decay of solutions, Indiana University Math. J., 29 (1980), 639-681. doi: 10.1512/iumj.1980.29.29048
    [27] Y. N. He, Two-level method based on finite element and Crank-Nicolson extrapolation for the timedependent Navier-Stokes equations, SIAM J. Numer. Anal., 41 (2004), 1263-1285.
    [28] G. Galdi, An introduction to the mathematical theory of the Navier-Stokes equations: Steady-state problems, Springer, New Yourk, 2011.
    [29] Y. Achdou, J. L. Guermond, Convergence analysis of a finite element projection/Lagrange-Galerkin method for the incompressible Navier-Stokes equations, SIAM J. Numer. Anal., 37 (2000), 799-826. doi: 10.1137/S0036142996313580
    [30] R. Bermejo, P. Galán del Sastre, L. Saavedra, A second order in time modified Lagrange-Galerkin finite element method for the incompressible Navier-Stokes equations, SIAM J. Numer. Anal., 50 (2012), 3084-3109. doi: 10.1137/11085548X
    [31] Z. Y. Si, J. L. Wang, W. W. Sun, Unconditional stability and error estimates of modified characteristics FEMs for the Navier-Stokes equations, Numer. Math., 134 (2016), 139-161. doi: 10.1007/s00211-015-0767-9
    [32] F. Hecht, New development in FreeFem++, J. Numer. Math., 20 (2012), 251-265.
  • This article has been cited by:

    1. Ke Zhang, Haiyan Su, Xinlong Feng, Second order unconditional linear energy stable, rotational velocity correction method for unsteady incompressible magneto-hydrodynamic equations, 2022, 236, 00457930, 105300, 10.1016/j.compfluid.2021.105300
    2. Zhiyong Si, Mingyi Wang, Yunxia Wang, A projection method for the non-stationary incompressible MHD coupled with the heat equations, 2022, 428, 00963003, 127217, 10.1016/j.amc.2022.127217
    3. Huifang Zhang, Tong Zhang, Optimal error estimates of two‐level iterative finite element methods for the thermally coupled incompressible MHD with different viscosities, 2022, 0170-4214, 10.1002/mma.8800
    4. Zhiyong Si, Xiaonan Guo, Yunxia Wang, A modified characteristics projection finite element method for the non‐stationary incompressible thermally coupled MHD equations, 2023, 46, 0170-4214, 17422, 10.1002/mma.9508
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(3990) PDF downloads(293) Cited by(4)

Figures and Tables

Tables(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog