A defect-correction mixed finite element method for solving the time-dependent Johnson-Segalman viscoelastic equations in two dimensions is given. In the defect step, the constitutive equation is computed with the artificially reduced Weissenberg parameter for stability, and the resulting residual is corrected in the correction step on the same grid. A streamline upwind Petrov-Galerkin (SUPG) approximation is used to stabilize the hyperbolic character of the constitutive equation for the stress. We establish a priori error estimates for the defect step and the first correction step of the defect correction method. The derived theoretical results are supported by numerical tests.

1. Introduction

Time-dependent calculations of viscoelastic flows are important to the understanding of many problems in non-Newtonian fluid mechanics, particularity those related to flow instabilities. Error analysis of finite element approximations to time-dependent viscoelastic flow was first analyzed by Baranger and Wardi in [1], using an implicit Euler temporal discretization and a discontinuous Galerkin (DG) approximation for the hyperbolic constitutive equation. In [2, 3], Ervin and Miles have analyzed the problem using an implicit Euler time discretization and an SUPG discretization for the constitutive equation. The key of SUPG method, first introduced by Hughes and Brooks [4], is to stabilize the numerical solution by adding dissipative effects limited to the direction in which their influence is needed. Error analysis of a modified Euler-SUPG approximation to time-dependent viscoelastic flow problem was presented by Bensaada and Esselaoui in [5]. In [6], Ervin and Heuer have proposed a Crank-Nicolson time discretization scheme for this problem. A fractional step 𝜃-method for the time-dependent problem is described and analyzed in [7]. The reader may find more relevant work, for example, [8, 9] and the reference therein.

In addition, when the Weissenberg number increases, boundary layers for the stress develop. This will add to the difficulty of computing accurate numerical approximations. There are also many papers in developing stable numerical algorithms for high Weissenberg number flows (see [1015] and the references therein) over the years. In [1012], a defect-correction method has successfully applied in steady viscoelastic and Oseen-viscoelastic fluid flow model, aiming at high Weissenberg number. The basic idea of defect correction method for viscoelastic fluid flow is as follows. In the defect step, the Weissenberg number is reduced artificially by using a mesh-dependent parameter to obtain better convergence of an iteration scheme. Then, in a correction step, the initial approximation is improved using the residual correction. The defect correction method is an iterative improvement technique for increasing the accuracy of a numerical solution without applying a grid refinement. Due to its good efficiency, there are many works devoted to this method, for example, [1622].

Inspired by [2, 3, 1012, 21, 22], in this paper, we will use the defect correction method for solving the time-dependent Johnson-Segalman viscoelastic equations in two dimensions. Approximation in space is made by finite element (FE) method. Approximate stress, velocity, and pressure are, respectively, 𝑃1 continuous, 𝑃2 continuous, and 𝑃1 continuous. As the constitutive equation has hyperbolic character for the stress, we use streamline upwind Petrov-Galerkin (SUPG) formulation to stabilize the constitutive equation. Approximate in time is made by Euler scheme. Finally, we establish a priori error estimates for the defect step and the first correction step of the defect correction method. While in [22], the authors present this problem by DG stabilized style, we know that these two styles have many differences. Choice of the Johnson and Segalman model is not essential, and the results obtained can be extended to more realistic differential models likes those of Phan-Thien and Tanner, Giesekus and Oldroyd, provided a Newtonian viscosity is present.

Detailed descriptions, analysis, and numerical examples of the defect correction method for the viscoelastic fluid flow are presented in the rest of the paper as follows. In the next section, we present the time-dependent viscoelastic fluid flow of the Johnson-Segalman model and its variational formulation. In Section 3, we introduce finite element spaces and present the fully discrete approximation scheme (Euler in time, FE in space). Defect correction algorithm is given in Section 4. Then a priori error estimates for the defect step and the first correction step of the defect correction method are derived in Section 5. Numerical results are presented in Section 6. Finally, a summary and discussion of continuing work are presented in Section 7.

2. The Johnson-Segalman Model and Its Variational Formulation

In this section, we consider a fluid flowing in a connected, bounded polygonal domain Ω𝑑,𝑑=2 with the Lipschitz continuous boundary Γ.

The following notation will be used. 𝑅2 is equipped with cartesian coordinates 𝑥𝑖,𝑖=1,2. For a function 𝑢, 𝜕𝑢/𝜕𝑥𝑖 is written 𝑢,𝑖 and 𝜕𝑢/𝜕𝑡 is written 𝑢𝑡. Einstein's convention of summation is used. For a scalar function 𝑝, gradient of 𝑝 is a vector 𝑝 with (𝑝)𝑖=𝑝,𝑖. For a vector 𝐮, gradient of 𝐮 is a tensor 𝐮 with (𝐮)𝑖𝑗=𝐮𝑖,𝑗. For a vector function 𝐮, divergence of 𝐮 is a scalar 𝐮=𝐮𝑖,𝑖. with (𝐮)𝑖𝑗=𝐮𝑖,𝑗. For a tensor function 𝜏 and a vector function 𝐮, divergence of 𝜏 is a vector 𝜏 with (𝜏)𝑖=𝜏𝑖𝑗,𝑗 and 𝐮 denotes the operator 𝑢𝑖(𝜕/𝜕𝑥𝑖).

Let 𝑇>0 be a given final time, we consider the time-dependent Johnson-Segalman viscoelastic equations as follows:𝐮Re𝑡+𝑢𝑢𝜎2(1𝛼)𝐷(𝐮)+𝑝=𝐟,inΩ×(0,𝑇),𝜆𝜎𝑡+𝜎+𝜆(𝐮)𝜎+𝜆𝑔𝑎(𝜎,𝐮)2𝛼𝐷(𝐮)=0,inΩ×(0,𝑇),div𝐮=0,inΩ×(0,𝑇),𝐮=𝟎,onΓ×(0,𝑇),𝐮(0,𝑥)=𝐮0(𝑥),𝜎(0,𝑥)=𝜎0(𝑥),onΩ×{0},(2.1) where 𝜆 is the Weissenberg number, Re is the Reynolds number, 0<𝛼<1 may be considered as the fraction of viscoelastic viscosity (𝛼=1 for Maxwell's model and this case is excluded here, see reference [23] for a description of various models), and 𝐟 the body forces. The unknowns are the fluid velocity vector 𝐮=(𝑢1,𝑢2), the pressure 𝑝, and the stress 𝜎, which is the viscoelastic part of the total stress tensor 𝜎tot=𝜎+2(1𝛼)𝐷(𝐮)𝑝𝐈. In (2.1), 𝐷(𝐮)=(1/2)(𝐮+𝐮𝑇) is the rate of the strain tensor, for all 𝑎[1,1], 𝑔𝑎(𝜎,𝐮) is defined by𝑔𝑎(𝜎,𝐮)=1𝑎2𝜎𝐮+(𝐮)𝑇𝜎1+𝑎2(𝐮)𝜎+𝜎(𝐮)𝑇.(2.2) Note that an Oldroyd B constitutive model is obtained when 𝑎=1 in 𝑔𝑎(𝜎,𝐮).

The 𝐿2(Ω) inner product is denoted by (,), and the 𝐿𝑝(Ω) norm by 𝐿𝑝, with the special cases of 𝐿2(Ω) and 𝐿(Ω) norms being written as and . For 𝑘𝑁, we denote the norm associated with the Sobolev space 𝑊𝑚,𝑝(Ω) by 𝑊𝑚,𝑝, with the special case 𝑊𝑚,2(Ω) being written as 𝐻𝑚(Ω) with the norm 𝑚 and seminorm ||𝑚. In order to introduce a variational formulation, we set the spaces 𝑋, 𝑄, 𝑆, 𝑉 as follows:𝑋=𝐻10(Ω)2=𝑣𝐻1(Ω)2,𝐯=𝟎onΓ𝑄=𝐿20(Ω)=𝑞𝐿2(Ω);Ω,𝜏𝑞𝑑𝑥=0𝑆=𝜏=𝑖𝑗;𝜏𝑖𝑗=𝜏𝑗𝑖;𝜏𝑖𝑗𝐿2𝜏(Ω);𝑖,𝑗=1,2𝜏=𝑖𝑗;𝐯𝜏𝐿2(,Ω),𝐯𝑋𝑉=𝐯𝑋;Ω.𝑞(𝐯)𝑑𝑥=0,𝑞𝑄(2.3)

The corresponding weak formulation of problem (2.1) is then obtained: find (𝜎,𝐮,𝑝)(𝑆,𝑋,𝑄) such that for all (𝜏,𝐯,𝑞)(𝑆,𝑋,𝑄)Re𝐮𝑡+𝜆𝜎+Re𝐮𝐮,𝐯(𝜎,𝐷(𝐯))+2(1𝛼)(𝐷(𝐮),𝐷(𝐯))+(𝑝,𝐯)=(𝐟,𝐯),𝑡𝑔,𝜏+(𝜎,𝜏)+𝜆((𝐮)𝜎,𝜏)+𝜆𝑎(𝜎,𝐮),𝜏2𝛼(𝐷(𝐮),𝜏)=𝟎,(𝑞,𝐮)=0.(2.4)

Using the weak divergence free space 𝑉, the weak formulation (2.4) can be written as follows: find (𝜎,𝐮)(𝑆,𝑉) such that for all (𝜏,𝐯)(𝑆,𝑉)Re𝐮𝑡+𝜆𝜎+Re𝐮𝐮,𝐯(𝜎,𝐷(𝐯))+2(1𝛼)(𝐷(𝐮),𝐷(𝐯))=(𝐟,𝐯),𝑡𝑔,𝜏+(𝜎,𝜏)+𝜆((𝐮)𝜎,𝜏)+𝜆𝑎(𝜎,𝐮),𝜏2𝛼(𝐷(𝐮),𝜏)=𝟎.(2.5)

Existence results for problem (2.1) are proved in [24] for the "slow flow" model of (2.1) (i.e., the 𝐮𝐮 term in momentum equation is ignored). They are of two kinds: local existence in time of strong solutions in space variable in a 𝐶3 domain, and global existence (in time) of a unique solution for 𝐮 and 𝜎 under a small data assumption on 𝐟, 𝐟𝑡, 𝐮0, 𝜎0. For a more complete discussion of existence and uniqueness issues, see [25].

3. Finite Element Approximation

In this section, we present a fully discrete approximation to (2.5). We begin by describing the finite element approximation framework.

Suppose 𝑇 is a uniformly regular triangulation of Ω such that Ω={𝐾𝐾𝑇} and assume that there exist positive constants 𝐶1,𝐶2 such that 𝐶1𝐾𝐶2𝜌𝐾, where 𝐾 is the diameter of 𝐾, 𝜌𝐾 is the diameter of the greatest ball included in 𝐾, and =max𝐾𝑇𝐾. Throughout the paper, the constants 𝐶,𝐶1,𝐶2, denote different constants which are independent of mesh size and time step 𝑘. We use the classical Taylor-Hood FE for the approximation in space of (𝐮,𝑝)𝑃2-continuous in velocity, 𝑃1-continuous in pressure; we consider 𝑃1-continuous approximation of the stress 𝜎. The corresponding FE spaces are 𝑋=𝐯𝑋𝐶0(Ω)2;𝐯|𝐾𝑃2(𝐾)2,𝐾𝑇,𝑆=𝜏𝑆𝐶0(Ω)2×2;𝜏|𝐾𝑃1(𝐾)2×2;𝐾𝑇,𝑄=𝑞𝑄𝐶0(Ω);𝑞|𝐾𝑃1(𝐾);𝐾𝑇,𝑉=𝐯𝑋;(𝑞,𝑣)=0,𝑞𝑄,(3.1) where 𝑃𝑤(𝐾) denotes the space of polynomials of degree 𝑤 on 𝐾𝑇. It is well known [26, 27] that the Taylor-Hood pair (𝑋,𝑄) satisfies the discrete inf-sup (or LBB) condition.

To obtain the fully discrete approximation, the time derivatives are replaced by backward Euler differences, and the nonlinear terms are lagged. In order to stabilize the hyperbolic constitutive equation, an SUPG formulation is used to avoid spurious oscillations in the numerical approximation. Let 𝑁 be an integer, we divide the interval [0,𝑇] into 𝑁 intervals of equal length 𝑘 and denote 𝑡𝑛=𝑛𝑘. For notational convenience, we denote 𝑣𝑛=𝑣(,𝑡𝑛) and 𝑑𝑡𝑔𝑠𝑔𝑡=𝑠𝑡𝑔𝑠1𝑘,𝑏(𝐮,𝜎,𝜏)=(𝐮𝜎,𝜏),𝑐(𝐮,𝐯,𝐰)=(𝐮𝐯,𝐰).(3.2) The following norms are used in the analysis: ||||𝑔,𝑚=max1𝑛𝑁𝑔𝑛𝑚,||||𝑔0,𝑚=𝑁𝑛=1𝑘𝑔𝑛2𝑚1/2,(3.3) and when 𝜓(𝐱,𝑡) is defined on the entire time interval (0,𝑇), we denote 𝜓,𝑚=sup0𝑡𝑇𝜑(,𝑡)𝑚,𝜓20,𝑚=𝑇0𝜓(,𝑡)2𝑚𝑑𝑡,𝜓(𝑡)=𝜓(,𝑡).(3.4)

Then, based on SUPG formulation, the fully discrete approximating systems of (2.5) is the following;

given 𝐮,0=𝐮0,,𝐮,𝑛; 𝜎,0=𝜎0,,𝜎,𝑛, 𝑛=0,1,2,,𝑁1, find (𝑢,𝑛+1,𝜎,𝑛+1)(𝑉,𝑆) such that for all (𝐯,𝜏)(𝑉,𝑆)𝑑Re𝑡𝐮,𝑛+1𝐮,𝐯+Re𝑐,𝑛,𝐮,𝑛+1𝐷𝐮,𝐯+2(1𝛼),𝑛+1+𝜎,𝐷(𝐯),𝑛+1=𝐟,𝐷(𝐯)𝑛+1,𝜆𝑑,𝐯𝑡𝜎,𝑛+1+𝜎,𝜏,𝑛+1,𝜏𝐮+𝜆𝑏,𝑛,𝜎,𝑛+1,𝜏𝐷𝐮2𝛼,𝑛+1,𝜏𝑔=𝜆𝑎𝜎,𝑛,𝐮,𝑛,𝜏,(3.5) where 𝜏=𝜏+𝜈𝜏𝑛,𝑢,𝜏𝑛,𝑢=𝑢,𝑛𝜏, 𝜈 is a small positive constant.

The goal of the parameter 𝜈 is used to suppress the production of spurious oscillations in the approximation. The discretization of the constitutive equation is the usual Galerkin finite element method if 𝜈=0.

The existence of a unique solution and a priori error estimate to (3.5) can be found in [2, 3, 5].

4. Defect Correction Method

In this section, our defect correction method used in computing the solution to (3.5) is described as follows.

Algorithm 1 (Defect-correction method). Step 1. Solve the defected problem: find (𝐮1,𝑛+1,𝜎1,𝑛+1)(𝑉,𝑆) such that for all (𝐯,𝜏)(𝑉,𝑆)𝑑Re𝑡𝐮1,𝑛+1𝐮,𝐯+Re𝑐1,𝑛,𝐮1,𝑛+1𝐷𝐮,𝐯+2(1𝛼)1,𝑛+1+𝜎,𝐷(𝐯)1,𝑛+1=𝐟,𝐷(𝐯)𝑛+1,𝜆𝑑,𝐯(4.1)𝑡𝜎1,𝑛+1+𝜎,𝜏1,𝑛+1,𝜏+𝜆1𝑏𝐮1,𝑛,𝜎1,𝑛+1,𝜏𝐷𝐮2𝛼1,𝑛+1,𝜏𝑔=𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜏,(4.2) where 𝐸 is chosen such that 𝜆1=𝜆𝐸>0.Step 2. For 𝑖=1,2,, solve the correction problem: find (𝐮,𝑛+1𝑖+1,𝜎,𝑛+1𝑖+1)(𝑉,𝑆) such that for all (𝐯,𝜏)(𝑉,𝑆)𝑑Re𝑡𝐮,𝑛+1𝑖+1𝐮,𝐯+Re𝑐,𝑛𝑖+1,𝐮,𝑛+1𝑖+1𝐷𝐮,𝐯+2(1𝛼),𝑛+1𝑖+1+𝜎,𝐷(𝐯),𝑛+1𝑖+1=𝐟,𝐷(𝐯)𝑛+1,𝜆𝑑,𝐯(4.3)𝑡𝜎,𝑛+1𝑖+1+𝜎,𝜏,𝑛+1𝑖+1,𝜏+𝜆1𝑏𝐮,𝑛𝑖+1,𝜎,𝑛+1𝑖+1,𝜏𝐷𝐮2𝛼,𝑛+1𝑖+1,𝜏𝑔=𝜆𝑎𝜎,𝑛𝑖+1,𝐮,𝑛𝑖+1,𝜏𝜆𝜆1𝑏𝐮,𝑛𝑖+1,𝜎𝑖,𝑛+1,𝜏.(4.4) The initial value approximations are taken to be 𝐮1,0=𝐮,0𝑖+1=𝐮0, and 𝜎1,0=𝜎,0𝑖+1=𝜎0. In order to ensure computability of the algorithm, we begin by showing that (4.1)-(4.2) and (4.3)-(4.4) are uniquely solvable for (𝐮𝑗,𝑛+1,𝜎𝑗,𝑛+1),𝑗=1,2 at each time step 𝑛+1. We use the following induction hypothesis, which simply states that the computed iterates (𝐮𝑗,𝑛+1,𝜎𝑗,𝑛+1), 𝑗=1,2 are bounded independent of and 𝑛: 𝑢(𝐈𝐇𝟏)𝑗,𝑠,𝜎𝑗,𝑠𝐾,𝑗=1,2,0𝑠𝑛.(4.5)

In next section, we will prove the induction hypothesis (𝐈𝐇𝟏) is right.

Lemma 4.1. Suppose (𝐈𝐇𝟏) is true. For sufficiently small step size 𝑘, then Step 1 of Algorithm 1 admits a unique solution (𝑢1,𝑛+1,𝜎1,𝑛+1)(𝑉,𝑆).

Proof. Letting 𝐯=𝐮1,𝑛+1,𝜏=𝜎1,𝑛+1, multiplying (4.1) by 2𝛼, and adding to (4.2), we get 𝐴𝐮1,𝑛+1,𝜎1,𝑛+1;𝐮1,𝑛+1,𝜎1,𝑛+1𝑓=2𝛼𝑛+1,𝐮1,𝑛+1+2𝛼Re𝑘𝐮1,𝑛,𝐮1,𝑛+1𝑔𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜎1,𝑛+1+𝜆𝑘𝜎1,𝑛,𝜎1,𝑛+1,(4.6) where the bilinear form 𝐴(𝐮,𝜎;𝐯,𝜏) is defined by 𝐴(𝐮,𝜎;𝐯,𝜏)=2𝛼Re𝑘(𝐮,𝐯)+2𝛼Re𝑐(𝐮𝑛,𝐮,𝐯)+4𝛼(1𝛼)(𝐷(𝐮),𝐷(𝐯))+2𝛼(𝜎,𝐷(𝐯))+𝜎,𝜏+𝜆𝑘(𝜎,𝜏)+𝜆1𝑏(𝐮𝑛,𝜎,𝜏)+𝜆1𝑏(𝐮𝑛,𝜎,𝜈𝑢𝑛𝜏)2𝛼(𝐷(𝐮),𝜎)2𝛼(𝐷(𝐮),𝜈𝑢𝑛𝜏).(4.7)
As (4.1)-(4.2) represent a finite system of liner equations, the positivity of 𝐴(𝐮1,𝑛+1,𝜎1,𝑛+1;𝐮1,𝑛+1,𝜎1,𝑛+1) is a sufficient condition for the existence and uniqueness of (𝐮1,𝑛+1,𝜎1,𝑛+1).
We now estimate each term of 𝐴(𝐮1,𝑛+1,𝜎1,𝑛+1;𝐮1,𝑛+1,𝜎1,𝑛+1). We have |||𝐮2𝛼Re𝑐1,𝑛,𝐮1,𝑛+1,𝐮1,𝑛+1|||𝜖1𝐷𝐮1,𝑛+12+𝛼2𝑑𝐾2𝐶2𝐾Re2𝜖1𝐮1,𝑛+12,|||𝜎1,𝑛+1,𝜈𝜎1𝑛,𝑢|||𝜎1,𝑛+120+𝜈24𝜎1𝑛,𝑢2,|||𝜆1𝑏𝐮1,𝑛,𝜎1,𝑛+1,𝜎1,𝑛+1|||𝜖2𝜎1𝑛,𝑢2+𝜆214𝜖2𝜎1,𝑛+12,𝜆1𝑏𝐮1,𝑛,𝜎1,𝑛+1,𝜈𝜎1𝑛,𝑢=𝜆1𝜈𝜎1𝑛,𝑢2,|||𝐷𝐮2𝛼1,𝑛+1,𝜈𝜎1𝑛,𝑢|||𝜖3𝐷𝐮1,𝑛+12+𝛼2𝜈2𝜖3𝜎1𝑛,𝑢2,(4.8) where we have used the induction hypothesis (𝐈𝐇𝟏)and the Korn's lemma, 𝐶𝐾 is Korn constant.
Substituting all above equations into 𝐴(,;,), we have 𝐴𝐮1,𝑛+1,𝜎1,𝑛+1,𝐮1,𝑛+1,𝜎1,𝑛+12𝛼Re𝑘𝛼2𝑑𝐾2𝐶2𝐾Re2𝜖1𝐮1,𝑛+12+𝜆𝑘𝜆214𝜖2𝜎1,𝑛+12+4𝛼(1𝛼)𝜖1𝜖3𝐮𝐷1,𝑛+12+𝜆1𝜈𝜖2𝜈24𝛼2𝜈2𝜖3𝜎1𝑛,𝑢2.(4.9) Choosing 𝜈(2𝜆1(1𝛼))/(1+3𝛼), 𝑘min{2𝜆𝜈/𝜆1,2(1𝛼)/𝑑𝐾2𝐶2𝐾Re} and 𝜖1=𝜖3=𝛼(1𝛼), 𝜖2=𝜆1𝜈/2, the bilinear form 𝐴(,;,) is positive definite. Hence, for 𝑘 chosen sufficiently small, Step 1 of Algorithm 1 admits a unique solution (𝑢1,𝑛+1,𝜎1,𝑛+1)(𝑉,𝑆).

Similarity as Lemma 4.1, we can derive that Step 2 of Algorithm 1 admits a unique solution (𝑢2,𝑛+1,𝜎2,𝑛+1)(𝑉,𝑆).

5. A Priori Error Estimates

In this section, we explore the error estimates in approximating the true solution (𝐮,𝑝,𝜎) of (2.1) by the defect step approximation (𝑢1,𝑝1,𝜎1) and one correction step approximation (𝑢2,𝑝2,𝜎2). The main results of this section are presented in the following theorem.

Theorem 5.1. Suppose that (2.1) has a solution (𝐮,𝜎,𝑝)𝐶2(0,𝑇;𝐻3)×𝐶2(0,𝑇;𝐻2)×𝐶1(0,𝑇;𝐻2). In addition, assume that 𝑘, 𝜈𝑐𝑑/2, and 𝐮, 𝐮, 𝜎, 𝜎𝑀 for all 𝑡[0,𝑇], (𝐮1,𝜎1) and (𝐮2,𝜎2) satisfy (4.1)-(4.2) and (4.3)-(4.4), respectively. Then, there exists a constant 𝐶=𝐶(Re,𝛼,Ω,𝑇,𝐮,𝑝,𝜎,𝐟,𝜆1,𝜆), such that ||𝐮𝐮1||,0+||𝜎𝜎1||,0+||𝐮𝐮1||0,1+||𝜎𝜎1||0,0𝐵(𝑘,,𝜈),(5.1)||𝐮𝐮2||,0+||𝜎𝜎2||,0+||𝐮𝐮2||0,1+||𝜎𝜎2||0,0𝐵(𝑘,,𝜈),(5.2) where 𝐵(𝑘,,𝜈)=𝐶2||||𝐮0,3+3||𝐮𝑡||0,3||||+𝜎0,2+2||𝜎𝑡||0,2+2||||𝑝0,2+𝐶3||||𝐮,3+2||||𝜎,2𝐮+𝐶𝑘𝑡0,1+𝐮𝑡𝑡0,0+𝜎𝑡0,1+𝜎𝑡𝑡0,0||𝜎+𝐶𝜈𝑡||0,1+||𝜎𝑡||,0+𝐶𝑇𝜆𝜆1.(5.3)

Remark 5.2. As mentioned in Algorithm 1, 𝜆𝜆1=𝐸, thus, for 𝑖=1,2, the essence of the estimates (5.1) and (5.2) is ||𝐮𝐮𝑖||,0+||𝜎𝜎𝑖||,0+||𝐮𝐮𝑖||0,1+||𝜎𝜎𝑖||0,0𝐶(+𝑘+𝜈).(5.4)

Before deriving the estimates (5.2), (5.9), we first introduce some notation and some approximating properties of finite element spaces. Let 𝐮𝑛=𝐮(𝑡𝑛),𝜎𝑛=𝜎(𝑡𝑛) represent the solution of (2.5) at time 𝑡𝑛=𝑛𝑘. Assume (𝐮,𝜎,𝑝)𝐶2(0,𝑇;𝐻3)×𝐶2(0,𝑇;𝐻2)×𝐶1(0,𝑇;𝐻2). Let (̃𝐮𝑛,̃𝑝𝑛) denote the Stokes projection of (𝐮𝑛,𝑝𝑛) into (𝑉,𝑄), and 𝜎𝑛 a Clement interpolant of 𝜎𝑛 [26, 28]. The following approximating properties are right.𝐮𝑛̃𝐮𝑛+(𝐮𝑛̃𝐮𝑛)𝐶3𝐮𝑛3,𝑝𝑛̃𝑝𝑛𝐶2𝑝𝑛2,𝜎𝑛𝜎𝑛𝜎+𝑛𝜎𝑛𝐶2𝜎𝑛2.(5.5)

Let (𝐮1,𝑛,𝜎1,𝑛) and (𝐮2,𝑛,𝜎2,𝑛) denote the solutions of (4.1)-(4.2) and (4.3)-(4.4), respectively. We denote Φ𝑛𝑗,𝑌𝑛𝑗,𝑒𝑛𝑗,Ψ𝑛𝑗,𝑈𝑛𝑗,𝜖𝑛𝑗, 𝑗=1,2 as follows:Φ𝑛𝑗=𝐮𝑛̃𝐮𝑛,𝑌𝑛𝑗=̃𝐮𝑛𝐮𝑗,𝑛,𝑒𝑛𝑗=Φ𝑛𝑗+𝑌𝑛𝑗=𝐮𝑛𝐮𝑗,𝑛,Ψ𝑛𝑗=𝜎𝑛𝜎𝑛,𝑈𝑛𝑗=𝜎𝑛𝜎𝑗,𝑛,𝜖𝑛𝑗=Ψ𝑛𝑗+𝑈𝑛𝑗=𝜎𝑛𝜎𝑗,𝑛.(5.6)

In order to establish Theorem 5.1, we need another induction hypothesis, that is, (𝐈𝐇𝟐)𝑠1𝑛=0𝑘𝑌𝑛𝑗1,𝑗=1,2.(5.7)

The induction hypothesis (𝐈𝐇𝟐) will be proved later.

Proof. We first give the profile of the proof. The proof of Theorem 5.1 is established in two steps.Step 1. Prove the error estimate (5.1) is right. We divide Step 1 into two substep.
Substep 1.1. Under the induction hypothesis (𝐈𝐇𝟏) and (𝐈𝐇𝟐) for 𝑗=1, we prove the error estimate (5.1) is true.Substep 1.2. We prove that the two induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐) for 𝑗=1 are right.
Step 2. Show that the error estimate (5.2) is true.
Now, we start to prove Theorem 5.1.
Step 1. Prove the error estimate (5.1) is right.
Substep 1.1. Under the induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐), we prove the error estimate (5.1) is true.
As (𝜎,𝐮,𝑝) being the exact solution of (2.1), it satisfies the following consistency equation: for all (𝜏,𝐯)𝑆×𝑉, in particular, at time 𝑡=𝑡𝑛+1 (noting that we denote 𝑣𝑛+1=𝑣(,𝑡𝑛+1)), 𝐮2𝛼Re𝑡𝑛+1+,𝑣𝜆𝜎𝑡𝑛+1,𝜏+2𝛼𝑐Re𝐮𝑛+1,𝐮𝑛+1𝜎,𝐯+2𝛼𝑛+1𝐷𝐮,𝐷(𝐯)+4𝛼(1𝛼)𝑛+1𝐮,𝐷(𝐯)+𝜆𝑏𝑛+1,𝜎𝑛+1,𝜏+𝜎𝑛+1,𝜏𝐷𝐮2𝛼𝑛+1,𝜏𝑓=2𝛼𝑛+1𝑝,𝐯+2𝛼𝑛+1𝑔,𝐯𝜆𝑎𝜎𝑛+1,𝐮𝑛+1,𝜏.(5.8) Equation (5.8) can also be written as 𝑑2𝛼Re𝑡𝐮𝑛+1+,𝐯𝜆𝑑𝑡𝜎𝑛+1,𝜏+2𝛼𝑐Re𝐮1,𝑛,𝐮𝑛+1𝐷𝐮,𝐯+4𝛼(1𝛼)𝑛+1𝜎,𝐷(𝐯)+2𝛼𝑛+1,𝐷(𝐯)+𝜆1𝑏𝐮1,𝑛,𝜎𝑛+1,𝜏+𝜎𝑛+1,𝜏𝐷𝑢2𝛼𝑛+1,𝜏𝐟=2𝛼𝑛+1𝑝,𝐯+2𝛼𝑛+1𝑔,𝐯𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜏𝑑+2𝛼Re𝑡𝐮𝑛+1𝐮,𝐯2𝛼Re𝑡𝑛+1,𝐯+2𝛼𝑐Re𝐮1,𝑛,𝐮𝑛+1,𝐯2𝛼𝑐Re𝐮𝑛+1,𝐮𝑛+1𝑔,𝐯+𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜏𝑔𝜆𝑎𝜎𝑛+1,𝐮𝑛+1,𝜏+𝜆𝑑𝑡𝜎𝑛+1,𝜏𝜆𝜎𝑡𝑛+1,𝜏+𝜆1𝑏𝐮1,𝑛,𝜎𝑛+1,𝜏𝐮𝜆𝑏𝑛+1,𝜎𝑛+1,𝜏.(5.9) Multiplying (4.1) by 2𝛼 and adding to (4.2), for all (𝜏,𝐯)𝑆×𝑉, we get 𝑑2𝛼Re𝑡𝐮1,𝑛+1+,𝐯𝜆𝑑𝑡𝜎1,𝑛+1,𝜏+2𝛼𝑐Re𝐮1,𝑛,𝐮1,𝑛+1𝜎,𝐯+2𝛼1,𝑛+1𝐷𝐮,𝐷(𝐯)+4𝛼(1𝛼)1,𝑛+1,𝐷(𝐯)+𝜆1𝑏𝐮1,𝑛,𝜎1,𝑛+1,𝜏+𝜎1,𝑛+1,𝜏𝐷𝐮2𝛼1,𝑛+1,𝜏𝑓=2𝛼𝑛+1𝑔,𝐯𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜏.(5.10) Subtracting (5.9) from (5.10), for all (𝜏,𝐯)𝑆×𝑉, we obtain the following equation for 𝑒1𝑛+1 and 𝜖1𝑛+1: 𝑑2𝛼Re𝑡𝑒1𝑛+1+,𝐯𝜆𝑑𝑡𝜖1𝑛+1,𝜏+2𝛼𝑐Re𝐮1,𝑛,𝑒1𝑛+1𝐷𝑒,𝑣+4𝛼(1𝛼)1𝑛+1𝜖,𝐷(𝐯)+2𝛼1𝑛+1,𝐷(𝐯)+𝜆1𝑏𝐮1,𝑛,𝜖1𝑛+1,𝜏+𝜖1𝑛+1,𝜏𝐷𝑒2𝛼1𝑛+1,𝜏𝑑=2𝛼Re𝑡𝐮𝑛+1𝐮𝑡𝑛+1𝐮,𝐯+2𝛼Re𝑐1,𝑛𝐮𝑛+1,𝑢𝑛+1,𝐯𝜆𝜎𝑡𝑛+1,𝜏+𝜆𝑑𝑡𝜎𝑛+1𝐮,𝜏𝜆𝑏𝑛+1,𝜎𝑛+1,𝜏+𝜆1𝑏𝐮1,𝑛,𝜎𝑛+1,𝜏𝑔+𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝜏𝑔𝜆𝑎𝜎𝑛+1,𝐮𝑛+1,𝜏𝑝+2𝛼𝑛+1.,𝐯(5.11) Substituting 𝑒1𝑛+1=𝑌1𝑛+1+Φ1𝑛+1,𝜖1𝑛+1=𝑈1𝑛+1+Ψ1𝑛+1, 𝐯=𝑌1𝑛+1 and 𝜏=𝑈1𝑛+1 into (5.11), we get 𝑑2𝛼Re𝑡𝑌1𝑛+1,𝑌1𝑛+1𝑑+𝜆𝑡𝑈1𝑛+1,𝑈1𝑛+1𝑈+2𝛼1𝑛+1𝑌,𝐷1𝑛+1+2𝛼𝑐Re𝐮1,𝑛,𝑌1𝑛+1,𝑌1𝑛+1𝐷𝑌+4𝛼(1𝛼)1𝑛+1𝑌,𝐷1𝑛+1+𝑈1𝑛+1,𝑈1𝑛+1+𝜆1𝑏𝐮1,𝑛,𝑈1𝑛+1,𝑈1𝑛+1𝐷𝑌2𝛼1𝑛+1,𝑈1𝑛+1𝑌=𝐅1𝑛+1,𝑈1𝑛+1,(5.12) where 𝐅𝑌1𝑛+1,𝑈1𝑛+1𝑝=2𝛼𝑛+1,𝑌1𝑛+1𝑑+2𝛼Re𝑡𝐮𝑛+1𝐮𝑡𝑛+1,𝑌1𝑛+1𝜆𝜎𝑛+1t,𝑈1𝑛+1𝐮+2𝛼Re𝑐1,𝑛𝐮𝑛+1,𝐮𝑛+1,𝑌1𝑛+1𝑔+𝜆𝑎𝜎1,𝑛,𝐮1,𝑛,𝑈1𝑛+1+𝜆𝑑𝑡𝜎𝑛+1,𝑈1𝑛+1𝑔𝜆𝑎𝜎𝑛+1,𝐮𝑛+1,𝑈1𝑛+12𝛼𝑐Re𝐮1,𝑛,Φ1𝑛+1,𝑌1𝑛+1𝐷Φ4𝛼(1𝛼)1𝑛+1𝑌,𝐷1𝑛+1Ψ2𝛼1𝑛+1𝑌,𝐷1𝑛+1Ψ1𝑛+1,𝑈1𝑛+1𝐮𝜆𝑏𝑛+1,𝜎𝑛+1,𝑈1𝑛+1+𝜆1𝑏𝐮1,𝑛,𝜎𝑛+1,𝑈1𝑛+1𝜆1𝑏𝐮1,𝑛,Ψ1𝑛+1,𝑈1𝑛+1𝐷Φ+2𝛼1𝑛+1,𝑈1𝑛+1𝑑2𝛼Re𝑡Φ1𝑛+1,𝑌1𝑛+1𝑑𝜆𝑡Ψ1𝑛+1,𝑈1𝑛+1.(5.13) Using the identity (𝑎𝑏,𝑎)=(1/2)(𝑎2𝑏2+𝑎𝑏2), we obtain 𝛼Re𝑘𝑌1𝑛+12𝑌𝑛12+𝑌1𝑛+1𝑌𝑛12𝐷𝑌+4𝛼(1𝛼)1𝑛+12+𝜆𝑈2𝑘1𝑛+12𝑈𝑛12+𝑈1𝑛+1𝑈𝑛12+𝑈1𝑛+12+𝜈𝜆1𝑈1𝑛,𝑢2|||𝑐𝐮2𝛼Re1,𝑛,𝑌1𝑛+1,𝑌1𝑛+1|||||𝑈+𝜈1𝑛+1,𝑈1𝑛,𝑢||+𝜆1|||𝑏𝐮1,𝑛,𝑈1𝑛+1,𝑈1𝑛+1|||||𝐷𝑌+2𝛼𝜈1𝑛+1,𝑈1𝑛,𝑢||𝑌+𝐅1𝑛+1,𝑈1𝑛+1.(5.14) First multiply (5.14) by 𝑘 and sum it form 𝑛=0 to 𝑠1, we obtain 𝛼Re𝑠1𝑛=0𝑌1𝑛+12𝑌𝑛12+𝑌1𝑛+1𝑌𝑛12+4𝛼(1𝛼)𝑘𝑠1𝑛=0𝐷𝑌1𝑛+12+𝜆2𝑠1𝑛=0U1𝑛+12𝑈𝑛12+𝑈1𝑛+1𝑈𝑛12+𝑘𝑠1𝑛=0𝑈1𝑛+12+𝜈𝜆1𝑘𝑠1𝑛=0𝑈1𝑛,𝑢2𝑘𝑠1𝑛=0|||𝑐𝐮2𝛼Re1,𝑛,𝑌1𝑛+1,𝑌1𝑛+1|||||𝑈+𝜈1𝑛+1,𝑈1𝑛,𝑢||+𝜆1|||𝑏𝐮1,𝑛,𝑈1𝑛+1,𝑈1𝑛+1|||||𝐷𝑌+2𝛼𝜈1𝑛+1,𝑈1𝑛,𝑢||+𝑘𝑠1𝑛=0𝐅𝑌1𝑛+1,𝑈1𝑛+1.(5.15) Then, we get, remarking that 𝑌01=0 and 𝑈01=0, 𝑌𝛼Re𝑠12+𝑠1𝑛=0𝑌1𝑛+1𝑌𝑛12+𝜆2𝑈𝑠12+𝑠1𝑛=0𝑈1𝑛+1𝑈𝑛12+𝑘𝑠1𝑛=0𝑈1𝑛+12+4𝛼(1𝛼)𝑘𝑠1𝑛=0𝐷𝑌1𝑛+12+𝜈𝜆1𝑘𝑠1𝑛=0𝑈1𝑛,𝑢2𝑘𝑠1𝑛=0|||𝑐𝐮2𝛼Re1,𝑛,𝑌1𝑛+1,𝑌1𝑛+1|||||𝑈+𝜈1𝑛+1,𝑈1𝑛,𝑢||+𝜆1|||𝑏𝐮1,𝑛,𝑈1𝑛+1,𝑈1𝑛+1|||||𝐷𝑌+2𝛼𝜈1𝑛+1,𝑈1𝑛,𝑢||+𝑘𝑠1𝑛=0𝐅𝑌1𝑛+1,𝑈1𝑛+1.(5.16)
For controlling each term on the right-hand side(RHS) of (5.16), the assumption (𝐈𝐇𝟏) and (𝐈𝐇𝟐) are needed here. We will prove the two induction hypotheses in the next subsection. Let us estimate each term of the RHS of (5.16). For details, please see [2, 3, 6]. We start from the first four terms of the RHS of (5.16); |||𝑐𝐮2𝛼Re1,𝑛,𝑌1𝑛+1,𝑌1𝑛+1|||4𝛼2Re2𝐶2𝐾𝜖1𝐷𝑌1𝑛+12+𝑑𝐾24𝜖1𝑌1𝑛+12,𝜈||𝑈1𝑛+1,𝑈1𝑛,𝑢||𝑈1𝑛+12+𝜈24𝑈1𝑛,𝑢2,|||𝜆1𝑏𝐮1,𝑛,𝑈1𝑛+1,𝑈1𝑛+1|||𝜆12𝑑𝑌𝑛1𝑈1𝑛+12+𝑀𝑈1𝑛+12,||𝐷𝑌2𝛼𝜈1𝑛+1,𝑈1𝑛,𝑢||4𝛼2𝜖2𝐷𝑌1𝑛+12+𝜈24𝜖2𝑈1𝑛,𝑢2.(5.17) Now, we estimate each term of 𝐅(𝑌1𝑛+1,𝑈1𝑛+1). As 𝐅 has many terms and in order to make the proof more clear, we divide the terms of 𝐅 into four parts, that is, 𝑐(,,) terms of 𝐅, 𝑏(,,) terms of 𝐅, 𝑔𝑎 terms of 𝐅, other terms of 𝐅. Firstly, we estimate 𝑐(,,) terms of 𝐅; |||𝑐𝐮2𝛼Re1,𝑛,Φ1𝑛+1,𝑌1𝑛+1|||2𝛼Re2𝑌1𝑛+12+𝐾2𝑑4Φ1𝑛+12.|||𝑐𝐮2𝛼Re1,𝑛𝐮𝑛+1,𝐮𝑛+1,𝑌1𝑛+1|||2𝛼Re2𝑑2𝑀22𝑌𝑛12+32𝛼𝑌1𝑛+12+Re2𝑑2𝑀22Φ𝑛12+Re22𝑑2𝑀2𝑘𝑡𝑛+1𝑡𝑛𝐮𝑡2.𝑑𝑡(5.18) For the 𝑏(,,) terms of 𝐅(𝑌1𝑛+1,𝑈1𝑛+1), we have 𝜆1|||𝑏𝑢1,𝑛,Ψ1𝑛+1,𝑈1𝑛+1|||𝜆21𝑈1𝑛+12+𝜈2𝑈1𝑛,𝑢2+𝐾2𝑑2Ψ1𝑛+12,𝜆1|||𝑏𝐮1,𝑛𝐮𝑛+1,𝜎𝑛+1,𝑈1𝑛+1|||𝜆21𝑈1𝑛+12+𝜆21𝜈2𝑈1𝑛,𝑢2+1.5𝑑3𝑀2𝑌𝑛12+1.5𝑑3𝑀2Φ𝑛12+1.5𝑑3𝑀2𝑘𝑡𝑛+1𝑡𝑛𝐮2|||𝑑𝑡,𝜆𝜆1𝑏𝐮𝑛+1,𝜎𝑛+1,𝑈1𝑛+1|||12𝜆𝜆12𝑀4+𝑈1𝑛+12+𝜈2𝑈1𝑛,𝑢2.(5.19) For the 𝑔𝑎 terms of 𝐅(𝑌1𝑛+1,𝑈1𝑛+1), we can obtain 𝜆|||𝑔𝑎𝜎1,𝑛,𝐮1,𝑛𝑔𝑎𝜎𝑛+1,𝐮𝑛+1,𝑈1𝑛+1|||8𝑑2𝐾2𝜆2𝜖6𝑈1𝑛+12+𝜖3𝑌𝑛12+8𝑑2𝐾2𝜆2𝜖3𝜈2𝑈1𝑛,𝑢2+8𝑑2𝐾2Φ𝑛12+5𝜆2𝑈1𝑛+12+5𝜆2𝜈2𝐹1𝑛,𝑢2+8𝑑2𝐾2𝜆2𝑘𝑡𝑛+1𝑡𝑛𝐮𝑡2𝑑𝑡+8𝑑2𝑀2𝑈𝑛12+8𝑑2𝑀2Ψ𝑛12+8𝑑2𝑀2𝑘𝑡𝑛+1𝑡𝑛𝜎𝑡2𝑑𝑡.(5.20) Now, we estimate other terms of 𝐅(𝑌1𝑛+1,𝑈1𝑛+1); ||𝑝2𝛼𝑛+1,𝑌1𝑛+1||𝐶2𝛼2𝐾𝜖4𝐷𝑌1𝑛+12+𝑑4𝜖4𝑝𝑛+1̃𝑝𝑛+12,||𝑑2𝛼Re𝑡Φ1𝑛+1,𝑌1𝑛+1||2𝛼Re2𝑌1𝑛+12+14𝑑𝑡Φ1𝑛+12,||𝐷Φ4𝛼(1𝛼)1𝑛+1𝑌,𝐷1𝑛+1||𝜖2𝛼(1𝛼)5𝐷𝑌1𝑛+12+𝐷Φ1𝑛+12𝜖5,||Ψ2𝛼1𝑛+1𝑌,𝐷1𝑛+1||𝜖2𝛼6𝐷𝑌1𝑛+12+14𝜖6Ψ1𝑛+12,𝑑2𝛼Re𝑡𝐮𝑛+1𝐮𝑡𝑛+1,Yn1+12𝛼Re2𝑌1𝑛+12+14𝑑𝑡𝐮𝑛+1𝐮𝑡𝑛+12,𝜆||𝑑𝑡Ψ1𝑛+1,𝑈1𝑛+1||𝜆2𝑈1𝑛+12+14𝑑𝑡Ψ1𝑛+12,|||𝐷Φ2𝛼1𝑛+1,𝑈1𝑛+1|||𝑈1𝑛+12+𝜈2𝑈1𝑛,𝑢2+2𝛼2Φ1𝑛+12,|||Ψ1𝑛+1,𝑈1𝑛+1|||𝑈1𝑛+12+𝜈2𝑈1𝑛,𝑢2+12Ψ1𝑛+12,𝜆||𝑑𝑡𝜎𝑛+1𝜎𝑡𝑛+1,𝑈1𝑛+1||𝜆2𝑈1𝑛+12+14𝑑𝑡𝜎𝑛+1𝜎𝑡𝑛+12,||𝜆𝜎𝑡𝑛+1,𝜈𝑈1𝑛,𝑢||𝜆22+𝑑𝑌𝑛1𝑈1𝑛+12+𝜈24𝑑2𝑀2+𝑑𝑌𝑛1𝜎𝑡𝑛+12+𝜈24𝐾2𝑑𝜎𝑡𝑛+12.(5.21) Combining the estimates in (5.18)–(5.21), we obtain the estimate for 𝐅(𝑌1𝑛+1,𝑈1𝑛+1): 𝐅𝑌1𝑛+1,𝑈1𝑛+1𝐶2𝛼2𝐾𝜖4+(1𝛼)𝜖5+𝜖6𝐷𝑌1𝑛+12+𝐶2𝐾𝜖3𝐷𝑌1𝑛+12+2𝛼3Re2𝑌+1.51𝑛+12+3𝑑𝑀22+𝛼Re2𝑑2𝑀2𝑌𝑛12+𝑑𝑌𝑛1+3+9𝜆2+2𝜆21+8𝑑2𝐾2𝜆2𝜖3𝑈1𝑛+12+8𝑑2𝑀2𝑈𝑛12+𝜈25𝜆2+𝜆21+4+8𝑑2𝐾2𝜆2𝜖3𝑈1𝑛,𝑢2+2𝛼𝑑4𝜖4𝑝𝑛+1̃𝑝𝑛+12+𝛼Re2𝑑2𝑀2+3𝑑3𝑀22Φ𝑛12𝐾+2𝛼2𝑑4+1𝛼𝜖5+2𝛼2Φ1𝑛+12+8𝑑2𝐾2Φ𝑛12+𝛼2𝑑𝑡Φ1𝑛+12+14𝑑𝑡Ψ1𝑛+12+𝛼2𝜖6+12Ψ1𝑛+12+𝑑𝐾22Ψ1𝑛+12+8𝑑2𝑀2Ψ𝑛12+𝛼2𝑑𝑡𝐮𝑛+1𝐮𝑡𝑛+12+14𝑑𝑡𝜎𝑛+1𝜎𝑡𝑛+12+𝜈24𝑑𝑌𝑛1+𝑑2𝑀2𝜎𝑡𝑛+12+𝛼Re2𝑑2𝑀2𝑘+1.5𝑑3𝑀2𝑘×𝑡𝑛+1𝑡𝑛𝐮𝑡2𝜈𝑑𝑡+24𝐾2𝑑𝜎𝑡𝑛+12+8𝑑2𝑀2𝑘𝑡𝑛+1𝑡𝑛𝜎𝑡2𝑑𝑡+8𝑑2𝐾2𝑘𝑡𝑛+1𝑡𝑛𝐮𝑡21𝑑𝑡+2𝜆𝜆12𝑀4.(5.22) Choosing 𝜖1=(1𝛼/12𝐶2𝐾Re3𝛼),𝜖2=(1𝛼/12𝛼),𝜖3=(2𝛼(1𝛼)/6𝐶2𝐾),𝜖4=(1𝛼/6𝐶2𝐾),𝜖5=1/6,𝜖6=(1𝛼)/6, and substituting (5.17) and (5.22) into (5.16) yields 𝑌𝛼Re𝑠12+𝜆2𝑈𝑠12+2𝛼(1𝛼)𝑠1𝑛=0𝑘𝐷𝑌1𝑛+12+𝑘𝑠1𝑛=0𝑈1𝑛+12+𝜈𝜆1𝜈27𝛼+2(1𝛼)28𝑑2𝐶2𝐾𝐾2𝜆2+𝛼(1𝛼)174+5𝜆2+𝜆21𝑘𝑠1𝑛=0𝑈1𝑛,𝑢2𝐶3𝑘𝑠1𝑛=0𝑌1𝑛+12+𝑘𝑠1𝑛=0𝐶41+𝑌𝑛1𝑈1𝑛+12+𝐶5𝑘𝑠1𝑛=0Φ1𝑛+12+𝐶6𝑘𝑠1𝑛=0Φ1𝑛+12+𝛼2𝑘𝑠1𝑛=0𝑑𝑡Φ1𝑛+12+14𝑘𝑠1𝑛=0𝑑𝑡Ψ1𝑛+12+𝐶7𝑘𝑠1𝑛=0Ψ1𝑛+12+𝐶8𝑘𝑠1𝑛=0Ψ1𝑛+12+𝛼2𝑘𝑠1𝑛=0𝑑𝑡𝑢𝑛+1𝑢𝑡𝑛+12+14𝑘𝑠1𝑛=0𝑑𝑡𝜎𝑛+1𝜎𝑡𝑛+12+3𝛼𝑑𝐶2𝐾𝑘1𝛼𝑠1𝑛=0𝑝𝑛+1̃𝑝𝑛+12+𝑘𝑠1𝑛=0𝜈4𝑑2𝑀2+𝑑𝑌𝑛1𝜎𝑡𝑛+12+𝑘2𝛼Re2𝑑2𝑀2+1.5𝑑3𝑀2𝑢𝑡20,0+𝐾2𝜈𝑑24||𝜎𝑡||20,1+8𝑑2𝑘2𝑀2𝜎𝑡20,0+𝐾2𝑢𝑡20,1+𝑘𝑠1𝑛=012𝜆𝜆12𝑀4.(5.23) According to the approximating properties and the definition (3.5), we can get 𝑘𝑠1𝑛=0Φ1𝑛+12+𝑘𝑠1𝑛=0Ψ1𝑛+12𝑘𝐶𝑠1𝑛=04𝑢𝑛+123+2𝜎𝑛+122𝐶4||||𝑢20,3+2||||𝜎20,2,𝑘(5.24)𝑠1𝑛=0Φ1𝑛+12+𝑘𝑠1𝑛=0Ψ1𝑛+12+𝑘𝑠1𝑛=0𝑝𝑛+1̃𝑝𝑛+12𝐶6𝑘𝑠1𝑛=0𝑢𝑛+123+4𝑘𝑠1𝑛=0𝜎𝑛+122+4𝑘𝑠1𝑛=0𝑝𝑛+122𝐶6||||𝑢20,3+4||||𝜎20,2+4||||𝑝20,2,𝑘(5.25)𝑠1𝑛=0𝑑𝑡Φ1𝑛+12𝐶6||𝑢𝑡||20,3,𝑘(5.26)𝑠1𝑛=0𝑑𝑡Ψ1𝑛+12𝐶4||𝑢𝑡||20,2,𝑘(5.27)𝑠1𝑛=0𝑑𝑡𝑢𝑛+1𝑢𝑡𝑛+1213𝑘2𝑢𝑡𝑡20,0,𝑘(5.28)𝑠1𝑛=0𝑑𝑡𝜎𝑛+1𝜎𝑡𝑛+1213𝑘2𝜎𝑡𝑡20,0.(5.29) Combining the inequalities (5.24)–(5.29) with (5.23), we yield 𝑌𝛼Re𝑠12+𝜆2𝑈𝑠12+𝑘𝑠1𝑛=0𝐷𝑌2𝛼(1𝛼)1𝑛+12+𝜈2𝜆1𝑈1𝑛,𝑢2𝐶𝑘𝑠1𝑛=0𝑌1𝑛+12+𝑈1𝑛+12+𝑌𝑛1𝑈1𝑛+12+𝐶𝜈2||𝜎𝑡||20,1+𝜈24𝑘𝑠1𝑛=0𝑌𝑛1𝜎𝑡𝑛+12+𝐶𝑘2𝑢𝑡20,1+𝜎𝑡20,0+𝑢𝑡𝑡20,0+𝜎𝑡𝑡20,0+𝐶6||||𝑢20,3+𝐶4||||𝜎20,2+𝐶4||||𝑝20,2+𝐶4||||𝑢20,3+𝐶6||𝑢𝑡||20,3+𝐶2||||𝜎20,2+𝐶4||𝜎𝑡||20,2+12𝜆𝜆12𝑇𝑀4,(5.30) with 𝜈 chosen such that 𝜆𝜈127𝛼2+(1𝛼)28𝑑2𝐶2𝐾𝐾2𝜆2𝛼+(1𝛼)174+5𝜆2+𝜆211.(5.31) Applying the induction hypothesis (𝐈𝐇𝟐) and the discrete Gronwall's lemma [29] to (5.30), we have 𝑌𝑠12+𝑈𝑠12𝐻(𝑘,,𝜈),(5.32) where 𝐻(𝑘,,𝜈)=𝐶4||||𝑢20,3+6||𝑢𝑡||20,3+2||||𝜎