A Defect-Correction Method for Time-Dependent Viscoelastic Fluid Flow Based on SUPG Formulation

Abstract

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 [10–15] and the references therein) over the years. In [10–12], 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, [16–22].

Inspired by [2, 3, 10–12, 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, 𝑃₁ continuous, 𝑃₂ continuous, and 𝑃₁ 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. 𝑅² 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:𝜆R𝜎e_𝑡𝐮_𝑡+𝑢⋅∇𝑢−∇⋅𝜎−2(1−𝛼)∇⋅𝐷(𝐮)+∇𝑝=𝐟,inΩ×(0,𝑇),𝐮𝐮di(=v0,+𝐮𝟎𝑥,=𝜎)=+0o,n𝐮𝜆₀(Γi𝐮n×⋅Ω∇(0×),𝜎𝑇(0+),,𝑇𝜆𝑔)_𝑎,(𝜎,∇𝐮)−2𝛼𝐷(𝐮)=0,inΩ×(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 𝐮=(𝑢₁,𝑢₂), 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 𝐿²(Ω) inner product is denoted by (⋅,⋅), and the 𝐿^𝑝(Ω) norm by ‖⋅‖_𝐿^𝑝, with the special cases of 𝐿²(Ω) 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:𝑋𝑄∶∶==𝐻𝐿¹⁰₂₀(Ω)²∶=𝑣∈𝐻¹(Ω)²∶𝐯=𝟎onΓ,𝑆∶=𝜏(=Ω)𝜏=_𝑖𝑗𝑞∈𝐿²(Ω);_Ω𝑞𝑑𝑥=0,∩𝜏=𝜏_𝑖𝑗;𝜏_𝑖𝑗=𝜏_𝑗𝑖;𝜏_𝑖𝑗∈𝐿²(Ω);𝑖,𝑗=1,2𝑉∶=𝐯∈𝑋;_Ω;𝐯⋅∇𝜏∈𝐿²(Ω),∀𝐯∈𝑋,𝑞(∇⋅𝐯)𝑑𝑥=0,∀𝑞∈𝑄.(2.3)

The corresponding weak formulation of problem (2.1) is then obtained: find (𝜎,𝐮,𝑝)∈(𝑆,𝑋,𝑄) such that for all (𝜏,𝐯,𝑞)∈(𝑆,𝑋,𝑄)𝜆R𝜎e𝐮_𝑡+Re𝐮⋅∇𝐮,𝐯+(𝜎,𝐷(𝐯))+2(1−𝛼)(𝐷(𝐮),𝐷(𝐯))+(𝑝,∇⋅𝐯)=(𝐟,𝐯),_𝑡(𝑞,∇,⋅𝜏𝐮)+=(𝜎0.,𝜏)+𝜆((𝐮⋅∇)𝜎,𝜏)+𝜆𝑔_𝑎(𝜎,∇𝐮),𝜏−2𝛼(𝐷(𝐮),𝜏)=𝟎,(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 𝐶³ domain, and global existence (in time) of a unique solution for 𝐮 and 𝜎 under a small data assumption on 𝐟, 𝐟_𝑡, 𝐮₀, 𝜎₀. 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 𝐶₁,𝐶₂ such that 𝐶₁ℎ⩽ℎ_𝐾⩽𝐶₂𝜌_𝐾, where ℎ_𝐾 is the diameter of 𝐾, 𝜌_𝐾 is the diameter of the greatest ball included in 𝐾, and ℎ=max_𝐾∈𝑇^ℎℎ_𝐾. Throughout the paper, the constants 𝐶,𝐶₁,𝐶₂,… 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 (𝐮,𝑝)∶𝑃₂-continuous in velocity, 𝑃₁-continuous in pressure; we consider 𝑃₁-continuous approximation of the stress 𝜎. The corresponding FE spaces are 𝑋𝑆^ℎ_ℎ=𝐯∈𝑋∩𝐶⁰(Ω)²;𝐯_|𝐾∈𝑃₂(𝐾)²,∀𝐾∈𝑇^ℎ,𝑄_ℎ=𝜏∈𝑆∩𝐶⁰(Ω)^2×2;𝜏_|𝐾∈𝑃₁(𝐾)^2×2;∀𝐾∈𝑇^ℎ,𝑉_ℎ=𝑞∈𝑄∩𝐶⁰(Ω);𝑞_|𝐾∈𝑃₁(𝐾);∀𝐾∈𝑇^ℎ,=𝐯∈𝑋^ℎ;(𝑞,∇⋅𝑣)=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: ||‖𝑔‖||_∞,𝑚∶=max_{1⩽𝑛⩽𝑁}‖𝑔^𝑛‖_𝑚,||‖𝑔‖||_0,𝑚∶=^𝑁_𝑛=1𝑘‖𝑔^𝑛‖^2𝑚^1/2,(3.3) and when 𝜓(𝐱,𝑡) is defined on the entire time interval (0,𝑇), we denote ‖𝜓‖_∞,𝑚∶=sup_{0⩽𝑡⩽𝑇}‖𝜑(⋅,𝑡)‖_𝑚,‖𝜓‖^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,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,𝜏−2𝛼𝐷𝐮^1ℎ,𝑛+1,𝜏𝜎^1ℎ,𝑛,∇𝐮^1ℎ,𝑛,𝜏,(4.2) where 𝐸 is chosen such that 𝜆₁=𝜆−𝐸ℎ>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},𝜏−2𝛼𝐷𝐮^{ℎ,𝑛+1𝑖+1},𝜏𝜎^{ℎ,𝑛𝑖+1},∇𝐮^{ℎ,𝑛𝑖+1},𝜏−𝜆−𝜆₁𝑏𝐮^{ℎ,𝑛𝑖+1},𝜎^{𝑖ℎ,𝑛+1},𝜏.(4.4) The initial value approximations are taken to be 𝐮^1ℎ,0=𝐮^ℎ,0𝑖+1=𝐮₀, and 𝜎^1ℎ,0=𝜎^ℎ,0𝑖+1=𝜎₀. 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𝑐(𝐮^𝑛+2𝛼(𝜎,𝐷(𝐯))+𝜎,,𝐮,𝐯)+4𝛼(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‖‖²+𝛼²𝑑′𝐾²𝐶^2𝐾Re²𝜖₁|||𝜎_1ℎ,𝑛+1‖‖𝐮^1ℎ,𝑛+1‖‖²,,𝜈𝜎^1𝑛,𝑢|||≤‖‖𝜎^1ℎ,𝑛+1‖‖²⁰+𝜈²4‖‖𝜎^1𝑛,𝑢|||𝜆₁‖‖²,𝑏𝐮^1ℎ,𝑛,𝜎^1ℎ,𝑛+1,𝜎^1ℎ,𝑛+1|||≤𝜖²‖‖𝜎^1𝑛,𝑢‖‖²+𝜆²¹4𝜖₂𝜆₁‖‖𝜎^1ℎ,𝑛+1‖‖²,|||2𝛼𝑏𝐷𝐮𝐮^1ℎ,𝑛,𝜎^1ℎ,𝑛+1_1ℎ,𝑛+1,𝜈𝜎^1𝑛,𝑢=𝜆¹𝜈‖‖𝜎^1𝑛,𝑢‖‖²,,𝜈𝜎^1𝑛,𝑢|||≤𝜖³‖‖𝐷𝐮^1ℎ,𝑛+1‖‖²+𝛼²𝜈²𝜖₃‖‖𝜎^1𝑛,𝑢‖‖²,(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,2𝜎𝛼Re^1ℎ,𝑛+1,𝐮^1ℎ,𝑛+1,𝜎^1ℎ,𝑛+1𝑘−𝛼²𝑑^𝐾²𝐶^2𝐾Re²𝜖₁‖‖𝐮^1ℎ,𝑛+1‖‖²+𝜆𝑘−4𝜆𝜖²¹₂‖‖𝜎^1ℎ,𝑛+1‖‖²+4𝛼(1−𝛼)−𝜖₁−𝜖₃‖𝐷𝐮^1ℎ,𝑛+1‖²+𝜆₁𝜈−𝜖₂−𝜈²4−𝛼²𝜈²𝜖₃‖𝜎^1𝑛,𝑢‖².(4.9) Choosing 𝜈≤(2𝜆₁(1−𝛼))/(1+3𝛼), 𝑘≤min{2𝜆𝜈/𝜆₁,2(1−𝛼)/𝑑^𝐾²𝐶^2𝐾Re} and 𝜖₁=𝜖₃=𝛼(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 (𝑢₁,𝑝₁,𝜎₁) and one correction step approximation (𝑢₂,𝑝₂,𝜎₂). The main results of this section are presented in the following theorem.

Theorem 5.1. Suppose that (2.1) has a solution (𝐮,𝜎,𝑝)∈𝐶²(0,𝑇;𝐻³)×𝐶²(0,𝑇;𝐻²)×𝐶¹(0,𝑇;𝐻²). 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,𝛼,Ω,𝑇,𝐮,𝑝,𝜎,𝐟,𝜆₁,𝜆), 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 𝐵(𝑘,ℎ,𝜈)=+𝐶𝐶ℎ²ℎ||‖₃𝐮‖||_0,3+ℎ³||‖‖𝐮_𝑡‖‖||_0,3+ℎ||‖𝜎‖||_0,2+ℎ²||‖‖𝜎_𝑡‖‖||_0,2+ℎ²||‖𝑝‖||_0,2||‖𝐮‖||_∞,3+ℎ²||‖𝜎‖||_∞,2+𝐶𝑘‖‖𝐮^𝑡‖‖_0,1+‖‖𝐮_𝑡𝑡‖‖_0,0+‖‖𝜎_𝑡‖‖_0,1+‖‖𝜎_𝑡𝑡‖‖_0,0+𝐶𝜈||‖‖𝜎_𝑡‖‖||_0,1+||‖‖𝜎_𝑡‖‖||_∞,0+𝐶√𝑇𝜆−𝜆₁.(5.3)

Remark 5.2. As mentioned in Algorithm 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 (𝐮,𝜎,𝑝)∈𝐶²(0,𝑇;𝐻³)×𝐶²(0,𝑇;𝐻²)×𝐶¹(0,𝑇;𝐻²). Let (𝐮̃^𝑛,𝑝^𝑛̃) denote the Stokes projection of (𝐮^𝑛,𝑝^𝑛) into (𝑉^ℎ,𝑄^ℎ), and 𝜎^𝑛 a Clement interpolant of 𝜎^𝑛 [26, 28]. The following approximating properties are right.‖𝐮^𝑛−𝐮̃^𝑛‖‖+𝑝_𝑛ℎ‖∇(𝐮^𝑛−𝐮̃^𝑛)‖⩽𝐶ℎ³‖𝐮^𝑛‖₃,‖‖𝜎_𝑛−𝑝^𝑛̃‖⩽𝐶ℎ²‖𝑝^𝑛‖₂,−𝜎^𝑛‖‖+ℎ‖‖∇𝜎^𝑛−𝜎^𝑛‖‖⩽𝐶ℎ²‖𝜎^𝑛‖₂.(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+𝐮2𝛼^𝑡𝑛+1,𝜎𝑣+_𝑛+1𝜆𝜎^𝑡𝑛+1,𝜏+2𝛼𝑐Re𝐮^𝑛+1,𝐮^𝑛+1,𝐯+𝜆𝑏𝐮_𝑛+1,𝐷(𝐯)+4𝛼(1−𝛼)𝐷𝐮^𝑛+1,𝐷(𝐯)=2𝛼𝑓_𝑛+1,𝜎^𝑛+1,𝜏+𝜎^𝑛+1,𝜏−2𝛼𝐷𝐮^𝑛+1,𝜏,𝐯+2𝛼𝑝^𝑛+1,∇⋅𝐯−𝜆𝑔_𝑎𝜎^𝑛+1,∇𝐮^𝑛+1,𝜏.(5.8) Equation (5.8) can also be written as 2𝛼Re+𝑑_𝑡4𝐮𝛼(^𝑛+11,𝐯−𝛼+)𝐷𝜆𝑑_𝑡𝐮𝜎^𝑛+1_𝑛+1,𝜏+2𝛼𝑐Re𝐮^1ℎ,𝑛,𝐮^𝑛+1,𝐯+𝜎_𝑛+1,𝐷(𝐯)+2𝛼𝜎^𝑛+1,𝐷(𝐯)+𝜆₁𝑏𝐮^1ℎ,𝑛,𝜎^𝑛+1,𝜏=2𝛼𝐟,_𝑛+1𝜏−2𝛼𝐷𝑢^𝑛+1,𝜏+2𝛼Re,𝐯𝑑_𝑡+2𝛼𝑝^𝑛+1,∇⋅𝐯−𝜆𝑔_𝑎𝜎^1ℎ,𝑛,∇𝐮^1ℎ,𝑛,𝜏−2𝛼𝑐Re𝐮𝐮^𝑛+1_𝑛+1,𝐯−2𝛼Re𝐮^𝑡𝑛+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+𝑑2_𝑡𝛼𝐮𝜎^1ℎ,𝑛+1,𝐯_1ℎ,𝑛+1+𝜆𝑑_𝑡𝜎^1ℎ,𝑛+1,𝜏+2𝛼𝑐Re𝐮^1ℎ,𝑛,𝐮^1ℎ,𝑛+1,𝐯+𝜆₁,𝐷(𝐯)+4𝛼(1−𝛼)𝐷𝐮^1ℎ,𝑛+1,𝐷(𝐯)=2𝛼𝑏𝑓𝐮^1ℎ,𝑛_𝑛+1,𝜎^1ℎ,𝑛+1,𝜏+𝜎^1ℎ,𝑛+1,𝜏−2𝛼𝐷𝐮^1ℎ,𝑛+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+𝑑_𝑡4𝑒𝛼(^1𝑛+11,𝐯−𝛼+)𝐷𝜆𝑑_𝑡𝜖𝑒^1𝑛+1_1𝑛+1,𝜏+2𝛼𝑐Re𝐮^1ℎ,𝑛,𝑒^1𝑛+1,𝑣+𝜖_1𝑛+1,𝐷(𝐯)+2𝛼𝜖^1𝑛+1,𝐷(𝐯)+𝜆₁𝑏𝐮^1ℎ,𝑛,𝜖^1𝑛+1,𝜏=2𝛼Re,𝑑𝜏_𝑡−2𝛼𝐷𝑒^1𝑛+1,𝜏−𝜆𝜎_𝑡𝑛+1𝐮^𝑛+1−𝐮^𝑡𝑛+1,𝐯+2𝛼Re𝑐𝐮^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+𝑑2_𝑡𝑌𝛼𝑐^1𝑛+1R,𝑌e𝐮^1𝑛+1_1ℎ,𝑛+𝜆𝑑_𝑡𝑈^1𝑛+1,𝑈^1𝑛+1+2𝛼𝑈^1𝑛+1,𝐷𝑌^1𝑛+1+𝑈,𝑌^1𝑛+1,𝑌^1𝑛+1+4𝛼(1−𝛼)𝐷𝑌^1𝑛+1,𝐷𝑌^1𝑛+1_1𝑛+1=𝐅𝑌_1𝑛+1,𝑈^1𝑛+1+𝜆¹𝑏𝐮^1ℎ,𝑛,𝑈^1𝑛+1,𝑈^1𝑛+1−2𝛼𝐷𝑌^1𝑛+1,𝑈^1𝑛+1,𝑈^1𝑛+1,(5.12) where 𝐅𝑌=^1𝑛+12,𝛼𝑈𝑝^1𝑛+1_𝑛+1+2𝛼Re𝑐,∇𝐮⋅𝑌^1𝑛+1_1ℎ,𝑛+2𝛼Re𝑑_𝑡𝐮^𝑛+1−𝐮^𝑡𝑛+1,𝑌^1𝑛+1−𝜆𝜎^𝑛+1_t,𝑈^1𝑛+1+𝜆𝑑_𝑡−𝐮^𝑛+1,𝐮^𝑛+1,𝑌^1𝑛+1+𝜆𝑔^𝑎𝜎^1ℎ,𝑛,∇𝐮^1ℎ,𝑛,𝑈^1𝑛+1−4𝛼(1𝜎−𝛼^𝑛+1,)𝑈𝐷^1𝑛+1Φ−𝜆_1𝑛+1𝑔_𝑎𝜎^𝑛+1,∇𝐮^𝑛+1,𝑈^1𝑛+1−2𝛼𝑐Re𝐮^1ℎ,𝑛,Φ^1𝑛+1,𝑌^1𝑛+1−𝜆𝑏𝐮_𝑛+1,𝐷𝑌^1𝑛+1−2𝛼Ψ^1𝑛+1,𝐷𝑌^1𝑛+1−Ψ^1𝑛+1,𝑈^1𝑛+1+2𝛼𝐷Φ,𝜎^𝑛+1_1𝑛+1,𝑈^1𝑛+1+𝜆¹𝑏𝐮^1ℎ,𝑛,𝜎^𝑛+1,𝑈^1𝑛+1−𝜆¹𝑏𝐮^1ℎ,𝑛,Ψ^1𝑛+1,𝑈^1𝑛+1,𝑈^1𝑛+1−2𝛼Re𝑑^𝑡Φ^1𝑛+1,𝑌^1𝑛+1−𝜆𝑑_𝑡Ψ^1𝑛+1,𝑈^1𝑛+1.(5.13) Using the identity (𝑎−𝑏,𝑎)=(1/2)(‖𝑎‖²−‖𝑏‖²+‖𝑎−𝑏‖²), we obtain 𝛼Re𝑘‖‖𝑌^1𝑛+1‖‖²−‖‖𝑌^𝑛1‖‖²+‖‖𝑌^1𝑛+1−𝑌^𝑛1‖‖²+4𝛼(1−𝛼)‖‖𝐷𝑌^1𝑛+1‖‖²+𝜆2𝑘‖‖𝑈^1𝑛+1‖‖²−‖‖𝑈^𝑛1‖‖²+‖‖𝑈^1𝑛+1−𝑈^𝑛1‖‖²+‖‖𝑈^1𝑛+1‖‖²+𝜈𝜆₁‖‖𝑈^1𝑛,𝑢‖‖²⩽2𝛼Re|||𝑐𝐮^1ℎ,𝑛,𝑌^1𝑛+1,𝑌^1𝑛+1|||+𝜈||𝑈^1𝑛+1,𝑈^1𝑛,𝑢||+𝜆₁|||𝑏𝐮_1ℎ,𝑛+2𝛼𝜈||𝐷𝑌_1𝑛+1,𝑈^1𝑛+1,𝑈^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𝑛+1‖‖²−‖‖𝑌^𝑛1‖‖²+‖‖𝑌^1𝑛+1−𝑌^𝑛1‖‖²+4𝛼(1−𝛼)𝑘^𝑠−1_𝑛=0‖‖𝐷𝑌^1𝑛+1‖‖²+𝜆2^𝑠−1_𝑛=0‖‖U^1𝑛+1‖‖²−‖‖𝑈^𝑛1‖‖²+‖‖𝑈^1𝑛+1−𝑈^𝑛1‖‖²+𝑘^𝑠−1_𝑛=0‖‖𝑈^1𝑛+1‖‖²+𝜈𝜆₁𝑘^𝑠−1_𝑛=0‖‖𝑈^1𝑛,𝑢‖‖²⩽𝑘^𝑠−1_𝑛=02𝛼Re|||𝑐𝐮^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 𝑌⁰¹=0 and 𝑈⁰¹=0, 𝛼Re‖‖𝑌^𝑠1‖‖²+^𝑠−1_𝑛=0‖‖𝑌^1𝑛+1−𝑌^𝑛1‖‖²+𝜆2‖‖𝑈^𝑠1‖‖²+^𝑠−1_𝑛=0‖‖𝑈^1𝑛+1−𝑈^𝑛1‖‖²+𝑘_𝑠−1_𝑛=0‖‖𝑈^1𝑛+1‖‖²+4𝛼(1−𝛼)𝑘^𝑠−1_𝑛=0‖‖𝐷𝑌^1𝑛+1‖‖²+𝜈𝜆₁𝑘^𝑠−1_𝑛=0‖‖𝑈^1𝑛,𝑢‖‖²⩽𝑘^𝑠−1_𝑛=02𝛼Re|||𝑐𝐮^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𝛼Re|||𝑐𝐮^1ℎ,𝑛,𝑌^1𝑛+1,𝑌^1𝑛+1|||⩽4𝛼²Re²𝐶^2𝐾𝜖₁‖‖𝐷𝑌^1𝑛+1‖‖²+𝑑^𝐾²4𝜖₁𝜈||𝑈_1𝑛+1‖‖𝑌^1𝑛+1‖‖²,,𝑈^1𝑛,𝑢||⩽‖‖𝑈^1𝑛+1‖‖²+𝜈²4‖‖𝑈^1𝑛,𝑢|||𝜆₁‖‖²,𝑏𝐮^1ℎ,𝑛,𝑈^1𝑛+1,𝑈^1𝑛+1|||⩽𝜆¹2𝑑^2𝛼𝜈||𝐷𝑌_1𝑛+1‖∇𝑌^𝑛1‖_∞‖𝑈^1𝑛+1‖²+𝑀‖𝑈^1𝑛+1‖²,,𝑈^1𝑛,𝑢||⩽4𝛼²𝜖₂‖‖𝐷𝑌^1𝑛+1‖‖²+4𝜈𝜖²₂‖‖𝑈^1𝑛,𝑢‖‖².(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𝛼Re|||𝑐𝐮^1ℎ,𝑛,Φ^1𝑛+1,𝑌^1𝑛+1|||⩽2𝛼Re²‖‖𝑌^1𝑛+1‖‖²+𝐾²𝑑^4‖‖∇Φ^1𝑛+1‖‖²2𝛼Re|||𝑐𝐮_1ℎ,𝑛.−𝐮^𝑛+1,𝐮^𝑛+1,𝑌^1𝑛+1|||⩽2𝛼Re²𝑑^2𝑀²2‖‖𝑌^𝑛1‖‖²+32𝛼‖‖𝑌^1𝑛+1‖‖²+Re²𝑑′²𝑀²2‖‖Φ^𝑛1‖‖²+Re²2𝑑′²𝑀²𝑘^𝑡𝑛+1_𝑡𝑛‖‖𝐮_𝑡‖‖²𝑑𝑡.(5.18) For the 𝑏(⋅,⋅,⋅) terms of 𝐅(𝑌^1𝑛+1,𝑈^1𝑛+1), we have 𝜆₁|||𝑏𝑢_1ℎ,𝑛,Ψ^1𝑛+1,𝑈^1𝑛+1|||⩽𝜆²¹‖‖𝑈^1𝑛+1‖‖²+𝜈²‖‖𝑈^1𝑛,𝑢‖‖²+𝐾²𝑑′2‖‖∇Ψ^1𝑛+1𝜆₁‖‖²,|||𝑏𝐮_1ℎ,𝑛−𝐮^𝑛+1,𝜎^𝑛+1,𝑈^1𝑛+1|||⩽𝜆²¹‖‖𝑈^1𝑛+1‖‖²+𝜆²¹𝜈²‖‖𝑈^1𝑛,𝑢‖‖²+1.5𝑑′³𝑀²‖‖𝑌^𝑛1‖‖²|||𝜆−𝜆₁+1.5𝑑′³𝑀²‖‖Φ^𝑛1‖‖²+1.5𝑑′³𝑀²𝑘^𝑡𝑛+1_𝑡𝑛‖𝐮‖²𝑑𝑡,𝑏𝐮^𝑛+1,𝜎^𝑛+1,𝑈^1𝑛+1|||⩽12𝜆−𝜆¹²𝑀⁴+‖‖𝑈^1𝑛+1‖‖²+𝜈²‖‖𝑈^1𝑛,𝑢‖‖².(5.19) For the 𝑔_𝑎 terms of 𝐅(𝑌^1𝑛+1,𝑈^1𝑛+1), we can obtain 𝜆|||𝑔_𝑎𝜎^1ℎ,𝑛,∇𝐮^1ℎ,𝑛−𝑔^𝑎𝜎^𝑛+1,∇𝐮^𝑛+1,𝑈^1𝑛+1⩽8𝑑′₂|||𝐾²𝜆²𝜖₆‖‖𝑈^1𝑛+1‖‖²+𝜖₃‖‖∇𝑌^𝑛1‖‖²+8𝑑′²𝐾²𝜆²𝜖₃𝜈²‖‖𝑈^1𝑛,𝑢‖‖²+8𝑑′²𝐾²‖‖∇Φ^𝑛1‖‖²+5𝜆²‖‖𝑈^1𝑛+1‖‖²+5𝜆²𝜈²‖‖𝐹^1𝑛,𝑢‖‖²+8𝑑′²𝐾²𝜆²𝑘^𝑡𝑛+1_𝑡𝑛‖‖∇𝐮_𝑡‖‖²𝑑𝑡+8𝑑′²𝑀²‖‖𝑈^𝑛1‖‖²+8𝑑′²𝑀²‖‖Ψ^𝑛1‖‖²+8𝑑′²𝑀²𝑘^𝑡𝑛+1_𝑡𝑛‖‖𝜎_𝑡‖‖²𝑑𝑡.(5.20) Now, we estimate other terms of 𝐅(𝑌^1𝑛+1,𝑈^1𝑛+1); 2𝛼||𝑝^𝑛+1,∇⋅𝑌^1𝑛+1||⩽2𝛼𝐶^2𝐾𝜖₄‖‖𝐷𝑌^1𝑛+1‖‖²+𝑑′4𝜖₄‖‖𝑝^𝑛+1−𝑝̃^𝑛+1‖‖²2𝛼Re||𝑑_𝑡,Φ^1𝑛+1,𝑌^1𝑛+1||⩽2𝛼Re²‖‖𝑌^1𝑛+1‖‖²+14‖‖𝑑^𝑡4𝛼(1−𝛼)||𝐷Φ_1𝑛+1Φ^1𝑛+1‖‖²,,𝐷𝑌^1𝑛+1||⩽2𝛼(1−𝛼)𝜖₅‖‖𝐷𝑌^1𝑛+1‖‖²+‖‖𝐷Φ^1𝑛+1‖‖²2𝛼||Ψ_1𝑛+1𝜖₅,,𝐷𝑌^1𝑛+1||⩽2𝛼𝜖₆‖‖𝐷𝑌^1𝑛+1‖‖²+14𝜖₆‖‖Ψ^1𝑛+1‖‖²2𝛼Re𝑑_𝑡,𝐮^𝑛+1−𝐮^𝑡𝑛+1,Yⁿ¹⁺¹⩽2𝛼Re²‖‖𝑌^1𝑛+1‖‖²+14‖‖𝑑^𝑡𝜆||𝑑_𝑡𝐮^𝑛+1−𝐮^𝑡𝑛+1‖‖²,Ψ^1𝑛+1,𝑈^1𝑛+1||⩽𝜆²‖‖𝑈^1𝑛+1‖‖²+14‖‖𝑑^𝑡2𝛼|||𝐷Φ_1𝑛+1Ψ^1𝑛+1‖‖²,|||Ψ,_1𝑛+1𝑈^1𝑛+1|||⩽‖‖𝑈^1𝑛+1‖‖²+𝜈²‖‖𝑈^1𝑛,𝑢‖‖²+2𝛼²‖‖∇Φ^1𝑛+1‖‖²,,𝑈^1𝑛+1|||⩽‖‖𝑈_1𝑛+1‖‖²+𝜈²‖‖𝑈^1𝑛,𝑢‖‖²+12‖‖Ψ^1𝑛+1𝜆||𝑑_𝑡‖‖²,𝜎^𝑛+1−𝜎^𝑡𝑛+1,𝑈^1𝑛+1||⩽𝜆²‖‖𝑈^1𝑛+1‖‖²+14‖‖𝑑^𝑡||𝜆𝜎_𝑡𝑛+1𝜎^𝑛+1−𝜎^𝑡𝑛+1‖‖²,,𝜈𝑈^1𝑛,𝑢||⩽𝜆²2+𝑑′‖‖∇𝑌^𝑛1‖‖_∞‖‖𝑈^1𝑛+1‖‖²+𝜈²4𝑑′²𝑀²+𝑑′‖‖∇𝑌^𝑛1‖‖_∞‖‖𝜎^𝑡𝑛+1‖‖²+𝜈²4𝐾²𝑑′‖‖∇𝜎^𝑡𝑛+1‖‖².(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𝐾𝜖₄+(1−𝛼)𝜖₅+𝜖₆‖‖𝐷𝑌^1𝑛+1‖‖²+𝐶^2𝐾𝜖₃‖‖𝐷𝑌^1𝑛+1‖‖²+2𝛼3Re²+1.5‖‖𝑌^1𝑛+1‖‖²+3𝑑′𝑀²2+𝛼Re²𝑑′²𝑀²‖‖𝑌^𝑛1‖‖²+𝑑′‖‖∇𝑌^𝑛1‖‖_∞+3+9𝜆²+2𝜆²¹+8𝑑′²𝐾²𝜆²𝜖₃‖‖𝑈^1𝑛+1‖‖²+8𝑑′²𝑀²‖‖𝑈^𝑛1‖‖²+𝜈²5𝜆²+𝜆²¹+4+8𝑑′²𝐾²𝜆²𝜖₃‖‖𝑈^1𝑛,𝑢‖‖²+2𝛼𝑑′4𝜖₄‖‖𝑝^𝑛+1−𝑝̃^𝑛+1‖‖²+𝛼Re²𝑑′²𝑀²+3𝑑′³𝑀²2‖‖Φ^𝑛1‖‖²+2𝛼𝐾²𝑑′4+1𝜖−₅𝛼+2𝛼²‖‖∇Φ^1𝑛+1‖‖²+8𝑑′²𝐾²‖‖∇Φ^𝑛1‖‖²+𝛼2‖‖𝑑^𝑡Φ^1𝑛+1‖‖²+14‖‖𝑑^𝑡Ψ^1𝑛+1‖‖²+𝛼2𝜖₆+12‖‖Ψ^1𝑛+1‖‖²+𝑑′𝐾²2‖‖∇Ψ^1𝑛+1‖‖²+8𝑑′²𝑀²‖‖Ψ^𝑛1‖‖²+𝛼2‖‖𝑑^𝑡𝐮^𝑛+1−𝐮^𝑡𝑛+1‖‖²+14‖‖𝑑^𝑡𝜎^𝑛+1−𝜎^𝑡𝑛+1‖‖²+𝜈²4𝑑′‖‖∇𝑌^𝑛1×_{𝑡𝑛+1𝑡𝑛}‖‖_∞+𝑑′²𝑀²⋅‖‖𝜎^𝑡𝑛+1‖‖²+𝛼Re²𝑑′²𝑀²𝑘+1.5𝑑′³𝑀²𝑘‖‖𝐮_𝑡‖‖²𝑑𝑡+𝜈²4𝐾²+8𝑑′₂𝑑′‖‖∇𝜎^𝑡𝑛+1‖‖²+8𝑑′²𝑀²𝑘^𝑡𝑛+1_𝑡𝑛‖‖𝜎_𝑡‖‖²𝑑𝑡𝐾²𝑘^𝑡𝑛+1_𝑡𝑛‖‖∇𝐮_𝑡‖‖²𝑑𝑡+12𝜆−𝜆¹²𝑀⁴.(5.22) Choosing 𝜖₁=(1−𝛼/12𝐶^2𝐾Re³𝛼),𝜖₂=(1−𝛼/12𝛼),𝜖₃=(2𝛼(1−𝛼)/6𝐶^2𝐾),𝜖₄=(1−𝛼/6𝐶^2𝐾),𝜖₅=1/6,𝜖₆=(1−𝛼)/6, and substituting (5.17) and (5.22) into (5.16) yields 𝛼Re‖‖𝑌^𝑠1‖‖²+𝜆2‖‖𝑈^𝑠1‖‖²+2𝛼(1−𝛼)^𝑠−1_𝑛=0𝑘‖‖𝐷𝑌^1𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0‖‖𝑈^1𝑛+1‖‖²+𝜈𝜆₁−𝜈²2(17−𝛼𝛼)+28𝑑′²𝐶^2𝐾𝐾²𝜆²𝛼(1−𝛼)+147+5𝜆²+𝜆²¹𝑘^𝑠−1_𝑛=0‖‖𝑈^1𝑛,𝑢‖‖²⩽𝐶₃𝑘^𝑠−1_𝑛=0‖‖𝑌^1𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0𝐶₄1+‖‖∇𝑌^𝑛1‖‖_∞‖‖𝑈^1𝑛+1‖‖²+𝐶₅𝑘^𝑠−1_𝑛=0‖‖Φ^1𝑛+1‖‖²+𝐶₆𝑘^𝑠−1_𝑛=0‖‖∇Φ^1𝑛+1‖‖²+𝛼2𝑘^𝑠−1_𝑛=0‖‖𝑑_𝑡Φ^1𝑛+1‖‖²+14𝑘^𝑠−1_𝑛=0‖‖𝑑_𝑡Ψ^1𝑛+1‖‖²+𝐶₇𝑘^𝑠−1_𝑛=0‖‖Ψ^1𝑛+1‖‖²+𝐶₈𝑘^𝑠−1_𝑛=0‖‖∇Ψ^1𝑛+1‖‖²+𝛼2𝑘^𝑠−1_𝑛=0‖‖𝑑_𝑡𝑢^𝑛+1−𝑢^𝑡𝑛+1‖‖²+14𝑘^𝑠−1_𝑛=0‖‖𝑑_𝑡𝜎^𝑛+1−𝜎^𝑡𝑛+1‖‖²+3𝛼𝑑′𝐶^2𝐾1−𝛼𝑘^𝑠−1_𝑛=0‖‖𝑝^𝑛+1−𝑝̃^𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0𝜈4𝑑′²𝑀²+𝑑′‖‖∇𝑌^𝑛1‖‖_∞‖‖𝜎^𝑡𝑛+1‖‖²+𝑘²𝛼Re²𝑑′²𝑀²+1.5𝑑′³𝑀²‖‖𝑢_𝑡‖‖^20,0+𝐾²𝑑′𝜈²4||‖‖𝜎^𝑡‖‖||^20,1+8𝑑′²𝑘²𝑀²‖‖𝜎_𝑡‖‖^20,0+𝐾²‖‖𝑢_𝑡‖‖^20,1+𝑘^𝑠−1_𝑛=012𝜆−𝜆¹²𝑀⁴.(5.23) According to the approximating properties and the definition (3.5), we can get 𝑘^𝑠−1_𝑛=0‖‖∇Φ^1𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0‖‖∇Ψ^1𝑛+1‖‖²⩽𝐶𝑘^𝑠−1_𝑛=0ℎ⁴‖‖𝑢^𝑛+1‖‖²³+ℎ²‖‖𝜎^𝑛+1‖‖²²⩽𝐶ℎ₄||‖𝑢‖||^20,3+ℎ²||‖𝜎‖||^20,2𝑘,(5.24)_𝑠−1_𝑛=0‖‖Φ^1𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0‖‖Ψ^1𝑛+1‖‖²+𝑘^𝑠−1_𝑛=0‖‖𝑝^𝑛+1−𝑝̃^𝑛+1‖‖²⩽𝐶ℎ⁶𝑘^𝑠−1_𝑛=0‖‖𝑢^𝑛+1‖‖²³+ℎ⁴𝑘^𝑠−1_𝑛=0‖‖𝜎^𝑛+1‖‖²²+ℎ⁴𝑘^𝑠−1_𝑛=0‖‖𝑝^𝑛+1‖‖²²⩽𝐶ℎ₆||‖𝑢‖||^20,3+ℎ⁴||‖𝜎‖||^20,2+ℎ⁴||‖𝑝‖||^20,2𝑘,(5.25)_𝑠−1_𝑛=0𝑘‖‖𝑑_𝑡Φ^1𝑛+1‖‖²⩽𝐶ℎ⁶||‖‖𝑢_𝑡‖‖||^20,3,(5.26)_𝑠−1_𝑛=0𝑘‖‖𝑑_𝑡Ψ^1𝑛+1‖‖²⩽𝐶ℎ⁴||‖‖𝑢_𝑡‖‖||^20,2,(5.27)_𝑠−1_𝑛=0‖‖𝑑_𝑡𝑢^𝑛+1−𝑢^𝑡𝑛+1‖‖²⩽13𝑘²𝑘‖𝑢_𝑡𝑡‖^20,0,(5.28)_𝑠−1_𝑛=0‖‖𝑑_𝑡𝜎^𝑛+1−𝜎^𝑡𝑛+1‖‖²⩽13𝑘²‖𝜎_𝑡𝑡‖^20,0.(5.29) Combining the inequalities (5.24)–(5.29) with (5.23), we yield 𝛼Re‖‖𝑌^𝑠1‖‖²+𝜆2‖‖𝑈^𝑠1‖‖²+𝑘^𝑠−1_𝑛=02𝛼(1−𝛼)‖‖𝐷𝑌^1𝑛+1‖‖²+𝜈2𝜆¹‖‖𝑈^1𝑛,𝑢‖‖²⩽𝐶𝑘_𝑠−1_𝑛=0‖‖𝑌^1𝑛+1‖‖²+‖‖𝑈^1𝑛+1‖‖²+‖‖∇𝑌^𝑛1‖‖_∞‖‖𝑈^1𝑛+1‖‖²+𝐶𝜈²||‖‖𝜎_𝑡‖‖||^20,1+𝜈²4𝑘^𝑠−1_𝑛=0+𝐶ℎ₆‖‖∇𝑌^𝑛1‖‖_∞‖‖𝜎^𝑡𝑛+1‖‖²+𝐶𝑘²‖‖𝑢^𝑡‖‖^20,1+‖‖𝜎_𝑡‖‖^20,0+‖‖𝑢_𝑡𝑡‖‖^20,0+‖‖𝜎_𝑡𝑡‖‖^20,0||‖𝑢‖||^20,3+𝐶ℎ⁴||‖𝜎‖||^20,2+𝐶ℎ⁴||‖𝑝‖||^20,2+𝐶ℎ⁴||‖𝑢‖||^20,3+𝐶ℎ⁶||‖‖𝑢_𝑡‖‖||^20,3+𝐶ℎ²||‖𝜎‖||^20,2+𝐶ℎ⁴||‖‖𝜎_𝑡‖‖||^20,2+12𝜆−𝜆¹²𝑇𝑀⁴,(5.30) with 𝜈 chosen such that 𝜈⩽𝜆¹22(17−𝛼𝛼)+28𝑑′²𝐶^2𝐾𝐾²𝜆²𝛼(1−𝛼)+147+5𝜆²+𝜆²¹⁻¹.(5.31) Applying the induction hypothesis (𝐈𝐇𝟐) and the discrete Gronwall's lemma [29] to (5.30), we have ‖‖𝑌^𝑠1‖‖²+‖‖𝑈^𝑠1‖‖²⩽𝐻(𝑘,ℎ,𝜈),(5.32) where 𝐻(𝑘,ℎ,𝜈)=+𝐶𝐶𝑘ℎ⁴₂||‖𝑢‖||^20,3+ℎ⁶||‖‖𝑢_𝑡‖‖||^20,3+ℎ²||‖𝜎‖||^20,2+ℎ⁴||‖‖𝜎_𝑡‖‖||^20,2+ℎ⁴||‖𝑝‖||^20,2+𝐶𝜈₂‖‖𝑢^𝑡‖‖^20,1+‖‖𝑢_𝑡𝑡‖‖^20,0+‖‖𝜎_𝑡‖‖^20,1+‖‖𝜎_𝑡𝑡‖‖^20,0||‖‖𝜎^𝑡‖‖||^20,1+||‖‖𝜎_𝑡‖‖||^2∞,0+𝐶𝜆−𝜆¹²𝑇.(5.33)
That is to say, under the induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐), we establish that the inequality (5.32) is right for all 1⩽𝑠⩽𝑛, and consequently for all 𝑛∶0⩽𝑛⩽𝑇/𝑘, the inequality (5.32) is right. In the next subsection, we will prove that the induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐) are right.
Let us continue to prove the inequality (5.1).
Using triangle inequality, the estimate (5.32), approximation properties (5.5), and (5.24)-(5.25), we have ||‖‖𝐮−𝐮^ℎ1‖‖||^2∞,0+||‖‖𝜎−𝜎^ℎ1‖‖||^2∞,0⩽||‖‖𝑌₁‖‖||^2∞,0+||‖‖Φ₁‖‖||^2∞,0+||‖‖𝑈₁‖‖||^2∞,0+||‖‖Ψ₁‖‖||^2∞,0⩽𝐻(𝑘,ℎ,𝜈)+𝐶ℎ⁶||‖𝑢‖||^2∞,3+ℎ⁴||‖𝜎‖||^2∞,2.(5.34) Similarly, by means of triangle inequality, the estimate (5.32) and (5.30), approximation properties (5.5) and (5.24)-(5.25), we can obtain ||‖‖𝐮−𝐮^ℎ1‖‖||^20,1+||‖‖𝜎−𝜎^ℎ1‖‖||^20,0⩽||‖‖𝑌₁‖‖||^20,1+||‖‖Φ₁‖‖||^20,1+||‖‖𝑈₁‖‖||^20,0+||‖‖Ψ₁‖‖||^20,0⩽𝐶(𝑇)𝐻(𝑘,ℎ,𝜈)+ℎ⁶||‖𝐮‖||^20,3+ℎ⁴||‖𝜎‖||^20,2.(5.35) Thus, we complete the proof of Substep 1.1. That is the inequality (5.1) is right if the induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐) are right.
Substep 1.2. In this subsection, we will prove that the two induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐) for 𝑗=1 are right. We first verify the induction hypothesis (𝐈𝐇𝟏).
Assume (𝐈𝐇𝟏) is right for any 𝑛=0,1,…,𝑠−1, we will prove that (𝐈𝐇𝟏) is right for 𝑛=𝑠. By approximating properties (5.5), inverse inequality [26], and (5.32), we have ‖‖𝑢^1ℎ,𝑠‖‖_∞⩽⩽𝐶‖‖𝑢ℎ^1ℎ,𝑠−_−𝑑′/2𝑢^𝑠‖‖_∞+‖𝑢^𝑠‖_∞⩽‖‖𝑌^𝑠1‖‖_∞+‖‖Φ^𝑠1‖‖_∞+𝑀⩽𝐶𝑘ℎ‖‖𝑌^𝑠1_−𝑑′/2‖‖+𝐶ℎ^−𝑑′/2‖‖Φ^𝑠1‖‖+𝑀+𝜈ℎ^−𝑑′/2+ℎ^{(1−𝑑′)/2}+ℎ^{(2−𝑑′)/2}+ℎ^{(3−𝑑′)/2}+ℎ^−𝑑′/2𝜆−𝜆₁+𝑀.(5.36) As 𝐶(𝑘ℎ^−𝑑′/2+𝜈ℎ^−𝑑′/2+ℎ^{(1−𝑑′)/2}+ℎ^{(2−𝑑′)/2}+ℎ^{(3−𝑑′)/2}+ℎ^−𝑑′/2(𝜆−𝜆₁)) independent of 𝑠, we can choose 𝑘,𝜈⩽ℎ^𝑑′/2/𝐶,𝜆−𝜆₁⩽ℎ^𝑑′/2/𝐶, and get ‖𝑢^1ℎ,𝑠+1‖_∞⩽𝑀+6.
We can get ‖𝜎^1ℎ,𝑠+1‖_∞⩽𝑀+6 by the same method. Define 𝐾=𝑀+6, we confirm (𝐈𝐇𝟏) is right for 𝑛=𝑠.
Now we establish that the induction hypothesis (𝐈𝐇𝟐) is right.
Assume (𝐈𝐇𝟐) is right for any 𝑛=0,1,…,𝑠−1, we will prove that (𝐈𝐇𝟐) is right for 𝑛=𝑠. Using inverse inequality [26], Korn's inequality, and (5.30), we have ^𝑠_𝑛=0𝑘‖∇𝑌^𝑛1‖_∞=^𝑠−1_𝑛=0⩽𝐶ℎ𝑘‖‖∇Y^1𝑛+1_{−𝑑′/2𝑠−1}‖‖_∞as𝑌⁰¹=0_𝑛=0⩽𝐶ℎ_−𝑑′/2𝑘‖‖∇𝑌^1𝑛+1‖‖√𝑘𝑠^𝑠−1_𝑛=0𝑘‖‖∇𝑌^1𝑛+1‖‖²^1/2⩽𝐶√𝑇𝑘ℎ^−𝑑′/2+𝜈ℎ^−𝑑′/2+ℎ^{(1−𝑑′)/2}+ℎ^{(2−𝑑′)/2}+ℎ^−𝑑′/2𝜆−𝜆₁⩽1,(5.37) with the choosing 𝑘,𝜈,𝜆−𝜆₁⩽ℎ^𝑑′/2/6𝐶√𝑇 and ℎ^{(2−𝑑′)/2},ℎ^{(1−𝑑′)/2},ℎ^{(2−𝑑′)/2}⩽1/6𝐶√𝑇.
Step 2. We will show that the inequality (5.2) is true.
In order to get the inequality (5.2), we also need induction hypotheses (𝐈𝐇𝟏) and (𝐈𝐇𝟐) for 𝑗=2. As the procedure of proof is almost same as Step 1, we only give the different places with Step 1.
Now, we combine the correction problem (4.3)-(4.4) with (5.8) and introduce the approximation error 𝑒^2𝑛+1,𝜖^2𝑛+1. This gives 2𝛼Re+𝑑4_𝑡𝑒𝛼(1^2𝑛+1,−𝐯𝛼+)𝐷𝜆𝑑_𝑡𝑒𝜖^2𝑛+1_2𝑛+1,𝜏+2𝛼𝑐Re𝐮^2ℎ,𝑛,𝑒^2𝑛+1,𝐯+2𝛼𝜖^2𝑛+1,𝐷(𝐯)=2𝛼𝑝_𝑛+1,𝐷(𝐯)+𝜆₁𝑏𝐮^2ℎ,𝑛,𝜖^2𝑛+1,𝜏+𝜖^2𝑛+1,𝜏−2𝛼𝐷𝑒^2𝑛+1,𝜏−𝜆𝜎_𝑡𝑛+1,∇⋅𝐯+2𝛼Re𝑑_𝑡𝑢^𝑛+1−𝑢^𝑡𝑛+1,𝐯+2𝛼Re𝑐𝐮^2ℎ,𝑛−𝐮^𝑛+1,𝑢^𝑛+1,𝐯+𝜆−𝜆₁,𝜏+𝜆𝑑_𝑡𝜎^𝑛+1,𝜏+𝜆−𝜆₁𝑏𝐮^2ℎ,𝑛,𝜎^1ℎ,𝑛+1−𝜎^𝑛+1,𝜏𝑏𝐮^𝑛+1,𝜎^𝑛+1,𝜏+𝜆𝑔_𝑎𝜎^2ℎ,𝑛,∇𝐮^2ℎ,𝑛,𝜏−𝜆𝑔_𝑎𝜎^𝑛+1,∇𝐮^𝑛+1,𝜏.(5.38) Substituting 𝑒^2𝑛+1=𝑌^2𝑛+1+Φ^2𝑛+1,𝜖^2𝑛+1=𝑈^2𝑛+1+Ψ^2𝑛+1, 𝑣=𝑌^2𝑛+1 and 𝜏=𝑈^2𝑛+1 into (5.38), we obtain 𝛼Re𝑘‖‖𝑌^2𝑛+1‖‖²−‖‖𝑌^𝑛2‖‖²+2𝜆𝑘‖‖𝑈^2𝑛+1‖‖²−‖‖𝑈^𝑛2‖‖²+4𝛼(1−𝛼)𝐷𝑌_2𝑛+1+𝜆₁,𝐷𝑌^2𝑛+1+2𝛼𝑈^2𝑛+1,𝐷𝑌^2𝑛+1+𝑈^2𝑛+1,𝑈^2𝑛+1⩽2𝛼𝑏𝑝𝐮^2ℎ,𝑛_𝑛+1,𝑈^2𝑛+1,𝑈^2𝑛+1−2𝛼𝐷𝑌^2𝑛+1,𝑈^2𝑛+1+2𝛼𝑐Re𝐮^2ℎ,𝑛,𝑌^2𝑛+1,𝑌^2𝑛+1+2𝛼Re𝑐,∇𝐮⋅𝑌^2𝑛+1_2ℎ,𝑛+2𝛼Re𝑑_𝑡𝐮^𝑛+1−𝐮^𝑡𝑛+1,Y^2𝑛+1+𝜆𝑑_𝑡𝜎^𝑛+1,𝑈^2𝑛+1−𝜆𝜎_𝑡𝑛+1−𝐮^𝑛+1,𝐮^𝑛+1,𝑌^2𝑛+1−𝜆¹𝑏𝐮^2ℎ,𝑛,Ψ^2𝑛+1,𝑈^2𝑛+1−2𝛼𝑐R,e𝑈𝐮^2𝑛+1_2ℎ,𝑛+𝜆𝑔^𝑎𝜎^2ℎ,𝑛,∇𝐮^2ℎ,𝑛,𝑈^2𝑛+1−𝜆𝑔^𝑎𝜎^𝑛+1,∇𝐮^𝑛+1,𝑈^2𝑛+1−2𝛼Ψ_2𝑛+1,Φ^2𝑛+1,𝑌^2𝑛+1−4𝛼(1−𝛼)𝐷Φ^2𝑛+1,𝐷𝑌^2𝑛+1+𝜆−𝜆₁,𝐷𝑌^2𝑛+1−𝜆−𝜆₁𝑏𝐮^2ℎ,𝑛,Ψ^1𝑛+1,𝑈^2𝑛+1−2𝛼Re𝑑_𝑡𝑏𝐮^𝑛+1,𝜎^𝑛+1,𝑈^2𝑛+1+2𝛼𝐷Φ^2𝑛+1,𝑈^2𝑛+1−Ψ^2𝑛+1,𝑈^2𝑛+1Φ^2𝑛+1,𝑌^2𝑛+1−𝜆𝑑_𝑡Ψ^2𝑛+1,𝑈^2𝑛+1−𝜆−𝜆₁𝑏𝐮^2ℎ,𝑛,𝑈^1𝑛+1,𝑈^2𝑛+1.(5.39) Comparing (5.39) with (5.12), we find that they are only different on the following terms: −𝜆−𝜆₁𝑏𝐮^2ℎ,𝑛,Ψ^1𝑛+1,𝑈^2𝑛+1−𝜆−𝜆¹𝑏𝐮^2ℎ,𝑛,𝑈^1𝑛+1,𝑈^2𝑛+1.(5.40) So we will only deal with these terms. We have |||𝑏𝐮_2ℎ,𝑛,Ψ^1𝑛+1,𝑈^2𝑛+1|||⩽𝑑^𝑘²2‖‖∇Ψ^1𝑛+1‖‖²+‖‖𝑈^2𝑛+1‖‖²+𝑣²‖‖𝑈^𝑛,u₂‖‖²⩽𝑑^𝑘²2𝑐ℎ²|||𝑏𝐮_2ℎ,𝑛‖‖𝜎^𝑛+1‖‖²²+‖‖𝑈^2𝑛+1‖‖²+𝑣²‖‖𝑈^𝑛,u₂‖‖²,,U^1𝑛+1,𝑈^2𝑛+1|||⩽𝑑^𝑘²2‖‖∇𝑈^1𝑛+1‖‖²+‖‖𝑈^2𝑛+1‖‖²+𝑣²‖‖𝑈^𝑛,u₂‖‖²⩽𝑑^𝑘²𝑐ℎ⁻²2+‖‖𝑈^1𝑛+1‖‖²+‖‖𝑈^2𝑛+1‖‖²+𝑣²‖‖𝑈^𝑛,u₂‖‖²,(5.41)
To conclude, repeat the proof of the first statement of Theorem 5.1, replacing 𝐮^1ℎ,𝑛, 𝜎^1ℎ,𝑛, 𝑌^𝑛1, Φ^𝑛1, 𝑈^𝑛1 and Ψ^𝑛1 by 𝐮^2ℎ,𝑛, 𝜎^2ℎ,𝑛, 𝑌^𝑛2, Φ^𝑛2, 𝑈^𝑛2, and Ψ^𝑛2, respectively, using Korn inequality, approximation properties (5.5), (5.24)–(5.29), and the bound ‖𝑌^1𝑛+1‖²+‖𝑈^1𝑛+1‖²⩽𝐻(𝑘,ℎ,𝜈). Hence, we can obtain ||‖‖𝐮−𝐮^ℎ2‖‖||_∞,0+||‖‖𝜎−𝜎^ℎ2‖‖||_∞,0+||‖‖𝐮−𝐮^ℎ2‖‖||_0,1+||‖‖𝜎−𝜎^ℎ2‖‖||_0,0⩽𝐵(𝑘,ℎ,𝜈).(5.42)

6. Numerical Results

In this section, numerical results for the defect correction method applied to viscoelastic fluid flow are presented using two test problems. The first example is a known analytical solution to verify numerical convergence rates for the defect correction method. The second example simulates viscoelastic flow through a four-to-one contraction flow, a prototypical problem for viscoelastic fluid flow. As mentioned above, continuous piecewise quadratic elements were used for modeling the velocity, and continuous piecewise linear elements were used for the pressure and stress. The constitutive equation was stabilized using an SUPG discretization with parameter 𝜈. In this paper, we will not investigate the influence of the parameter 𝜈, thus, we set 𝜈=0.6 h.

The defect correction algorithms are implemented using public domain finite element software [30]. Linear systems are solved using the UMFPACK solver. We use the stopping criterion defined by max{‖𝐮^ℎ𝑖−𝐮^ℎ𝑖−1‖,‖𝜎^ℎ𝑖−𝜎^ℎ𝑖−1‖}⩽10⁻⁸ for the iterative solver in both the defect step and the correction step of the method. We also set the maximum number of iteration equal to 15.

Example 6.1 (Analytical solution). The theoretical convergence rates were verified by considering fluid flow across a unit square with a known solution. As in [7, 31], Ω=[0,1]², and chosen functions are added to the right-hand sides of (2.1) such that the true solution to the problem is given by 𝑢=10𝑥²(𝑥−1)²𝑦(𝑦−1)(2𝑦−1)𝑒^−𝑡−1𝑝0=𝑥(1𝑥0−(21𝑥)−(21𝑥)−(21𝑦)−𝑦²1(𝑦)𝑒−_−𝑡1)²𝑒^−𝑡,𝜎=2𝛼𝐷(𝐮).,(6.1)
Let 𝑟 be the experimental global rate of convergence given by 𝑟=log(𝐸𝑟/𝐸𝑟′)/log(ℎ/ℎ′), where ℎ and ℎ′ denote two consecutive mesh sizes with corresponding global errors 𝐸𝑟 and 𝐸𝑟′. In this example, we select Re=1, 𝑎=0, 𝛼=0.5. The numerical results for Example 6.1 are presented in Tables 1, 2, 3, and 4.

Table 1 

Defect step approximation for 𝜆=5.0.

Table 2 

Correction step approximation for 𝜆=5.0.

Table 3 

Defect step approximation for 𝜆=1000.

Table 4 

Correction step approximation for 𝜆=1000.

To reduce the influence of the time discretization error, the time step is taken to be very small: 𝑘=𝑂(ℎ²).

For 𝜆=5, 𝜆=1000 and the final time 𝑇=0.1, the calculated convergence rates in Tables 1–4 confirm what is predicted by Theorem 5.1 for continuous (𝑃₂,𝑃₁,𝑃₁) elements in space. In fact, our numerical convergence rates are better than the theoretical ones. We will find the reason in future work.

Example 6.2 (Four-to-one contraction flow). Numerical simulations of viscoelastic flow through a planar or axisymmetric contraction have been widely studied [32, 33]. Here the case of planar flow through a contraction geometry with a ratio of 4 : 1 with respect to upstream and downstream channel widths is considered. The contraction angle is a fixed 3𝜋/2 and the channel lengths are sufficiently long to impose a fully developed Poiseuille flow in the inflow and outflow channels. The geometry of the computational domain is illustrated in Figure 1. The lower left corner of the domain corresponds to 𝑥=𝑦=0. The computations of the mesh are also shown in Figure 1 with Δ𝑥_min=0.0625 and Δ𝑦_min=0.015625. We denote Γ_in={(𝑥,𝑦)∶𝑥=0,0⩽𝑦⩽1} and Γ_out={(𝑥,𝑦)∶𝑥=8,0⩽𝑦⩽0.25}. On this domain the velocity boundary conditions are 𝑢₁=(1/32)(1−𝑦²),𝑢₂=0,onΓ_in and 𝑢₁=2((1/16)−𝑦²),𝑢₂=0,onΓ_out. On Γ_in, specified boundary conditions for 𝜎 are given as follows: 𝜎₁₁=(−𝛼𝜆(𝑎+1)(−𝑦/16)²)/((𝑎²−1)𝜆²(−𝑦/16)²−1), 𝜎₁₂=𝜎₂₁=−𝛼(−𝑦/16)/((𝑎²−1)𝜆²(−𝑦/16)²−1), 𝜎₂₂=−𝛼𝜆(𝑎−1)(−𝑦/16)²/((𝑎²−1)𝜆²(−𝑦/16)²−1). Symmetry conditions are imposed on the bottom of the computational domain. In this example, the parameters Re, 𝛼, 𝜆, and 𝑎 are set to 1, 8/9, 1.3, and 0, respectively.

Figure 1 

Plot of 4∶1 contraction domain geometry and sample mesh.

We performed the following study: starting from rest, we measured the time that the approximation solution reaches a steady state. The criterion to stop this process is the following: max‖‖𝐮^ℎ𝑛+1−𝐮^𝑛ℎ‖‖_𝐿2(Ω)‖‖𝐮^ℎ𝑛+1‖‖_𝐿2(Ω),‖‖𝜎^ℎ𝑛+1−𝜎^𝑛ℎ‖‖_𝐿2(Ω)‖‖𝜎^ℎ𝑛+1‖‖_𝐿2(Ω)⩽10⁻⁵,(6.2) where 𝑛+1,𝑛 denote 𝑡_𝑛+1,𝑡_𝑛, respectively.

In Figure 2, we plot the evolution of the kinetic energy ‖𝐮^ℎ𝑛+1‖²⁰/2 and ‖𝜎^ℎ𝑛+1‖²⁰/2 using time step 𝑘=0.1 until it reaches its steady state, where we observe they convergence towards a steady state and also the absence of oscillations along the process.

(a)

(b)

Figure 2 

Evolution of ‖𝐮ℎ𝑛+1‖20/2 (a) and ‖𝜎ℎ𝑛+1‖20/2 (b) in time with time step 𝑘=0.1.

Figure 3 presents the horizontal and vertical velocities near the reentrance corner along the vertical line 𝑥=4.0625 for 𝜆=1.3. We observe that the horizontal velocity is almost continuous, while the vertical velocity has high gradients near 𝑦=0.11 and 𝑦=0.23 from Figure 3. However, we find that the solutions of the time-dependent problem can converge to the solutions of the steady problem. Figure 4 present the streamlines of the fluid with 𝜆=1.3. We plot the streamlines for both the steady problem and the time-dependent problem at final time 𝑡=24.5. It is easy to observe that these two figures are almost alike.

(a)

(b)

Figure 3 

Horizontal (a) and vertical (b) velocity near reentrant corner. The marks “+” indicate results for steady problem and “–” indicate results for the time-dependent problem at time 𝑡=24.5 using time step 𝑘=0.1.

(a)

(b)

Figure 4 

Streamlines and magnitude of velocity contours for 𝐮. (a): steady problem, (b): time-dependent problem at time 𝑡=24.5 using 𝑘=0.1.

7. Conclusions

In this paper, we present a defect correction mixed finite element method for solving the time-dependent Johnson-Segalman viscoelastic equations. A priori error estimates for the defect step and the first correction step of the defect correction method are derived. Finally, we present two numerical examples. One is a known problem and the other is a benchmark problem.

For simplicity, in this paper, we choose the parameter 𝜈=0.6 h. More appropriate choice of parameter 𝜈 is currently under investigation. Further developments should be extending the method to other non-Newtonian flow or other discrete scheme.

Acknowledgments

The authors would like to thank the editor and the anonymous referees for their criticism and valuable comments, which led to the improvement of this paper. The authors want to thank Professor J. S. Howell for meaty the discussion of part program code. This paper was supported by the National Natural Science Foundation of China (no. 10871156, 10861014) and the Fund of Xi'an Jiaotong University (no. 2009xjtujc30).